L Najman
email: l.najman@esiee.fr

M Couprie
email: m.couprie@esiee.fr

Building the component tree in quasi-linear time

Keywords: Component tree, connected operators, mathematical morphology, classification, disjoint sets, union-find, image and signal processing, filtering

The level sets of a map are the sets of points with level above a given threshold. The connected components of the level sets, thanks to the inclusion relation, can be organized in a tree structure, that is called the component tree. This tree, under several variations, has been used in numerous applications. Various algorithms have been proposed in the literature for computing the component tree. The fastest ones (considering the worst-case complexity) have been proved to run in O(n ln(n)).

In this paper, we propose a simple to implement quasi-linear algorithm for computing the component tree on symmetric graphs, based on Tarjan's union-find procedure. We also propose an algorithm that computes the n most significant lobes of a map.

I. INTRODUCTION

The level sets of a map are the sets of points with level above a given threshold. The connected components of the level sets, thanks to the inclusion relation, can be organized in a tree structure, that is called the component tree. The component tree captures some essential features of the map. It has been used (under several variations) in numerous applications among which we can cite: image filtering and segmentation [START_REF] Hanusse | Sémantique des images par analyse dendronique[END_REF], [START_REF] Guillataud | Contribution à l'analyse dendroniques des images[END_REF], [START_REF] Couprie | Topological grayscale watershed transform[END_REF], [START_REF] Jones | Component trees for image filtering and segmentation[END_REF], video segmentation [START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF], image registration [START_REF] Mattes | Tree representation for image matching and object recognition[END_REF], [START_REF] Monasse | Morphological representation of digital images and application to registration[END_REF], image compression [START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF] and data visualization [START_REF] Chiang | Simple and optimal outputsensitive construction of contour trees using monotone paths[END_REF]. This tree is also fundamental for the efficient computation of the topological watershed introduced by M. Couprie and G. Bertrand [START_REF] Couprie | Topological grayscale watershed transform[END_REF], [START_REF] Couprie | Quasi-linear algorithms for the topological watershed[END_REF], [START_REF] Bertrand | On topological watersheds[END_REF].

While having been (re)discovered by several authors for image processing applications, the component tree concept was first introduced in statistics [START_REF] Wishart | Mode analysis: A generalization of the nearest neighbor which reduces chaining effects[END_REF], [START_REF] Hartigan | Statistical theory in clustering[END_REF] for classification and clustering. For image processing, the use of this tree in order to represent the "meaningful" information contained in a numerical function can be found in particular, in a paper by Hanusse and Guillataud [START_REF] Hanusse | Sémantique des images par analyse dendronique[END_REF], [START_REF] Guillataud | Contribution à l'analyse dendroniques des images[END_REF]; the authors claim that this tree can play a central role in image segmentation, and suggest a way to compute it, based on an immersion simulation. Several authors, such as Vachier [START_REF] Vachier | Extraction de caractéristiques, segmentation d'images et Morphologie Mathématique[END_REF], Breen and Jones [START_REF] Breen | Attribute openings, thinnings and granulometries[END_REF], Salembier et al. [START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF] have used some variations of this structure in order to implement efficiently some morphological operators (e.g. connected operators [START_REF] Salembier | Flat zones filtering, connected operators and filter by reconstruction[END_REF], granulometries, extinction functions, dynamics [START_REF] Bertand | On the dynamics[END_REF]).

Let us describe informally an "emergence" process that will later help us designing an algorithm for building the component tree. Using topographical references, we see the map as the surface of a relief, with the level of a point corresponding to its altitude. Imagine that the surface is completely covered by water, and that the level of water slowly decreases. Islands (regional maxima) appear. These islands form the leafs of the component tree. As the level of water decreases, islands grow, building the branches of the tree. Sometimes, at a given level, several islands merge into one connected piece. Such pieces are the forks of the tree. We stop when all the water has disappeared. The emerged area forms a unique component: the root of the tree.

Various algorithms have been proposed in the literature for computing the component tree [START_REF] Breen | Attribute openings, thinnings and granulometries[END_REF], [START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF], [START_REF] Mattes | Efficient algorithms to implement the confinement tree[END_REF], the latter reference also contains a discussion about time complexity of the different algorithms. The fastest ones (considering the worst-case complexity) have been proved to run in O(n ln(n)), where n denotes the number of pixels of the image. In this paper 1 , we propose a quasi-linear algorithm for computing the component tree of functions defined on general symmetric graphs, based on Tarjan's union-find [START_REF] Tarjan | Efficiency of a good but not linear set union algorithm[END_REF] procedure. More precisely, our algorithm runs in O(N × α(N )) where N denotes the size of the graph (number of vertices + number of edges) and α is a very slow-growing "diagonal inverse" of the Ackermann's function (we have α(10 80 ) ≈ 4). We would like to emphasize that this algorithm is simple to implement.

The paper is organised as follows: we first recall the definitions of some basic graph notions and define the component tree in this framework. We explain the disjoint set problem, together with the solution proposed by Tarjan. Using a disjoint set fomulation, we present our component tree algorithm, and we describe its execution on an example. We then show that the proposed algorithm is quasi-linear with respect to the size of the graph, and compare it to one of the most cited component tree algorithm. We illustrate the use of the component tree for automatic detection of some image features, based on a unique parameter which is the number of features that we expect to find in the image.

II. VERTEX-WEIGHTED GRAPH AND COMPONENT TREE

A. Basic notions for graphs

Let V be a finite set of vertices (or points), and let P(V ) denote the set of all subsets of V . Throughout this paper, E denotes a binary relation on V (that is, a subset of the cartesian product V × V ) which is anti-reflexive ((x, x) / ∈ E) and symmetric ((x, y) ∈ E ⇔ (y, x) ∈ E). We say that the pair (V, E) is a graph, and the elements of E are called edges. We denote by Γ the map from V to P(V ) such that, for all x ∈ V , Γ(x) = {y ∈ V |(x, y) ∈ E}. For any point x, the set Γ(x) is called the neighborhood of x. If y ∈ Γ(x) then we say that y is a neighbor of x and that x and y are adjacent.

Let X ⊆ V . Let x 0 , x n ∈ X. A path from x 0 to x n in X is a sequence π = (x 0 , x 1 , . . . , x n ) of points of X such that x i+1 ∈ Γ(x i ), with i = 0 . . . n -1. Let x, y ∈ X, we say that x and y are linked for X if there exists a path from x to y in X. We say that X is connected if any x and y in X are linked for X. We say that Y ⊆ V is a connected component of X if Y ⊆ X, Y is connected, and Y is maximal for these two properties (i.e., Y = Z whenever Y ⊆ Z ⊆ X and Z is connected).

In the following, we assume that the graph (V, E) is connected, that is, V is made of exactly one connected component.

B. Basic notions for vertex-weighted graphs

We denote by F(V, D), or simply by F, the set composed of all maps from V to D, where D can be any finite set equipped with a total order (e.g., a finite subset of the set of rational numbers or of the set of integers). For a map F ∈ F, the triplet

(V, E, F ) is called a (vertex-)weighted graph. For a point p ∈ V , F (p) is called the weight or level of p. Let F ∈ F, we define F k = {x ∈ V |F (x) ≥ k} with k ∈ D; F k is called a (cross-)section of F . A connected component of a section F k is called a (level k) component of F . A level k component of F that does not contain a level (k + 1
) component of F is called a (regional) maximum of F . We define k min = min {F (x)|x ∈ V } and k max = max {F (x)|x ∈ V }, which represent respectively, the minimum and the maximum level in the map F .

Although the notions we are dealing with in this paper are defined for general graphs, we are going to illustrate our work with the case of 2D images that we model by weighted graphs. Let Z denote the set of integers. We choose for V a subset of Z 2 . A point x ∈ V is defined by its two coordinates (x 1 , x 2 ). We choose for E the 4-connected adjacency relation defined by

E = {(x, y) ∈ V × V ; |x 1 -y 1 | + |x 2 -y 2 | = 1}.
Fig. 1.a shows a weighted graph (V, E, F ) and four crosssections of F , between the level k min = 1 and the level k max = 4. The set F 4 is made of two connected components which are regional maxima of F .

C. Component Tree

From the example of Fig. 1.a, we can see that the connected components of the different cross-sections may be organized, thanks to the inclusion relation, to form a tree structure (see also [START_REF] Bertand | On the dynamics[END_REF]).

Let F ∈ F. For any component c of F , we set h(c) = max{k|c is a level k component of F }. Note that h(c) = min{F (x)|x ∈ c}. We define C(F ) as the set composed of all the pairs [k, c], where c is a component of F and k = h(c). We call altitude of

[k, c] the number k. Remark that [k 1 , c] ∈ C(F ) and [k 2 , c] ∈ C(F ) implies k 1 = k 2 ,
in other words, any two distinct elements of C(F ) correspond to distinct sets of points. Let

1 1 1 1 1 1 1 1 3 3 2 3 4 1 1 3 3 2 3 4 1 1 1 1 1 1 3 1 1 3 3 2 1 1 1 1 4 3 2 2 2 1 1 F 1 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 0 0 0 0 1 0 0 1 1 1 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 F 2 0 0 0 0 0 0 0 0 1 1 0 1 1 0 0 1 1 0 1 1 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 F 3 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 F 4 (a) α1 α2 α3 α4 α5 α6 α8 α7 Level 1 Level 2 Level 3 Level 4 (b) α 1 α 1 α 1 α 1 α 1 α 1 α 1 α 1 α 4 α 4 α 2 α 5 α 7 α 1 α 1 α 4 α 4 α 2 α 5 α 7 α 1 α 1 α 1 α 1 α 1 α 1 α 5 α 1 α 1 α 6 α 6 α 3 α 1 α 1 α 1 α 1 α 8 α 6 α 3 α 3 α 3 α 1 α 1 α 1 α 1 α 1 α 1 α 1 α 1 (c)
F ∈ F, let [k 1 , c 1 ], [k 2 , c 2 ] be distinct elements of C(F ). We say that [k 1 , c 1 ] is the parent of [k 2 , c 2 ] if c 2 ⊂ c 1 and if there is no other [k 3 , c 3 ] in C(F ) such that c 1 ⊂ c 3 ⊂ c 2 .
In this case we also say that

[k 2 , c 2 ] is a child of [k 1 , c 1 ].
With this relation "parent", C(F ) forms a directed tree that we call the component tree of F , and that we will also denote by C(F ) by abuse of terminology. Any element of C(F ) is called a node. An element of C(F ) which has no child (a maximum of F ) is called a leaf, the node which has no parent (i.e., [k min , V ]) is called the root.

We define the component mapping M as the map which associates to each point p ∈ V the node [k, c] of C(F ) such that p ∈ c and F (p) = k. The component mapping is necessary for using the component tree in applications. Fig. 1.b shows the component tree of the weighted graph depicted in Fig. 1.a, and Fig. 1.c shows the associated component mapping. The component at level 1 is called α 1 , the two components at level 2 are called α 2 and α 3 (according to the usual scanning order), and so on.

III. COMPONENT TREE QUASI-LINEAR ALGORITHM

A. Disjoint Sets

The disjoint set problem consists in maintaining a collection Q of disjoint subsets of a set V under the operation of union. Each set X in Q is represented by a unique element of X, called the canonical element. In the following, x and y denote two distinct elements of V . The collection is managed by three operations:

• MakeSet(x): add the set {x} to the collection Q, provided that the element x does not already belongs to a set in Q. • Find(x): return the canonical element of the set in Q which contains x. • Link(x, y): let X and Y be the two sets in Q whose canonical elements are x and y respectively (x and y must be different). Both sets are removed from Q, their union Z = X ∪ Y is added to Q and a canonical element for Z is selected and returned. Tarjan [START_REF] Tarjan | Efficiency of a good but not linear set union algorithm[END_REF] proposed a very simple and very efficient algorithm called union-find to achieve any intermixed sequence of such operations with a quasi-linear complexity. More precisely, if m denotes the number of operations and n denotes the number of elements, the worst-case complexity is O(m × α(m, n)) where α(m, n) is a function which grows very slowly, for all practical purposes α(m, n) is never greater than four 2 .

The implementation of this algorithm is given below in procedure MakeSet and functions Link and Find. Each set of the collection is represented by a rooted tree, where the canonical element of the set is the root of the tree. To each element x is associated a parent Par(x) (which is an element) and a rank Rnk(x) (which is an integer). The mappings 'Par' and 'Rnk' are represented by global arrays in memory. One of the two key heuristics to reduce the complexity is a technique called path compression, that is aimed at reducing, in the long run, the cost of Find. It consists, after finding the root r of the tree which contains x, in considering each element y of the parent path from x to r (including x), and setting the parent of y to be r. The other key technique, called union by rank, consists in always choosing the root with the greatest rank to be the representative of the union while performing the Link operation. If the two canonical elements x and y have the same rank, then one of the elements, say y, is chosen arbitrarily to be the canonical element of the union: y becomes the parent of x; and the rank of y is incremented by one. The rank Rnk(x) is a measure of the depth of the tree rooted in x, and is exactly the depth of this tree if the path compression technique is not used jointly with the union by rank technique. Union by rank avoids creating degenerate trees, and helps keeping the depth of the trees as small as possible. For a more detailed explanation and complexity analysis, see Tarjan's paper [START_REF] Tarjan | Efficiency of a good but not linear set union algorithm[END_REF].

Procedure We can illustrate the use of the union-find algorithm on the classical problem of finding the connected components of a subset X of a graph (V, E). Algorithm 1 (ConnectedComponents) is given below. For a set X, this algorithm returns a map M that gives for each point p, the canonical element M (p) of the connected component of X which contains p.

Algorithm 1: ConnectedComponents

Data: During the first pass (loop 1), for each point p of the set X, the set {p} is added to the collection Q of disjoint subsets. Then, loop 2 processes all points of X in an arbitrary order. For each point p, we first find the canonical element of the set it belongs to (line 3). Then, for each neighbor q of p such that q ∈ X (line 4), we find the canonical element of the set which contains q (line 5). If p and q are not already in the same set, that is if the two canonical elements differ (line 6), then the corresponding sets are merged (line 7), and one of the two canonical elements is chosen to be the canonical element of the merged set. At the end, a simple pass on all the elements of X (loop 8) builds the map M . Note that, if the vertices can be processed in some very specific order (as the scanline order), the ConnectedComponents algorithm becomes linear [START_REF] Ch | Two linear Union-Find strategies for image processing[END_REF], [START_REF] Dillencourt | A general approach to connected-component labeling for arbitrary image representations[END_REF]. Unfortunately, such a specific strategy is not applicable for the component tree algorithm, where the scanning order depends on the altitudes of the vertices.

(V, E) -graph Data: A set X ⊆ V Result: M -map from X to V foreach p ∈ X do MakeSet(p); 1 foreach p ∈ X do 2 compp := Find(p); 3 foreach q ∈ Γ(p) ∩ X do

C. Component tree algorithm: high-level description

We are now ready to introduce our quasi-linear algorithm for building the component tree C(F ) from a weighted graph G = (V, E, F ).

The algorithm simulates the emergence process described in the introduction, and maintains several data structures. The main one is a forest, which initially consists of a set of mutually disconnected nodes, each node being associated (initially) to a single vertex of the graph G. During the emergence process, which is realized by scanning all the vertices of G by decreasing order of altitude, the vertices which belong to a same component and have the same altitude are grouped together thanks to a disjoint set collection called Q node . The canonical element of such a set is called a canonical node. Notice that the disjoint set collection Q node has essentially the same function as the disjoint set collection used by algorithm ConnectedComponents (sec. III-B).

Simultaneously, the canonical nodes are progressively linked together to form partial trees, each partial tree represents intuitively an emerged island. At the end of the execution, a unique tree groups all the canonical nodes, each one of these nodes represents a component of G, and the whole tree constitutes the component tree of G. To reach a quasilinear time complexity, we have to maintain another collection Q tree of disjoint sets, and an auxiliary map called lowestNode. Given an arbitrary node P , the collection Q tree allows to find, in quasi-constant time, a node T which "represents" the partial tree which contains P . Due to the particular management of Q tree , this node T cannot be guaranteed to be precisely the root of the partial tree, this is why we also need to maintain the map lowestNode which associates, to each canonical element of Q tree , the root of the corresponding partial tree.

D. Component tree algorithm: detailed view

Algorithm 2 (BuildComponentTree) is given below. It uses two auxiliary functions MakeNode and MergeNodes. To represent a node of C(F ), we use a structure called node containing the level of the node, and the list of nodes which are children of the current node. For building the component tree, we do not need the reverse link, that is we do not need to know the parent of a given node, but let us note that such information is useful for applications, and can easily be obtained in a linear-time post-processing step. In what follows, we are going to show how to compute some attributes associated to each node of the component tree; we thus need that the structure node contains some fields that store those attributes, namely level, area and highest. We defer both the precise definition of the attributes and the explaination on how they are computed until section VI, in order to concentrate on the component tree itself.

Function node MakeNode(int level)

Allocate a new node n with an empty list of children; n→ level := level; n→ area := 1; n→ highest := level; return n; After a preprocessing (line 1, achievable in linear time for short integers [START_REF] Cormen | Introduction to Algorithms[END_REF]) which sorts the points by decreasing order of level and which prepares the two union-find implementations (line 2), we process the points, starting with the highest ones. Let us suppose that we have processed a number of levels. We have built all nodes of the component tree that are above the current level, and we are building the nodes with exactly the current level. For a given point p of the current level (line 3), we know (through the collection Q tree ) the partial tree the node p belongs to (line 4). In each partial tree, there is only one node with the current level, that we can obtain through the auxiliary map lowestN ode. We then find the associated canonical node (line 5).

We then look at each neighbor q of p with a level greater or equal to the current one (loop 6). Note that, as the graph is symmetric, the "linking operations" between two points are done when one of the two points is processed as a neighbor of the other. Thus, we can use the order of scanning of the points, and we only need to examine the "already processed" neighbors of p. Such a neighbor q satisfies F (q) ≥ F (p).

Exactly as we have done for the point p, we search for the canonical node corresponding to the point q (lines 7-8). If the canonical node of p and the canonical node of q differ, that is if the two points are not already in the same node, we have two possible cases:

• either the two canonical nodes have the same level; this means that these two nodes are in fact part of the same component, and we have to merge the two nodes (line 9 and function MergeNodes). The merging of nodes of same level is done through the collection Q node of disjoint sets. The merging relies on the fact that the Link node function always chooses one of the two canonical elements of the sets that are to be merged as the canonical element of the merged set. This fact is used in the sequel of the function.

Once the merging has been done, one of the nodes is chosen to be the canonical element of the disjoint set.

Observe that the other node is not needed anymore. Indeed, we only have to know to which disjoint set this last node belongs to, and the answer to this question is given by the Find node function.

• or the canonical node of q is strictly above the current level, and thus this node becomes a child of the current Sort the points in decreasing order of level for F ; node (line 10). In both cases, we have to link the two partial trees, this is done using the collection Q tree (line 13). We also have to keep track of the node of lowest level for the union of the two partial trees, that we store in the array lowestN ode (line 14).

At the end of the algorithm, we have to do a post-processing to return the desired result. The root of the component tree can easily be found (line 15) using the array lowestN ode and the two disjoint set structures Q tree and Q node . The component mapping M can be obtained using the disjoint set Q node (loop 16).

IV. ILLUSTRATION OF THE ALGORITHM

Let us illustrate the work of the algorithm on an example. Consider the weighted graph of Fig. 2.a. The points are labelled according to their usual lexicographical order (Fig. 2

.b).

At the beginning of the sixth step, we have already constructed parts of the component tree (Fig. 3.b). We show in Fig. 3.a the maps Par tree , Par node , and lowestN ode. For the maps Par tree and Par node , the canonical elements appear in white. It should be noted that the lowestN ode mapping is only used for the canonical elements of Par tree : this explains why the values of lowestN ode for other elements (in grey) are not updated.

We are going to process nodes at level 50. The first node at level 50 is node 3. Node 0 is a neighbor of node 3. The canonical node corresponding to 0 is node 1, the level of which to the usual lexicographic order, but they will be processed by decreasing level (that is: 12, 0, is 90. Thus node 3 becomes the parent of node 1. Then, node 3 is linked for Q node succesively with nodes 4, 5 and 8. Then node 9 is examined, and is linked for Q node with node 10, the node 9 being chosen as the canonical one. Node 9 is a neighbour of node 12, the canonical element of which is node 13 (level 70). Thus, node 13 becomes a child of node 9. We are then at the beginning of step 11, and this is illustrated on Fig. 4.

Node 11 is a neighbor of both nodes 8 and 10. The canonical node of node 8 is node 3 at level 50. Thus, node 11 and node 3 are linked for Q node , and node 3 is chosen as the canonical one. The canonical node of node 10 is node 9 at level 50. Thus, nodes 9 and 3 are merged, that is, the corresponding partial trees are merged into a single tree. Node 9 is chosen as the canonical element of the level 50 component, and the children of node 3 are transfered to node 9. We are in the situation depicted in Fig. 5.

We then process node 6 at level 40, which becomes the parent of node 9 at level 50. Node 9 and node 6 are linked for Q tree , and node 9 is chosen as the canonical element for the partial tree. The lowest node in this partial tree is node 6 at level 40. We use the map lowestN ode to store that information, by setting lowestN ode[9] := 6. Then we process node 7 at level 20, which becomes the parent of node 6. Node 9 is chosen as the canonical element for the partial tree, and thus we have to store the lowest node by setting lowestN ode[9] := 7. There is no node lower than 20, and thus, the component tree is built. The final situation is depicted in Fig. 6.

The collection Q tree of disjoint sets is not useful anymore: indeed, each node of the graph has been examined, and they are all linked for Q tree , the canonical element being the node 9. The root of the component tree is the node 7. Each of the canonical elements of the collection Q node corresponds to a component of F : observe in particular the level 50, whose canonical node is node 9. The collection Q node can be used to compute the component mapping M .

V. COMPLEXITY ANALYSIS Let n denote the number of points in V , and let m denote the number of edges of the graph (V, E).

The sorting of the points (line 1) can be done in O(n) if the weigths are small integers (counting sort [START_REF] Cormen | Introduction to Algorithms[END_REF]), and in O(n log(log(n))) if each weight can be stored in a machine memory word (long integers or floating point numbers [START_REF] Andersson | Sorting in linear time?[END_REF]).

Loop 2 is the preparation for the union-find algorithm. It is obviously O(n).

In the function MergeNodes, the merging of the lists of children can be done in constant time, because we can merge two lists by setting the first member of one list to be the one that follows the last member of the other list. This requires the two lists to be disjoint, which is the case (we are dealing with disjoint sets), and an adequate representation for lists (chained structure with pointers on both first and last element).

The amortized complexity of line 6 is equal to the number m of edges of the graph (V, E). The amortized complexity of all calls to the union-find procedures is quasi-linear (in the sense explained in section III-A) with respect to m. The building of the component mapping M is obviously linear. Thus the complexity of the algorithm 2 (BuildComponent-Tree) is quasi-linear if the sorting step is linear.

Note that the memory for the lowestN ode array is not necessary: we can easily modify the code so that we store the content of lowestN ode as negative values in Par tree for the canonical element of Q tree . In this case, for an element x ∈ V , Find tree (x) still returns the canonical element c for Q tree , but lowestN ode(c) = -Par tree (c). The modifications that have to be made to MakeSet, to Find, and to BuildComponentTree are straightforward and do not change the complexity of the algorithm.

For comparison purpose, one can prove that the most cited component tree algorithm, the Salembier et al. algorithm [START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF] is quadratic. More precisely, although there is no complexity analysis in [START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF], one can verify that the Salembier et al. algorithm has a worst-case time complexity in O(n × h + m) where h is the number of levels of the image. The worst case can be attained using a series of artificially generated images such that half of the pixels are maxima of the images (an example of an image of the series is provided in Fig. 7). However, this worst case is rare in practice. We observe that, when the level of a point is a short integer (between 0 and 255), the Salembier et al. algorithm is generally twice as fast as our algorithm. This can be explained by the fact that, for each point of the image, we have to access the two union-find data structures, while this is not the case for the Salembier et al. algorithm.

VI. ATTRIBUTES

A major use of the component tree is for image filtering: for example, we may want to remove from an image the "lobes" that are not "important enough" or "negligible". Such an operation is easy to do by simply removing the "negligible" components of the component tree. To make such an idea practicable, it is necessary to quantify the relative importance of each node of the component tree. We can do that by computing some attributes for each node.

Among the numerous attributes that can be computed, three are natural: the height, the area, and the volume (Fig. 8).

Let The area is easy to compute while building the component tree. Each time two components merge (i.e. in the function MergeNodes) or each time a component is declared the parent of another one (i.e. line 11 of algorithm 2 BuildComponent-Tree), we keep as the new area the sum of the areas of the two components.

For computing the highest level in the component, we do as we did for the area, replacing the sum by the maximum (see line 12 of algorithmm 2 BuildComponentTree and the function MergeNodes). From this highest level, the height of a component n can easily be computed by setting height(n) = (n→highest) -(n→level)+1.

To compute the volume, we first need the area. We then apply the recursive function ComputeVolume on the root of the tree. The complexity of this function is linear with respect to the number of nodes.

VII. EXAMPLE OF APPLICATION AND CONCLUSION

We have mentioned a simple use of the component tree for filtration (removing nodes of the tree whose attribute is below a given threshold). A more advanced use consists in finding the most significant lobes of a given weighted graph F . More precisely, we want to find the N most significant components with respect to either the height, area or volume criterion. By using the tree, this task reduces to the search of the N nodes that have the largest attribute values and are not bound with each other (even transitively) by the inclusion relation. Algorithm 3 (Keep_N_Lobes) performs this task. Its time complexity is in O(sort(n) + m), where m is the number of vertices in the graph, n is the number of component tree nodes and sort(n) is the complexity of the sorting algorithm. At the end of the algorithm, the remaining leaves (more precisely, the pixels which are associated to these leaves) mark the desired significant lobes. For this algorithm, each node must include fields to store its parent and its number of children (but the list of children of a given node is not necessary).

Fig. 9 illustrates this algorithm. Fig. 9.a is an image of cell, in which we want to extract the ten bright lobes. Fig. 9.b shows

Fig. 1 .

 1 Fig. 1. (a) A vertex-weighted graph (V, E, F ) and its cross-sections at levels 1, 2, 3, 4. (b) The component tree of F . (c) The associated component mapping. The component at level 1 is called α 1 , the two components at level 2 are called α 2 and α 3 (according to the usual scanning order), and so on.

  MakeSet(element x) Par(x) := x; Rnk(x) := 0; Function element Find(element x) if (Par(x) = x) then Par(x) := Find(Par(x)); return Par(x); Function element Link(element x, element y) if (Rnk(x) > Rnk(y)) then exchange(x, y); if (Rnk(x) == Rnk(y)) then Rnk(y) := Rnk(y) + 1; Par(x) := y; return y; B. Illustration of union-find: labelling the connected components

4 compq

 4 := Find(q); 5 if (compp = compq) then 6 compp := Link(compq, compp); 7 foreach p ∈ X do M(p) := Find(p);8

Algorithm 2 :

 2 BuildComponentTree Data: (V, E, F ) -vertex-weighted graph with N points. Result: nodes -array [0 . . . N -1] of nodes. Result: Root -Root of the component tree Result: M -map from V to [0 . . . N -1] (component mapping). Local: lowestN ode -map from [0 . . . N -1] to [0 . . . N -1].

1 foreach 5 foreach

 15 p ∈ V do {MakeSet tree (p); MakeSet node (p); nodes[p]:= MakeNode(F (p)); lowestN ode[p] := p;}; 2 foreach p ∈ V in decreasing order of level for F do 3 curTree := Find tree (p);4curNode := Find node (lowestN ode[curTree]); already processed neighbor q of p with F (q) ≥ F (p) do 6 adjTree := Find tree (q); 7 adjNode := Find node (lowestN ode[adjTree]);

12 curTree : =

 12: Link tree (adjTree, curTree); 13 lowestN ode[curTree] := curNode; 14 Root := lowestN ode[Find tree (Find node (0))] ; 15 foreach p ∈ V do M (p) := Find node (p);16

Fig. 2 .

 2 Fig. 2. (a) Original vertex-weighted graph. (b) Points are labelled accordingto the usual lexicographic order, but they will be processed by decreasing level (that is: 12, 0, 2, 1,[START_REF] Jones | Component trees for image filtering and segmentation[END_REF][START_REF] Hartigan | Statistical theory in clustering[END_REF][START_REF] Bertrand | On topological watersheds[END_REF][START_REF] Breen | Attribute openings, thinnings and granulometries[END_REF][START_REF] Chiang | Simple and optimal outputsensitive construction of contour trees using monotone paths[END_REF][START_REF] Couprie | Quasi-linear algorithms for the topological watershed[END_REF][START_REF] Dillencourt | A general approach to connected-component labeling for arbitrary image representations[END_REF][START_REF] Ch | Two linear Union-Find strategies for image processing[END_REF][START_REF] Guillataud | Contribution à l'analyse dendroniques des images[END_REF][START_REF] Cormen | Introduction to Algorithms[END_REF][START_REF] Couprie | Topological grayscale watershed transform[END_REF].

Fig. 3 .

 3 Fig. 3. Beginning of step 6. (a) State of the maps Par tree , Par node and lowestN ode. (b) Partial trees constructed.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Beginning of step 11. (a) State of the maps Par tree , Par node and lowestN ode. (b) Partial trees constructed.

Fig. 6 .

 6 Fig. 6. End of step 14. (a) State of the map Par tree , Par node and lowestN ode. (b) Component tree.

Fig. 7 .

 7 Fig. 7. An example of an artificially generated image of size N*N, where values of pixel (x 1 , x 2 ) with x 1 + x 2 odd are uniformly distributed between 0 and N * N , and where the other half of the pixels are 0. Using a series of such images, one can verify that the component tree algorithm of Salembier et al. is quadratic.

Fig. 8 .

 8 Fig. 8. Illustration of the height, the area and the volume of a component.

Function

  int ComputeVolume(int n) vol := nodes[n]→area; foreach c child of nodes[n] do vol := vol + ComputeVolume(c) + c→area * (c→level -nodes[n]→level); nodes[n]→volume := vol ; return vol;

A preliminary and reduced version of this paper appeared in conference proceedings as[START_REF] Najman | Quasi-linear algorithm for the component tree[END_REF]. This work has been partially supported by the CNRS.

The precise definition of α, a "diagonal inverse" of the Ackermann's function, involves notions which are not in the scope of this paper, it can be found in[START_REF] Tarjan | Efficiency of a good but not linear set union algorithm[END_REF].

Algorithm 3: Keep_N_Lobes

Data: A vertex-weighted graph (V, E, F ), its component tree T with attribute value for each node, and the associated component mapping M Data: The number N of wanted lobes. Result: The filtered map F Sort the nodes of T by increasing order of

that the image 9.a contains numerous maxima. Fig. 9.c is the filtered image obtained by using algorithm 3 with the volume attribute and with parameter value 10, and Fig. 9.d shows the maxima of this filtered image. Note that a similar result could be obtained with this image by performing attribute based operations using several volume threshold values, following e.g. a dichotomic method, until the desired number of maxima is reached. This latter approach is not only less efficient than the proposed algorithm, but it may also fail to find the precise number of maxima required by the user, in the case of components having precisely the same attribute value. In such cases, the proposed algorithm always makes a choice in order to fulfill the user's requirement.

The component tree allows the efficient implementation of complex image and signal filtering, based for example on the use of criteria such as area, volume or depth, or even the use of non-increasing criteria [START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF]. Although some of these filters may be computed using specific and sometimes faster algorithm (in particular area filtering [START_REF] Meijster | A comparison of algorithms for connected set openings and closings[END_REF]), using the component tree is in general the simplest and the most efficient way to compute these filters. Moreover, once the component tree of a function is computed, any of these filters, with any parameter value, can be computed at a very low cost. The component tree is also a key element of an efficient algorithm for the topological watershed [START_REF] Couprie | Quasi-linear algorithms for the topological watershed[END_REF]. New classes of filters, such as second-order connected operators [START_REF] Serra | Connectivity on complete lattices[END_REF] have been recently introduced to generalize connected operators [START_REF] Salembier | Flat zones filtering, connected operators and filter by reconstruction[END_REF].

Those operators can also be efficiently implemented using the component tree [START_REF] Ouzounis | Second-order connected attribute filters using max-trees[END_REF]. In this paper, we have proposed a simple-to-implement quasi-linear algorithm for computing the component tree. We hope that such an algorithm will facilitate the extensive practical use of such operators.


L Najman
email: l.najman@esiee.fr

M Couprie
email: m.couprie@esiee.fr

Quasi-linear algorithm for the component tree

Keywords:

The level sets of a map are the sets of points with level above a given threshold. The connected components of the level sets, thanks to the inclusion relation, can be organized in a tree structure, that is called the component tree. This tree, under several variations, has been used in numerous applications. Various algorithms have been proposed in the literature for computing the component tree. The fastest ones have been proved to run in 0(n ln(n)) complexity. In this paper, we propose a simple to implement quasi-linear algorithm for computing the component tree on symmetric graphs, based on Tarjan's union-find principle.

INTRODUCTION

The level sets of a map are the sets of points with level above a given threshold. The connected components of the level sets, thanks to the inclusion relation, can be organized in a tree structure, that is called the component tree. The component tree captures some essential features of the map. Thus, it has been used (under several variations) in numerous applications among which we can cite: image filtering and segmentation, [START_REF] Hanusse | Sémantique des images par analyse dendronique[END_REF][START_REF] Guillataud | Contribution à l'analyse dendroniques des images[END_REF][START_REF] Couprie | Topological grayscale watershed transform[END_REF][START_REF] Jones | Component trees for image filtering and segmentation[END_REF] video segmentation, 5 image registration, [START_REF] Mattes | Tree representation for image matching and object recognition[END_REF][START_REF] Monasse | Morphological representation of digital images and application to registration[END_REF] image compression [START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF] and data visualization. [START_REF] Chiang | Simple and optimal output-sensitive construction of contour trees using monotone paths[END_REF] We also note that this tree is fundamental for the efficient computation of the topological watershed introduced by M. Couprie and G. Bertrand. [START_REF] Couprie | Topological grayscale watershed transform[END_REF][START_REF] Couprie | Quasi-linear algorithms for topological watershed[END_REF] While having been (re)discovered by several authors for image processing applications, the component tree concept has its root in statistics. [START_REF] Wishart | Mode analysis: A generalization of the nearest neighbor which reduces chaining effects[END_REF][START_REF] Hartigan | Statistical theory in clustering[END_REF] For image processing, the use of this tree in order to represent the "meaningful" information contained in a numerical function can be found in particular, in a paper by Hanusse and Guillataud 1, 2 ; the authors claim that this tree can play a central role in image segmentation, and suggest a way to compute it, based on an immersion simulation. Several authors, such as Vachier, [START_REF] Vachier | Extraction de caractéristiques, segmentation d'images et Morphologie Mathématique[END_REF] Breen and Jones, [START_REF] Breen | Attribute openings, thinnings and granulometries[END_REF] Salembier et al. [START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF] have used this structure in order to implement efficiently some morphological operators (e.g. connected operators, granulometries, extinction functions).

Let us describe informally an "emergence" process that will later help us designing an algorithm for the component tree. We are going to use topographical references. We can see the map as the surface of a relief with the level of a point corresponding to its altitude. Imagine the surface completely covered by water, and that the level of water slowly decreases. Islands (maxima) appear. These islands form the leafs of the component tree. As the level of water decreases, islands will grow, building the branches of the tree. Sometimes, at a given level, several islands will merge into one connected piece. Such pieces are the forks of the tree. We stop when all the water has disappeared.

Various algorithms have been proposed in the literature for computing the component tree, [START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF][START_REF] Breen | Attribute openings, thinnings and granulometries[END_REF][START_REF] Mattes | Efficient algorithms to implement the confinement tree[END_REF] the last reference also contains a discussion about time complexity of the different algorithms. The fastest ones run in 0(n ln(n)) complexity. In this paper, we propose a quasi-linear algorithm for computing the component tree on general symmetric graphs, based on Tarjan's union-find [START_REF] Tarjan | Efficiency of a good but not linear set union algorithm[END_REF] principle. We would like to emphasize that this algorithm is simple to implement. A proof of the complexity of this algorithm is given.

GRAPH STRUCTURE AND THE COMPONENT TREE

Basic notions for graphs

Let E be a finite set (set of vertices, or points). A graph (E, Γ) is defined through a mapping Γ from E to P(E), where P(E) denotes the set of all subsets of E. The mapping Γ associates to each point x of E, the set Γ(x) of points which are adjacent to x. A graph (E, Γ) is symmetric if for all x, y of E, y ∈ Γ(x) implies x ∈ Γ(y). If y ∈ Γ(x), the set {x, y} is called an edge of the graph, and y is called a neighbor of x.

Let (E, Γ) be a symmetric graph. Let X ⊆ E, and let x 0 , x n ∈ X. A path from x 0 to x n in X is an ordered family (x 0 , x 1 , . . . , x n ) of points of X such that x i+1 ∈ Γ(x i ), with i = 0 . . . n -1. Let x, y ∈ X, we say that x is X-connected to y if there exists a path from x to y in X. The relation "is X-connected to" is an equivalence relation. A connected component of X is an equivalence class for the relation "is X-connected to". We say that the a subset

X ⊆ E is connected if X is made of a single connected component. We say that the graph (E, Γ) is connected if E is connected.
In the sequel of the paper, (E, Γ) denotes a connected symmetric graph, such that for any x ∈ E, x does not belong to Γ(x).

Basic notions for weighted graphs

We denote by F(E) the set composed of all maps from E to D, where D can be the set of rational numbers or the set of integers. For a map F ∈ F(E), the triplet (E, Γ, F ) is called a weighted graph. For a point p ∈ E, F (p) is called the weight or level of p. When no ambiguity may occur, we will denote by F the weighted graph (E, Γ, F ).

Let F ∈ F(E), we define F k = {x ∈ E; F (x) ≥ k} with k ∈ Z; F k is called a cross-section of F . A connected component of a section F k is called a (level k) component of F . A level k component of F that does not contain a level (k + 1) component of F is called a (regional) maximum of F .
We define k min = min {F (x), x ∈ E} and k max = max {F (x), x ∈ E}, which represent respectively, the minimum and the maximum level in the map F .

While the notions we are dealing with in this paper are defined for general graphs, we are going to illustrate our work with the case of 2D images that we model by weighted graphs. Let Z denote the set of integers. We choose for E a subset of Z 2 . A point x ∈ E is defined by its two coordinates (x 1 , x 2 ). We choose for Γ the 4-connected adjacency relation defined by, for all x ∈ E, Γ 4 (x) = {y ∈ E \ {x}; |x 1y 1 | + |x 2y 2 | ≤ 1}. We assume that E is connected.

Fig. 1 shows a weighted graph F and four cross-sections of F , between the level k min = 1 and the level k max = 4. The set F 4 is made of two connected components which are regional maxima of F .

1 1 1 1 1 1 1 1 3 3 2 3 4 1 1 3 3 2 3 4 1 1 1 1 1 1 3 1 1 3 3 2 1 1 1 1 4 3 2 2 2 1 1 F 1 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 0 0 0 0 1 0 0 1 1 1 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 F 2 0 0 0 0 0 0 0 0 1 1 0 1 1 0 0 1 1 0 1 1 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 F 3 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 F 4 Figure 1.
A weighted graph F and its cross-sections at levels 1, 2, 3, 4

Component Tree

From the example of Fig. 1, we can see that the connected components of the different cross-sections may be organized, thanks to the inclusion relation, to form a tree structure.

Let F ∈ F(E), let k ∈ D, we denote by C k (F ) the set of all connected components of the cross-section F k , and we define C(F ) as the union of all the sets C k (F ) with k ∈ {k min . . . k max }. The elements of C(F ) are the connected components of the cross-sections of F , we will call them components.

For the sake of simplicity, we suppose that two components belonging to different cross-sections are always different sets of points, that is, are never strictly identicals. The extension of the following definitions to the general case is straightforward.

We define the component tree T (F ) as the rooted tree such that:

• the nodes of the component tree are the elements of C(F ),

• there is an arc from a component c ∈ C k (F ) to a component c ∈ C j (F ) if j = k + 1 and c ⊆ c . We then say that c is the father of c, and we also say that c is a son of c .

The components which have no sons are called leafs, the component which has no father is called the root.

We define the component mapping C as the map from E to C(F ) which associates to each point p ∈ E the node C(p) such that C(p) belongs to C k (F ) with k = F (p), and such that p ∈ C(p). Fig. 2.a shows the component tree of the weighted graph depicted in Fig. 1, and fig. 2.b shows the associated component mapping. The component at level 1 is called c 1 , the two components at level 2 are called c 2 and c 3 (according to the usual scanning order), and so on.

c1 c2 c3 c4 c5 c6 c8 c7 Level 1 Level 2 Level 3 Level 4 c1 c1 c1 c1 c1 c1 c1 c1 c4 c4 c2 c5 c7 c1 c1 c4 c4 c2 c5 c7 c1 c1 c1 c1 c1 c1 c5 c1 c1 c6 c6 c3 c1 c1 c1 c1 c8 c6 c3 c3 c3 c1 c1 c1 c1 c1 c1 c1 c1 (a) (b)

UNION-FIND ALGORITHM FOR THE COMPONENT TREE

Disjoint Sets

The disjoint set problem consists in maintaining a collection Q of disjoint subsets of a set E under the operation of union. Each set X in Q is represented by a unique element of X, called the canonical element. In the following, x and y denote two distinct elements of E. Three operations allow to manage the collection:

• MakeSet(x): add the set {x} to the collection Q, provided that the element x does not already belongs to a set in Q.

• Find(x): return the canonical element of the set in Q which contains x.

• Link(x, y): let X and Y be the two sets in Q whose canonical elements are x and y respectively (x and y must be different). Both sets are removed from Q, their union Z = X ∪ Y is added to Q and a canonical element for Z is selected and returned.

Tarjan [START_REF] Tarjan | Efficiency of a good but not linear set union algorithm[END_REF] proposed a very simple and very efficient algorithm called union-find to achieve any intermixed sequence of such operations with a quasi-linear complexity. More precisely, if m denotes the number of operations and n denotes the number of elements, the worst-case complexity is Θ(mα(m, n)) where α(m, n) is a function which grows very slowly, for all practical purposes α(m, n) is never greater than four.

The implementation of this algorithm is given below in procedure MakeSet and functions Link and Find. Each set of the collection is represented by a rooted tree, where the canonical element of the set is the root of the tree. To each element x is associated a father Fth(x) and a rank Rnk(x). The mappings 'Fth' and 'Rnk' are represented by global arrays in memory. One of the two key heuristics to reduce the complexity is a technique called path compression, that was used by Tarjan to reduce the cost of Find. It consists, while searching for the root r of the tree which contains of x, in considering each element y of the path from x to r (including x), and setting the parent of y to be r. The other key technique, called union by rank, consists in always choosing the root with the greatest rank to be the root representative of the union while performing the Link operation. If the two roots x and y have the same rank, then one of the roots, say y, is chosen arbitrarily to be the root of the union: y becomes the father x; and the rank of y is incremented by one. The rank Rnk(x) is a measure of the depth of the tree rooted in x, and is exactly the depth of this tree if the path compression technique is not used jointly with the union by rank technique. Union by rank avoids creating degenerate trees, and helps keeping the depth of the trees as small as possible. For a more detailed explanation and complexity analysis, see Tarjan's paper. [START_REF] Tarjan | Efficiency of a good but not linear set union algorithm[END_REF] Procedure MakeSet(element x)

Fth(x) := x; Rnk(x) := 0;

Function element Find (element x) if (Fth(x) = x) then Fth(x) := Find(Fth(x)); return Fth(x);

Function element Link (element x, element y) if (Rnk(x) > Rnk(y)) then exchange(x, y); if (Rnk(x) == Rnk(y)) then Rnk(y) := Rnk(y) + 1; Fth(x) := y; return y;

Illustration of union-find: building the connected components

We can illustrate the use of the union-find algorithm on the classical problem of finding the connected components of a subset X of a graph (E, Γ). The algorithm 1 (BuildConnectedComponents) is given below. For a set X, this algorithm returns a map M that gives for each point p, the canonical element M (p) of the connected component of X which contains p.

A first pass (loop 1) is made, for the union-find preprocessing; during this first pass, for each point p of the set X, the set {p} is added to the collection Q of disjoint subsets. Then, a scanning (loop 2) processes all points of X in any order. For each point p, we first find the canonical element of the disjoint set it belongs to (line 3). Then, for each neighbor q of p such that q ∈ X (line 4), we find the canonical element of the disjoint set which contains q (line 5). If p and q are not already in the same disjoint set, that is if the two canonical elements differ (line 6), then p and q are linked (line 7), and one of the two canonical elements is chosen to be the canonical element of their common disjoint set. At the end of the scanning, a simple pass on all the elements of X (loop 8) builds the map M .

Figure 2 .

 2 Figure 2. The component tree (a) of the weighted graph of Fig. 1 and the associated component mapping (b)

foreach q ∈ Γ(p) such that q ∈ X do

Quasi-linear component tree building algorithm

We are now ready to introduce our quasi-linear algorithm for building the component tree T (F ) from a weighted graph F .

To represent a node of T (F ), we use a structure called node containing the level of the node, and the list of nodes which are sons of the current node. For building the component tree, we do not need the reverse link, that is we do not need to know the father of a given node, but let us note that such an information is useful for applications, and can easily be obtained in a linear-time post-processing step.

The algorithm 2 (BuildComponentTree) is given below. We are going to use two collections Q 1 and Q 2 of disjoint sets. After a preprocessing (line 1) for sorting the points by decreasing order of level and for the preparation of the two union-find implementations (line 2), we process the points, starting with the highest ones.

Let us suppose that we have processed a number of levels, and thus that we have already built several parts of the component tree, that is several partial trees. Let us examine a given point p of the current level. The first collection Q 1 is used to tell which partial tree the point p belongs to (line 4). Using subtreeRoot, we can know the node that is the root of this partial tree; this node may result of the union of differents nodes of the same level. The second collection Q 2 of disjoint sets is used to obtain the canonical element of this union (line 5), that we call canonical node of p.

We then look at each neighbor q of level greater or equal to the one of the point p under examination (loop 6). Note that as the graph is symmetric, the Link between two points is done when one of the two points is processed as a neighbor of the other. Thus, we can use the order of scanning of the points, and we only need to examine the "already processed" neighbors of p. Such a neighbor q satisfies F (q) ≥ F (p).

Exactly as we have done for the point p, we search for the canonical node of q (lines 7-8). If the canonical node of p and the canonical node of q differ, that is if the two points are not already in the same node, we have two possible cases:

• either the two canonical nodes have the same level; in this case, that means that these two nodes are in fact part of the same component, and we have to merge the two nodes (line 9 and function MergeNodes).

The merging of nodes of same level is done by the second union-find implementation. The merging relies on the fact that the Link function always chooses one of the two canonical elements of the sets that are to be merged as the canonical element of the merged set.

Once the merging has been done, observe that we will not access anymore to the node that has not been chosen to be the canonical element of the set. We will only have to know to which set it belongs to, and the answer to this question is given by the Find part of the union-find algorithm.

• either the canonical node of q is of higher level, and thus this node becomes a son of the current node (line 10). N n := N n -1; return tmpNode;

In both cases, we have to link the two partial trees, this is done using the first Link implementation (line 11).

We also have to keep track of the node of lowest level for the union of the two partial trees, that we store in the array subtreeRoot (line 12).

At the end of the algorithm, we have to do a post processing to return the desired result. The root of the component tree can easily be found (line 13) using the array subtreeRoot and the two disjoint sets Q 1 and Q 2 . The component mapping can be obtained using the second disjoint set Q 2 (loop 14).

Complexity analysis

Let n denotes the number of points in E, and let m denotes the number of edges of the graph (E, Γ).

The sorting of the points (line 1) can be done in linear time if the weigths are small integers (counting sort [START_REF] Cormen | Introduction to Algorithms[END_REF] ), and in n log(log(n)) if the length of each weight can be stored in a machine memory word (long integers or floating point numbers [START_REF] Andersson | Sorting in linear time?[END_REF] ).

Loop 2 is the preparation for the union-find algorithm. It is obviously linear.

In the function MergeNodes, the merging of the lists of sons can be done in constant time, because we can merge two lists by setting the first member of one list to be the one that follows the last member of the other list. This requires the two lists to be disjoint, which is the case (we are dealing with disjoint sets), and an adequate representation for lists (chained structure with pointers both on first and last element).

The amortized complexity of line 6 is equal to the number m of edges of the graph (E, Γ). The amortized complexity of all calls to the union-find procedures is quasi-linear with respect to m. The building of the component mapping C is obviously linear.

Thus the complexity of the algorithm 2 (BuildComponentTree) is quasi-linear if the sorting step is linear.

EXAMPLE

We are going to illustrate the work of the algorithm on an example. Let us look at the weighted graph of fig. 3.a. The points are numbered (labeled) according to their usual lexicographical order (fig. 3.b). They are examined by the algorithm according to their level, starting with the highest ones. At the begining of the sixth step (end of fifth step), we have already constructed parts of the component tree (fig. 4.b). We have represented the maps Fth1, Fth2, and subtreeRoot (fig. 4.a). For the maps Fth1 and Fth2, the canonical elements appear in white. It should be noted that the subtreeRoot mapping is only used for the canonical elements of Fth1: this explains why the values of subtreeRoot for other elements (in grey) are not updated.

We are going to process points at level 50. The first point at level 50 is the number 3. The point 0 is a neighbor of 3. The node 0 belongs to the partial tree whose root is the node 1 at level 90. Thus the node 3 becomes the father of node 1. Then, the node 3 is linked succesively with nodes numbered 4, 5 and 8. Then point number 9 is examined, and its node is linked with node number 10, the node 9 being chosen as the canonical element. Then, through the point 12, the partial tree whose root is node 13 (level 70) is put under the node 9. We are then at the begining of step 11, and this is illustrated on figure 5.

The point 11 is a neighbor of both points 8 and 10. The node 8 is on the level 50 component whose root is 3. Thus, nodes 11 and nodes 3 are linked, and the node 3 is chosen as the canonical element. The node 10 is on the level 50 component whose root is 9. Thus, nodes 9 and 3 are merged, that is the corresponding partial trees become a single partial tree. The node 9 is chosen as the canonical element of the level 50 component, and the sons of node 3 are transfered to node 9. We are in the situation depicted in figure 6.

We then process the point 6 at level 40, whose node becomes the father of node 9 at level 50. But the union-find algorithm chooses 9 as the canonical element for the partial tree whose root is 6. We thus have to keep the root of the partial tree in subtreeRoot by setting subtreeRoot[9] := 6. Then we process the point 7 at level 20, and the node 7 becomes the father of node 6. Once again, the union-find algorithm choose 9 as the canoncial element, and thus we have to keep the root of the partial tree by setting subtreeRoot[9] := 7. There is no point higher than 20, and thus, the component tree is built. The final situation is depicted in figure 7. The first collection Q 1 of disjoint sets is not useful anymore: indeed, each node of the graph has been examined, and they are all linked, the canonical element being the node 9. The root of the component tree is the node 7. The second collection Q 2 gives a canonical element for each of the components of F : observe in particular the level 50, whose canonical element is the node 9. The collection Q 2 can be used to compute the component mapping C.


Laurent Najman
email: l.najman@esiee.fr

Michel Couprie
email: m.couprie@esiee.fr

Watershed algorithms and contrast preservation

Keywords: Mathematical Morphology, Watersheds, Contours Saliency, Topology

This paper is devoted to the study of watershed algorithms behavior. Through the introduction of a concept of pass value, we show that most classical watershed algorithms do not allow the retrieval of some important topological features of the image (in particular, saddle points are not correctly computed). An important consequence of this result is that it is not possible to compute sound measures such as depth, area or volume of basins using most classical watershed algorithms. Only one watershed principle, called topological watershed, produces correct watershed contours.

Introduction

This paper is a first of a series dedicated to the notion of watershed contour saliency. Using this concept, introduced in [START_REF] Najman | Morphologie Mathématique: de la Segmentation d'Images à l'Analyse Multivoque[END_REF][START_REF] Najman | Geodesic saliency of watershed contours and hierarchical segmentation[END_REF], we can sum up in one image all the contour information that we can obtain by filtering the image by attribute opening [START_REF] Vachier | Extraction de caractéristiques, segmentation d'images et Morphologie Mathématique[END_REF][START_REF] Breen | Attribute openings, thinnings and granulometries[END_REF][START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF][START_REF] Meijster | A comparison of algorithms for connected set openings and closings[END_REF] for all values of the parameter and applying a watershed on each of the filtered images. Several algorithms [START_REF] Najman | Geodesic saliency of watershed contours and hierarchical segmentation[END_REF][START_REF] Meyer | The dynamics of minima and contours[END_REF][START_REF] Lemaréchal | Comments on 'geodesic saliency of watershed contours and hierarchical segmentation[END_REF][START_REF] Schmitt | Response to the comment "geodesic saliency of watershed contours and hierarchical segmentation[END_REF] for computing saliency of watershed contours have been proposed. We expect to obtain the same result either by thresholding the saliency image at a given level k, or by filtering the original image using an attribute opening using k as parameter value and applying a watershed algorithm on the filtered image. None of the existing saliency algorithms computes this expected result. The goal of the series is to show why it is the case, and to propose a novel efficient algorithm that computes the expected result.

This paper is devoted to the study of watershed algorithms behavior with respect to what is needed to compute saliency of contours. For computing saliency, one needs a map of watershed basins neighborhood, with the altitude of their associated saddle points, and a valuation on each basin. A review of watershed algorithms and their associated results can be found in [START_REF] Roerdink | The watershed transform: Definitions, algorithms and parallelization strategies[END_REF]. This review does not study algorithms from the point of view of the preservation of important topological features of the original image; for instance, it does not consider the question: does the algorithm compute correct saddle points? We are going to tackle the difficult notion of saddle point through the introduction of a concept of "pass value". We demonstrate that watershed algorithms that are the most used in practice do not behave correctly with respect to the preservation of pass values, and thus cannot be used in a saliency algorithm. We show that the approach called topological watershed [START_REF] Couprie | Topological grayscale watershed transform[END_REF] (which is not mentioned in [START_REF] Roerdink | The watershed transform: Definitions, algorithms and parallelization strategies[END_REF]) provides the only existing algorithm that produces a correct entry point for a saliency algorithm.

2 Brief description of watershed algorithms

Intuitive notions for watershed

The intuitive idea underlying the watershed notion comes from the field of topography: a drop of water falling on a relief follows a descending path and eventually reaches a minimum. Watershed lines are the divide lines of the domains of attraction of drops of water. This intuitive approach is not well suited to practical implementations, and can yield biased results in some cases [START_REF] Vincent | Watersheds in digital spaces: An efficient algorithm based on immersion simulations[END_REF]. An alternative approach is to imagine the surface being immersed in a lake, with holes pierced in local minima. Water will fill up basins starting at these local minima, and, at points where waters coming from different basins would meet, dams are built. As a result, the surface is partitioned into regions or basins separated by dams, called watershed lines.

What is a watershed algorithm

This paper is not the place to describe in details the (large) family of watershed algorithms. Nevertheless, it is worthwhile to give a brief description of the main algorithms. Let E be a set of vertices (or points). Let P(E) denote the set of all subsets of E. Let G = (E, Γ ) be a (symmetric) graph, where Γ is a mapping from E into P(E), which associates to each point x of E, the set Γ (x) of points adjacent to x.

Let X ⊆ E, and let x 0 , x n ∈ X. A path from x 0 to x n in X is an ordered family (x 0 , x 1 , . . . , x n ) of points of X such that x i+1 ∈ Γ (x i ), with i = 0 . . . n-1. Let x, y ∈ X, we say that x is connected to y if there exists a path from x to y in X. The relation "is connected to" is an equivalence relation. A connected component of X is an equivalence class for the relation "is connected to".

Let E = Z 2 . We denote by F(E) the set composed of all functions from E to Z. Let F ∈ F(E). We denote by X the complement of X. We write

F k = {x ∈ E; F (x) ≥ k} with k ∈ Z; F k is called an upper (cross-) section of F , and F k is called a lower (cross-) section of F . A connected component of a section F k is called a (level k) lower-component of F . A level k lower-component of F that does not contain a level (k -1) lower-component of F is called a (regional) minimum of F .
Let us recall that a partition of a set S is a collection of non-empty disjoint subsets of S whose union is S.

A watershed algorithm builds a partition of the space:

it associates an influence zone B(M ) called catchment basin, to each minimum M of the image. The set B(M ) is connected and contains M ; it may produce a set of watershed lines which separates those catchment basins one from each other.

2.3 Vincent-Soille watershed algorithm [START_REF] Vincent | Watersheds in digital spaces: An efficient algorithm based on immersion simulations[END_REF] For any set A and any set B ⊂ A made of several connected components B i , the geodesic influence zone iz A (B i ) of B i in A is the locus of the points of A whose geodesic distance to B i is strictly smaller than their geodesic distance to any other component of B. We define the following recursion:

X hmin+1 = F hmin+1 = M IN hmin (1) X h+1 = M IN h ∪ IZ F h+1 (X h ) (2)
where h min is the lowest grey-value of F, where IZ F h+1 (X h ) is the union of the geodesic influence zones of the connected components of X h in F h+1 , and where M IN h is the union of minima of F with grey-level equal to h. The watershed lines are the complement of X hmax+1 . As noted in [START_REF] Roerdink | The watershed transform: Definitions, algorithms and parallelization strategies[END_REF], Vincent-Soille's algorithm does not implement exactly this recursion. Thanks to a fifo queue, it floods the catchment basins of the image, and to build the watershed lines, it associates a special value WSHED to the pixels where two different cacthment basins would merge. A point labelled WSHED by the algorithm is not considered again in the following iteration, as it should be the case. Furthermore, pixels labelled WSHED are propagated. This allows the detection of special thick watershed zones, like those called buttonholes (see fig. 2.a).

2.4

Meyer's watershed algorithm [START_REF] Meyer | Un algorithme optimal de ligne de partage des eaux[END_REF] Starting from a greyscale image F and a set M of markers with different labels (in our case, these will be the minima of F ), it expands as much as possible the set M , while preserving the number of connected components of M :

1. insert every neighbor x of every marked area in a hierarchical queue, with a priority level corresponding to the grey level F (x). Note that a point cannot be inserted twice in the queue; 2. extract a point x from the hierarchical queue, at the highest priority level, that is, the lowest grey level. If the neighborhood Γ (x) of x contains only points with the same label, then x is marked with this label, and its neighbors that are not yet marked are put into the hierarchical queue;

Step 2 must be repeated until the hierarchical queue is empty. The watershed lines set is the complement of the set of labeled points. Let us note that this algorithm does neither label nor propagate watershed pixels, which "stop" the flooding. Thus, the watershed lines produced by Meyer's algorithm are always thinner than lines produced by other watershed algorithms.

Cost-based watershed

The principle is to define a distance or a cost for travelling between pixels, and to define the influence zone of a minimum as the set of points which are strictly closer to this minimum than to any other minimum. Various costs or distances can be considered, the most popular ones being the topographical distance [START_REF] Meyer | Topographic distance and watershed lines[END_REF][START_REF] Najman | Watershed of a continuous function[END_REF], but other approaches exist [START_REF] Lotufo | Ift-watershed from gray-scale marker[END_REF] among which we can mention the max-arc path cost. The (so-called) topographical distance of an image F is a digital analogue to dt F (x, y) = inf π∈Π(x,y) π ||∇F (π(s))||ds. Let us note that if we are on a line of steepest slope between x and y, then dt F (x, y) = |F (x) -F (y)|. The catchment basin of a minimum m i is defined as the set of pixels x for which F (m i ) + dt F (m i , x) < F (m j ) + dt F (m j , x) for all minima m j = m i . The watershed lines set is the complementary of those catchment basins.

Another simple possible choice is the max-arc path cost [START_REF] Lotufo | Ift-watershed from gray-scale marker[END_REF] which assigns to a path the maximum of F for each pixel on the path. In this case, dm F (x, y) = inf π∈Π(x,y) max i F (π(i)), and the catchment basin of a minimum m i is defined as the set of pixels x for which dm F (m i , x) < dm F (m j , x) for all minima m j = m i . The watershed lines set is the complementary of those catchment basins.

2.6

Topological watershed [START_REF] Couprie | Topological grayscale watershed transform[END_REF] The idea is to define a transform that acts directly on the greyscale image, by lowering some points in such a manner that the connectivity of each lower cross-section F k is preserved. The regional minima of the result, which have been spread by this transform, can be interpreted as the catchment basins. The formal definition relies on a particular notion of simple point:

Definition 1. Let G = (E, Γ ) be a graph, and let X ⊂ E.

The point x ∈ X is simple (for X) if the number of connected components of X ∪ {x} equals the number of connected components of X. In other words, x is simple (for X) if x is adjacent to exactly one connected component of X.

We can now define a notion of destructible point, and the topological watershed:

Definition 2. Let F ∈ F(E), x ∈ E, and k = F (x). The point x is destructible (for F ) if x is simple for F k .
We say that W ∈ F(E) is a topological watershed of F if W may be derived from F by iteratively lowering destructible points by one until stability (that is, until all points of E be non-destructible for W ).

The catchment basins of the topological watershed W are the minima of W , and the watershed lines are the non-minima of W .

As a consequence of this definition, a topological watershed W of a function F is a function which has the same number of regional minima as F . Furthermore, the connectivity of any lower cross-section is preserved during this transformation.

Let us note that, in this case, and contrary to other watershed principles, the watershed lines are part of the definition: it does not exist a variation of this notion that does not build those lines. An efficient algorithm to compute the topological watershed has been proposed in [START_REF] Couprie | Topological grayscale watershed transform[END_REF].

Let us emphasize the essential difference between this notion of topological watershed and the notion of homotopic greyscale skeleton, pioneered by Goetcherian [START_REF] Goetcherian | From binary to grey tone image processing using fuzzy logic concepts[END_REF] and extensively studied in [START_REF] Bertrand | Image segmentation through operators based upon topology[END_REF][START_REF] Couprie | Topological operators for grayscale image processing[END_REF]. With the topological watershed, only the connected components of the lower cross-sections of the function are preserved, while the homotopic greyscale skeleton preserves both these components and the components of the upper cross-sections. As a consequence, an homotopic greyscale skeleton may be computed by using a purely local criterion for testing whether a point may be lowered or not, while computing a topological watershed requires the use of a global data structure [START_REF] Couprie | Topological grayscale watershed transform[END_REF].

Watershed algorithms comparison

Intuitively, for application to image analysis, the watershed lines represent the location of pixels which best separate the dark objects (regional minima), in terms of grey level difference (contrast). In order to evaluate the effectiveness of this separation, we have to consider the values of pixels along watershed lines. This motivates the following definition. Definition 3. The watershed contours of F is a grayscale image W such that W (x) = 0 for any x in a catchment basin, and W (x) = F (x) elsewhere.

Let us note that such a definition is not necessary for the topological watershed, which produces a function, and not a binary result.

Saddle point, pass value and the dynamics

To formalize the notion of contrast between two minima, we need to characterize first-contact points between basins. In the continuous framework, such points are called saddle points, but this notion is difficult to transfer to the digital grid. Furthemore, such a notion is not fundamental for contrast criteria. More precisely, for each couple of neighboring basins, we only need the altitude of the lowest contact point between them. This is the motivation for defining the pass value, a natural concept already used by several authors. Definition 4. Let M(F ) be the set of all minima of F . We define the pass value F (m 1 , m 2 ) between two minima m 1 and m 2 in M(F ) as

F (m 1 , m 2 ) = min π∈Π(m1,m2) max i F (π(i)) (3)
where Π(m 1 , m 2 ) is the set of all paths linking m 1 to m 2 .

For applications to image analysis like filtering, and especially for saliency, we want to compute all pass values; we would like watershed contours to have the same pass values as the original image;

and a measure of contrast or importance of each basin (minima) of the original image; such a measure should correspond to measure taken on lower cross-sections of the original image.

Various contrast measures can be computed, among which we can mention depth (dynamics [START_REF] Grimaud | A new measure of contrast: Dynamics[END_REF]), area and volume [START_REF] Vachier | Extraction de caractéristiques, segmentation d'images et Morphologie Mathématique[END_REF]. We are going to examine more particularly the case of the dynamics. We first recall the basic definitions introduced by Grimaud [START_REF] Grimaud | A new measure of contrast: Dynamics[END_REF] (in fact these definitions were proposed for 2D images, we extend them for arbitrary graphs). Let F ∈ F(E) and let X be a minimum for F . The attraction domain of X is the set composed of all point x such that there exist a descending path from x to X. The attraction domain of a minimum X is denoted by K(X). Let π be a path. The dynamics of π (for F ) is the value Dyn(π) = M ax{|F (x)-F (y)|; for all x, y in π}. Let x, y be two points. The dynamics between x and y (for F ) is the value Dyn(x, y) = M in{Dyn(π); for all π ∈ Π(x, y)}. Let X and Y be two subsets of E. The dynamics between X and Y (for F ) is the value Dyn(X, Y ) = M in{Dyn(x, y); for all x ∈ X, y ∈ Y }. Definition 5. Let X ∈ M(F ). The dynamics of X (for F ) is the number Dyn(X) such that:

-If F (X) = M in{F (Y ); Y ∈ M(F )}, then Dyn(X) = ∞; -Otherwise, Dyn(X) = M in{Dyn[X, K(Y )]; ∀Y ∈ M(F ), F (Y ) < F (X)}.

The case of the topological watershed

We can prove [START_REF] Najman | Topological watershed and contrast preservation[END_REF] that the topological watershed preserves the pass values.

Property 1 Let W be a topological watershed of F . For all (m 1 , m 2 ) ∈ M(F ) 2 , and for the corresponding minima (m 1 , m 2 ) ∈ M(W ) 2 , we have

F (m 1 , m 2 ) = W (m 1 , m 2 )
In the sequel of the paper, we are going to show that this property is neither true for the Vincent-Soille's algorithm, nor for Meyer's algorithm. An important consequence of this property is that measures (such as depth (dynamics [START_REF] Grimaud | A new measure of contrast: Dynamics[END_REF]), area or volume [START_REF] Vachier | Extraction de caractéristiques, segmentation d'images et Morphologie Mathématique[END_REF]) computed on the basins obtained by either Vincent-Soille's algorithm or Meyer's algorithm do not correspond to measures of connected components of lower-cross sections of the image. On the contrary, the topological watershed does allow such computations. In particular, in the case of the dynamics, a consequence of property 1 is the following result.

Property 2 Let F ∈ F(E) and let W be a topological watershed of F . Then the dynamics of a minimum for F is equal to the dynamics of the corresponding minimum for W .

We can also prove that, for suitable cost functions, cost based watersheds of an image F preserve the pass values of F . But, as we will see, cost based watershed produces very thick contours that prevent it to be used for a saliency algorithm.

Comparison and counter examples for other watershed algorithms

We are going to examine the behavior of watershed algorithms on several examples. In the sequel, the watershed examples are computed in 4-connectivity.

In particular, regional minima are 4-connected subsets of Z Z 2 . On all the pictures, the basins are labeled with letters, and the watershed pixels are given with their corresponding value in the original image. Similar configurations can be found for other connectivities. Let us emphasize that configurations similar to the examples presented in this paper were found in real images. Neither Vincent-Soille's nor Meyer's algorithm do preserve the pass values. A counter-example that illustrates this behavior is given in Figure 1. Figure 1.a presents a high contour at altitude 255. This contour is run over by the flooding principle of both Meyer and Vincent-Soille. This is especially visible on Meyer's algorithm, as in figure 1.c, the pass value between E and any other minima is 6 instead of 255. Vincent-Soille's algorithm, while having the same kind of problem, tries to detect special pixel configurations called buttonholes, and thus produces thick lines. But in this case, the Vincent-Soille's watershed is not thick enough, and the pass value between E and any other basin is 7 for the watershed contours, while it is 255 for the original image. The only correct result is produced both by the topological watershed and the cost based watershed, and is presented in figure 1.d.

Vincent-Soille's watershed algorithm aims at detecting watershed areas such as buttonholes. These areas are such that one cannot decide towards which minimum a drop falling on them will slide. Figure 2.a exhibits a particular case of buttonhole. Clearly, the pixels at altitude 20 are essential since they carry the pass value between the minimum A (level 2) and the minima B and C (level 1 and 0). We can observe on figures 2.b and 2.c that both Meyer's algorithm and Vincent-Soille's remove the contour at altitude 20; in fact, Meyer's algorithm does not "see" at all this buttonhole. In both cases, pass value between A and B or C is at an altitude of 10 instead of 20 for the watershed contours.

(a) A A A A A A A A A A A A A A A A A A A A A 30 A A

(b) A A A A A A A A A A A A A

In order to preserve pass values on the buttonhole, we have two possibilitites:

either keeping in the watershed lines all the pixels of the buttonhole: that is what is done by the cost-based watershed (Fig. 2.d), which produces contours that cover the whole buttonhole; or making a careful (but arbitrary) choice between all the contours possible in the buttonhole, the choice being such that it preserves the pass values. This is what is done by the topological watershed (Fig. 2.e).

On real images, both cost based and Vincent-Soille's watershed are very sensitive to buttonholes and the resulting watershed lines can cover a large part of the image [example not shown due to space constraint]. Meyer's algorithm and the topological watershed compute thinner lines.

Furthermore, cost-based watersheds produce very thick lines even in the absence of buttonhole, as noted in [START_REF] Roerdink | The watershed transform: Definitions, algorithms and parallelization strategies[END_REF], and tend to isolate basins. Figure 3 illustrates this problem. Indeed, those algorithms have been designed to compute basins, and not lines. Thus, they cannot be used as an entry for a saliency algorithm. Let us note that all watershed algorithms can produce thick watershed lines in some configurations (for instance, think of 4 lines crossing at one point).

Consequences and conclusion

In this paper, we have shown that -Meyer's and Vincent-Soille's algorithms do not preserve important topological features of the image; in particular, pass values are not correct. Only cost-based and topological watersheds are correct from this point of view; furthermore, cost-based watershed and Vincent-Soille's algorithm can produce very thick watershed lines.

Thus, only one watershed notion, the topological watershed, is suited to our task: the associated algorithm is the only one that produces a correct basin neighborhood map and correct pass values.

For computing saliency, we need a measure of contrast of the watershed basins, such as depth (dynamics [START_REF] Grimaud | A new measure of contrast: Dynamics[END_REF]), area or volume [START_REF] Vachier | Extraction de caractéristiques, segmentation d'images et Morphologie Mathématique[END_REF]. An important consequence of the results of this paper is that measures computed on the basins obtained by either Vincent-Soille's algorithm or Meyer's algorithm do not correspond to measures of connected components of lower-cross sections of the image. On the contrary, we have seen that the topological watershed allows such computations. Thus, it is not possible to use the propagation mechanism of the line-building versions of Meyer's or Vincent-Soille's algorithms to compute "on-the-fly" such a measure. Such a mechanism was implemented in the Najman-Schmitt's saliency algorithm [START_REF] Najman | Morphologie Mathématique: de la Segmentation d'Images à l'Analyse Multivoque[END_REF][START_REF] Najman | Geodesic saliency of watershed contours and hierarchical segmentation[END_REF], and has also been proposed in [START_REF] Andrade | A topological image segmentation method by attributes and applications[END_REF], leading to non-correct results.

We could think that we can correct past saliency algorithms by replacing their watershed operator by the topological watershed. Unfortunately, this is not enough. In future papers of the series, we are going to review past saliency algorithms, to show what hypothesis they rely on are wrong, and propose a novel efficient saliency algorithm.

Fig. 1 .

 1 Fig. 1. Counter-example to pass-values preservation. A greyscale image (a) and some results of watershed algorithms: (b) Vincent-Soille (c) Meyer and (d) Cost-based and Topological watershed. One can see that the pass value between E and any other basin is 6 in (c) and is 7 in (b). Both the cost-based and the topological watershed (d) do preserve the correct pass value of 255

Fig. 2 .

 2 Fig. 2. Another counter-example to pass values preservation. (a): Original "buttonhole" image, (b): Meyer's watershed contours, (c): Vincent-Soille's watershed contours, (d): cost-based watershed contours, (e): Topological watershed contours. One can note that the contour at altitude 20 is neither kept by Vincent-Soille's algorithm, nor by Meyer's algorithm. One can also note that both the cost-based and the topological watersheds preserve the pass values of the buttonhole (a), but the topological watershed (e) is thinner than the cost-based watershed (d)

Fig. 3 .

 3 Fig. 3. A greyscale image (a) and some results of watershed algorithms:(b) cost-based watershed and (c) result according to Vincent-Soille, Meyer or Topological watershed. Basin B is isolated in (b)