Thierry Géraud

Edwin Carlinet

Sébastien Crozet

Laurent Najman
email: l.najman@esiee.fr

A quasi-linear algorithm to compute the tree of shapes of nD images

To compute the morphological self-dual representation of images, namely the tree of shapes, the state-of-the-art algorithms do not have a satisfactory time complexity. Furthermore the proposed algorithms are only effective for 2D images and they are far from being simple to implement. That is really penalizing since a self-dual representation of images is a structure that gives rise to many powerful operators and applications, and that could be very useful for 3D images. In this paper we propose a simple-to-write algorithm to compute the tree of shapes; it works for nD images and has a quasi-linear complexity when data quantization is low, typically 12 bits or less. To get that result, this paper introduces a novel representation of images that has some amazing properties of continuity, while remaining discrete.

3 We can indifferently use the term "cut" or "threshold".

Introduction

The tree of shapes [START_REF] Monasse | Fast computation of a contrast invariant image representation[END_REF] is an important morphological structure that represents images in a self-dual way. Shortly put it can be seen as the result of merging the pair of dual component trees, min-tree and max-tree, into a single tree. Using the tree of shapes has many advantages. Since it is self-dual, it makes no assumption about the contrast of objects (either light object over dark background or the contrary). We only have one structure that represents the image contents so we do not have to juggle with the couple of dual trees. It intrinsically eliminates the redundancy of information contained in those trees. Last, it encodes the spatial inclusion of connected components in gray-level images so it is complementary to some other representations that focus on component (or region) adjacency. As a consequence the tree of shapes is not only an easy access to self-dual operators such as grain filters but it has many applications, as listed in [START_REF] Meinhardt-Llopis | Morphological and Statistical Techniques for the Analysis of 3D Images[END_REF] (pp. [START_REF] Meyer | Un algorithme optimal de ligne de partage des eaux[END_REF][START_REF] Monasse | Fast computation of a contrast invariant image representation[END_REF][START_REF] Najman | Discrete set-valued continuity and interpolation[END_REF], and some very recent works illustrate several powerful perspectives offered by that tree (see [START_REF] Xu | Context-based energy estimator: Application to object segmentation on the tree of shapes[END_REF][START_REF] Xu | Morphological filtering in shape spaces: Applications using tree-based image representations[END_REF][START_REF] Xu | Two applications of shape-based morphology: Blood vessel segmentation and generalisation of constrained connectivity[END_REF], and their bibliography).

In the following we consider a nD digital image u as a function defined on a regular cubical grid (precisely, u : Z n → Z), and to properly deal with some subsets of Z n and with their complementary, we consider the dual connectivities c 2n and c 3 n -1 . For any λ ∈ Z, the lower (strict) cuts 3 and upper (large) cuts of u are defined as [ u < λ ] = { x ∈ X | u(x) < λ } and [ u ≥ λ ] = { x ∈ X | u(x) ≥ λ }. From them we deduce two sets, T < (u) and T ≥ (u), composed of the connected components of respectively lower and upper cuts of u: T < (u) = { Γ ∈ CC c2n ([ u < λ ]) } λ and T ≥ (u) = { Γ ∈ CC c 3 n -1 ([ u ≥ λ ]) } λ , where CC denotes the operator that gives the set of connected components of a set. The elements of T < (u) and T ≥ (u) respectively give rise to two dual trees: the min-tree and the max-tree of u. We then define two other sets, S < (u) (set of lower shapes) and S ≥ (u) (set of upper shape), as the sets of components of resp. T < (u) and T ≥ (u) after having filled the cavities4 of those components. With the cavity-filling (or saturation) operator denoted by Sat, we have: The set of all shapes S(u) = S < (u) ∪ S ≥ (u) forms a tree, the so-called tree of shapes of u [START_REF] Monasse | Fast computation of a contrast invariant image representation[END_REF]. Indeed, for any pair of shapes Γ and Γ in S, we have Γ ⊂ Γ or Γ ⊂ Γ or Γ ∩ Γ = ∅. Actually, the shapes are the cavities of the elements of T < and T ≥ . For instance, if we consider a lower component Γ ∈ [ u < λ ] and a cavity H of Γ , this cavity is an upper shape, i.e., H ∈ S ≥ . Furthermore, in a discrete setting, H is obtained after having filled the cavities of a component of [ u ≥ λ ]. Figure 1 depicts on a sample image the three components trees (T < , T ≥ , and S). Just note that the Equations so far rely on the pair of dual connectivities, c 2n and c 3 n -1 , so discrete topological problems are avoided, and, in addition, we are forced to consider two kind of cuts: strict ones for c 2n and large ones for c 3 n -1 .

S < (u) = { Sat c 3 n -1 (Γ ); Γ ∈ T < (u) } and S ≥ (u) = { Sat c2n (Γ ); Γ ∈ T ≥ (u) }. D E B A C F O min-tree max-tree tree of shapes 3 2 0 1 A C, D, E O B F B A, F C D E O 3 2 0 1 A O F B C D E > 2 > 1 > 2 > 2 > 0 > 1 > 0 > 3
The state-of-the-art of tree of shapes computation (detailed in Section 5) suffers from two major flaws: existing algorithms have a time complexity of O(n 2 ) and they cannot easily be extended to nD images. Briefly put, this is due to the fact that either they follow shape contours or they have to know if a component has a cavity 4 . This paper presents an algorithm that can compute the tree of shapes with quasi-linear time complexity when image data are low quantized; furthermore this algorithm straightforwardly applies to nD images.

This paper is organized as follows. First we explain that a well-known algorithmic scheme can be reused to compute the tree of shapes (Section 2). Then this paper introduces a new discrete representation of images (Section 3) that has some properties borrowed from the continuous world. At that point we are ready to glue together the algorithmic scheme and the novel image representation to present a quasi-linear algorithm that compute the tree of shapes (Section 4). Related works about that tree computation is presented so that the reader can compare our approach to existing ones (Section 5). Last, we give a short conclusion (Section 6)5 .

Algorithmic scheme and the need for continuity

This section shows that the max-tree algorithm presented in [START_REF] Berger | Effective component tree computation with application to pattern recognition in astronomical imaging[END_REF] is actually an algorithmic "canvas" [START_REF] Géraud | Mathematical Morphology-From Theory to Applications[END_REF], that is, a kind of meta-algorithm that can be "filled in" so that it can serve different aims. In the present paper it gives an algorithm to compute the tree of shapes.

About union-find and component trees

An extremely simple union-find structure (attributed by Aho to McIlroy and Morris) was shown by Tarjan [START_REF] Tarjan | Efficiency of a good but not linear set union algorithm[END_REF] to be very efficient. This structure, also called disjoint-set data structure or merge-find set, has many advantages that are detailed in [START_REF] Carlinet | A (fair?) comparison of many max-tree computation algorithms[END_REF]; amongst them, memory compactedness, simplicity of use, and versatility. This structure and its related algorithms are of prime importance to handle connected operators [START_REF] Meijster | A comparison of algorithms for connected set openings and closings[END_REF][START_REF] Géraud | Ruminations on tarjan's union-find algorithm and connected operators[END_REF].

Let us denote by R the ancestor relationship in trees: we have a R p iff a is an ancestor of p. R can be encoded as an array of elements (nodes) so that a R p ⇔ index R (a) < index R (p); browsing that array thus corresponds to a downwards browsing of the tree, i.e., from root to leaves. To construct the maxtree of a given image, we rely on a rooted tree defined by a parenthood function, named parent, and encoded as an nD image (so parent(p) is an nD point). When a node of the max-tree contains several points, we choose its first point Algorithm 1: "Union-Find"-based computation of a morphological tree.

union find(R) : T begin for all p do zpar (p) ← undef for i ← N -1 to 0 do p ← R[i] parent(p) ← p zpar (p) ← p for all n ∈ N (p) such as zpar (n) = undef do r ← find root(zpar , n) if r = p then parent(r) ← p zpar (r) ← p return parent find root(zpar , x) : P begin if zpar (x) = x then return x else zpar (x) ← find root(zpar , zpar (x)) return zpar (x) compute tree (u) : Pair(Array[P], T) begin R ← sort(u) parent ← union find(R) canonicalize tree(u, R, parent) return (R, parent)
(with respect to R) as the representative for this node; that point is called a component "canonical point" or a "level root". Let Γ denote a component corresponding to a node of the max-tree, p Γ its canonical element, and p r the root canonical element. The parent function that we want to construct should verify the following four properties:

1. parent(p r ) = p r ; 2. ∀ p = p r , parent(p) R p ; 3. p is a canonical element iff p = p r ∨ u(parent(p)) = u(p) ; 4. ∀ p, p ∈ Γ ⇔ u(p) = u(p Γ ) ∧ ∃ i, parent i (p) = p Γ (therefore ∀ p ∈ Γ, p = p Γ ∨ p Γ R p).
The routine union find, given in Algorithm 1, is the classical "union-find" algorithm [START_REF] Tarjan | Efficiency of a good but not linear set union algorithm[END_REF] but modified so that it computes the expected morphological tree [START_REF] Berger | Effective component tree computation with application to pattern recognition in astronomical imaging[END_REF] while browsing pixels following R -1 , i.e., from leaves to root (let us recall that we do not feature here the union-by-rank version). Its result is a parent function that fulfills those first four properties. Obtaining the following extra property, "5. ∀p, parent(p) is a canonical element," is extremely interesting since it ensures that the parent function, when restricted to canonical elements only, gives a "compact" morphological tree such as the ones depicted in Figure 1. Precisely it allows to browse components while discarding their contents: a traversal is thus limited to one element (one pixel) per component, instead of passing through every image elements (pixels). Transforming the parent function so that property 5 is verified can be performed by a simple post-processing of the union-find computation. The resulting tree has now the simplest form that we can expect; furthermore we have an isomorphism between images and their canonical representations.

Computing the max-tree and the tree of shapes

The algorithm presented in [START_REF] Berger | Effective component tree computation with application to pattern recognition in astronomical imaging[END_REF] to compute the max-tree is surprisingly also able to compute the tree of shapes. The skeleton, or canvas, of this algorithm is the routine compute tree given in the right part of Algorithm 1; it is composed of three steps: sort the image elements (pixels); then run the modified unionfind algorithm to compute a tree encoded by a parent function; last modify the parent function to give that tree its canonical form.

In the case of the max-tree, the sorting step provides R encoded as an array of points sorted by increasing gray-levels in u, i.e., such that the array indices satisfy

i < i ⇒ u(R[i]) ≤ u(R[i ]
). When image data are low quantized, typically 12 bit data or less, then sorting points can be performed by a distribution sort algorithm. Last, the canonicalization post-processing is a trivial 5-line routine that the reader can find in [START_REF] Berger | Effective component tree computation with application to pattern recognition in astronomical imaging[END_REF]. In the case of the tree of shapes, it is also a tree that represents an inclusion relationship between connected components of the input image. As a consequence a first important idea to catch is that the tree of shapes can be computed with the exact same routine, union find, as the one used by max-tree.

What if...

The major and crucial difference between the max-tree and the tree of shapes computations is obviously the sorting step. For the union find routine to be able to compute the tree of shapes using R -1 , the sort routine has to sort the image elements so that R corresponds to a downward browsing of the tree of shapes. Schematically we expect that R contains the image pixels going from the "external" shapes to the "internal" ones (included in the former ones). Fig. 2. Tree computation of the max-tree (left) and of the tree of shapes (right). For both cases, the result R of the sorting step is given over the green arrow and the tree computation, browsing R -1 , is progressively depicted.

The similarity between the computations of both trees is illustrated in Figure 2. We can see that the modified union-find algorithm correctly computes both trees once R is properly defined. Therefore we "just" need to know how to compute R in the case of the tree of shapes to turn the canvas given in the previous Section 2.2 as the expected algorithm. Let us consider the image depicted on the left of Figure 3 with its tree of shapes. We can see that we need to reach the regions A and A' before the regions B and C in order to properly sort pixels, i.e., to compute R. It is only possible if we can pass "between" pixels. The representation depicted on the right of Figure 3 is well-suited for that since it contains some elements that materialize inter-pixel spaces. Furthermore, given a two adjacent pixels with respective values 0 and 3, the element in-between them has to bear all the "intermediate" values: not only 1 but also 2. Indeed, if we change the value of regions A and A' from 1 to 2, the tree structure is unchanged but inter-pixel elements between regions B and C have now to make A and A' connect with value 2. Eventually we need an image representation that is "continuous" in some way with respect to both the domain space and the value space.

Image representation

To be able to sort the image pixels so that R corresponds to a top-down browsing of tree of shapes elements, this paper introduces a novel representation of images 6 . It relies on a couple of theoretical tools briefly described hereafter7 .

Cellular complex and Khalimsky grid

From the sets H 1 0 = {{a}; a ∈ Z} and H 1 1 = {{a, a + 1}; a ∈ Z}, we can define H 1 = H 1 0 ∪ H 1 1 and the set H n as the n-ary Cartesian power of H 1 . If an element h ⊂ Z n is the Cartesian product of d elements of H 1 1 and n -d elements of H 1 0 , we say that h is a d-face of H n and that d is the dimension of h. The set of all faces, H n , is called the nD space of cubical complexes. Figure 4 depicts a set of faces {f, g, h} ⊂ H 2 where f = {0}×{1}, g = {0, 1}×{0, 1}, and h = {1}×{0, 1}; the dimension of those faces are respectively 0, 2, and 1. Let us write The set of faces of H n is arranged onto a grid, the so-called Khalimsky's grid, depicted in gray in Figure 4 (right); and inclusion between faces lead to a neighborhood relationship, depicted in gray and yellow. The set of 2-faces, the minimal open sets of H n , is the n-Cartesian product of H 1 and is denoted by

h ↑ = {h ∈ H n | h ⊆ h } and h ↓ = {h ∈ H n | h ⊆ h}. The pair (H n , ⊆) forms a poset and the set U = {U ⊆ H n | ∀h ∈ U, h ↑ ⊆ U } is a T0-Alexandroff topology on H n . With E ⊆ H n ,
H n 1 .

Set-valued maps

A set-valued map u : X Y is characterized by its graph, Gra(u) = { (x, y) ∈ X × Y | y ∈ u(x) }.
There are two different ways to define the "inverse" of a subset by a set-valued map:

u ⊕ (M ) = { x ∈ X | u(x) ∩ M = ∅ } is the inverse image of M by u, whereas u (M ) = { x ∈ X | u(x) ⊂ M } is the core of M by u.
Two distinct continuities are defined on set-valued maps. The one we are interested in is the "natural" extension of the continuity of a single-valued function. When X and Y are metric spaces and when u(x) is compact, u is said to be upper semi-continuous (u.s.c.) at

x if ∀ε > 0, ∃ η > 0 such that ∀ x ∈ B X (x, η), u(x ) ⊂ B Y (u(x), ε)
, where B X (x, η) denotes the ball of X of radius η centered at x. One characterization of u.s.c. maps is the following: u is u.s.c. if and only if the core of any open subset is open.

Interpolation

Following the conclusions of Section 2.3, we are going to immerse a discrete nD function defined on a cubical grid u : Z n → Z into some larger spaces in order to get some continuity properties. For the domain space, we use the subdivision

X = 1 2 H n of H n . Every element z ∈ Z n is mapped to an element m(z) ∈ 1 2 H n 1 with z = (z 1 , . . . , z n ) -→ m(z) = {z 1 , z 1 + 1 2 }×. . .×{z n , z n + 1 2 }.
The definition domain of u, D ⊆ Z n , has thus a counterpart in X, that will also be denoted D, and that is depicted in bold in Figure 5. For the value space, we immerse Z (the set of pixel values) into the larger space Y = 1 2 H 1 , where every integer becomes a closed singleton of H 1 0 . Thanks to an "interpolation" function, we can now define from u a set-valued map u = I(u). We have u : X Y and we set:

∀ h ∈ X, u(h) =    { u(m -1 (h)) } if h ∈ D max( u(h ) : h ∈ st(cl (h)) ∩ D ) if h ∈ 1 2 H n 1 \D span( u(h ) : h ∈ st(h) ∩ D ) if h ∈ X\ 1 2 H n 1 . (1)
An example of interpolation is given in Figure 5. Actually, whatever u, such a discrete interpolation I(u) can also be interpreted as a non-discrete set-valued map I R (u) : R n R (schematically To the authors knowledge the notion of cuts (or thresholds) have not been defined for set-valued maps. Since they are of prime importance for mathematical morphology, and for the tree of shapes in particular, we propose in this paper the following definitions. Given λ ∈ Y , let us state that [

I R (u)(x) = I(u)(h) with h such as x ∈ R n falls in h ∈ 1 2 H n ),
u λ ] = { x ∈ X | ∀ µ ∈ u(x), µ < λ } and [ u λ ] = { x ∈ X | ∀ µ ∈ u(x), µ > λ }.
We can show [START_REF] Najman | Discrete set-valued continuity and interpolation[END_REF] that, with those definitions, ∀u, ∀λ, [ I(u) λ ] and [ I(u) λ ] are well-composed [START_REF] Latecki | Well-composed sets[END_REF]. That is, strict cut components and their complementary sets can be handled both with the same unique connectivity, c 2n . As a consequence, the operators star and Sat commute on those sets, and we can prove [START_REF] Najman | Discrete set-valued continuity and interpolation[END_REF] that:

S I (u) = { Sat c2n (Γ ); Γ ∈ {CC c2n ([ I(u) λ ])} λ ∪ {CC c2n ([ I(u) λ ])} λ }
is a set of components that forms a tree. Moreover, we can also prove that

T < (u) = { Γ ∩ D; Γ ∈ { CC c2n ([ I(u) λ + 1/2 ]) } λ∈H0 }

Putting things altogether

About saturation and initialization

Classically the root node of the tree of shapes represents the whole image and, formally, the saturation operator is defined w.r.t. a point at infinity, p ∞ , located outside the image domain D. A rather natural idea is that the root level, ∞ , should only depend on the internal border of D (which is unfortunately not the case for the algorithms proposed in the literature). To that aim, before interpolating u, we add to this image an external border with a unique value, ∞ , set to the median value of the internal border. p ∞ is then one point from the added border.

Handling the hierarchical queue

To sort the faces of the domain X of U , we use a classical front propagation based on a hierarchical queue [START_REF] Meyer | Un algorithme optimal de ligne de partage des eaux[END_REF], denoted by q, the current level being denoted In [START_REF] Caselles | Geometric Description of Images as Topographic Maps[END_REF] the "Fast Level Set Transform" (flst) relies on a region-growing approach to decompose the image into shapes. It extracts each branch of the tree starting from the leaves and growing them up to the root until at least one saddle point is encountered. Each time a saddle point is encountered, the branch extraction procedure has to stop until every parallel branch meeting at this point is extracted. So each saddle point invalidates the shape currently being extracted, forcing the algorithm to visit its pixels again once a parallel branch is extracted. Since an image like a checkerboard contains O(n) saddle points meeting on O(n) pixels, the flst has a O(n 2 ) worst case time complexity.

Song [START_REF] Song | A topdown algorithm for computation of level line trees[END_REF] takes a top-down approach to build the tree of level lines in O(n+t) time, where t is the total length of all level lines (note that filling the interior of each level line allows for retrieving the tree of shapes). The algorithm is restricted to 2D images with hexagonal pixels. Its key idea is to perform a recursion (starting from the image boundary): for a given component, follow every contours of its holes, and repeat this procedure for each hole component. Since the total length of level lines of an image can be of order O(n 2 ), the worst case has a quadratic-time complexity.

Conclusion

In this paper, we have presented a new algorithm to compute the tree of shapes of an image which features a quasi-linear time complexity, runs on nD images, and benefits from a much simpler implementation than existing algorithms. We have also proposed a novel representation of images as set-valued maps which has some continuity properties while remaining discrete.

Actually we believe that this representation is a good start to get a "pure" self-duality for images and operators, that is, a way to get rid of the pair of dual connectivities c 2n and c 3 n -1 , and of the dissymmetry of cuts (strict and large cuts for respectively lower and upper cuts). In particular, replacing the maximum operator by the median operator in Equation 1 leads to a pure self-dual definition of the tree of shapes of 2D images [START_REF] Najman | Discrete set-valued continuity and interpolation[END_REF]. Furthermore the perspectives offered by that new representation might be far from being limited to the tree of shapes computation.

For our experiments we use our free software library [START_REF] Levillain | Why and how to design a generic and efficient image processing framework: The case of the Milena library[END_REF]; in particular, the fact that our tool makes it easy to write generic software in the case of mathematical morphology and discrete topology is discussed in [START_REF] Levillain | Writing reusable digital geometry algorithms in a generic image processing framework[END_REF]. The work presented here will be available in the next release of our software for we advocate reproducible research.

Fig. 1 .

 1 Fig. 1. Three morphological trees of the same image.

Fig. 3 .

 3 Fig. 3. A sample image and its tree of shapes (left); a step towards an ad-hoc image representation (right).

  we have a star operator st(E) = ∪ h∈E h ↑ and a closure operator cl (E) = ∪ h∈E h ↓ , that respectively gives the smallest open set and the smallest closed set of P(H n ) containing E.

Fig. 4 .

 4 Fig.[START_REF] Caselles | Constructing the tree of shapes of an image by fusion of the trees of connected components of upper and lower level sets[END_REF]. Three faces depicted as subsets of Z 2 (left) and as geometrical objects (middle); Khalimsky grid (right) with 0-to 2-faces respectively painted in red, blue, and green.

Fig. 5 .

 5 Fig. 5. The function u : Z 2 → Z (left) is transformed into the set-valued map u :

1 2 H 1

 11 (middle); d-faces with d ∈ {0, 1} are interval-valued in u with the span of their respective (d + 1)-face neighbors (right).

  and T ≥ (u) = { Γ ∩ D; Γ ∈ { CC c2n ([ I(u) λ + 1/2 ]) } λ∈H0 }. So eventually we have: S(u) = { Γ ∩ D; Γ ∈ S I (u)}. That final property means that strict cuts of the interpolation of u, considering only c 2n for the different operators, allows for retrieving the shapes of u, as defined with the pair of dual connectivities c 2n and c 3 n -1 .

  and we can show that I R (u) is an u.s.c. map.

						24	24	24	24	24	24	24	24	24	24	24
						24	24	24	24 24 24 24 24 24 24	24
						24	24	24	24	0	0	0	0	0	24	24	8	8	8	[0,8]	0
						24	24 24 24	6	8	8 8	0	24	24
	0 24 24 24 24 24 24 0 0 24 24 24	0 24 24 24 24 24 24	6 24	6 6	8 8	8 8	0 24 24 24 8 24	24 24	8 8	[8,24] [8,24]	24 [8,24]	24 [0,24]	24 [0,24]
	24	0	6	8 0	24	24	24	24	24	0	0	0	24	24	24	24	[0,8]	[0,24]	24	24	24
	0 24 24 24 24 24 24 24 0 24 24 24	24 24 24 24 24 24	24 24 24 24 24 24 24 24 24 24	24 24 24 24 24 24	0	[0,24]	24	24	24

In 2D, a cavity of a set S ∈ Ω is called a "hole"; in nD, it is a connected component of Ω\S which is not the "exterior" of S. Browsing the elements of S in nD, with n ≥ 3, does not allow to know whether S has a cavity or not[START_REF] Henle | A Combinatorial Introduction to Topology[END_REF].

Due to limited place, this paper does not contain the following topics (they will be included into an extended version of this paper). A comparison of execution times of existing algorithms. Actually it is possible to reduce the space complexity (i.e., memory usage) of the algorithm proposed in this paper so the shorter version presented here is not our "competitive" version. The union-by-rank procedure that guaranties quasi-linear complexity. So that the Union-Find routine (given in[START_REF] Berger | Effective component tree computation with application to pattern recognition in astronomical imaging[END_REF] and recalled in Algorithm 1) remains short, its code does not feature tree balancing; yet it is explained in[START_REF] Carlinet | A (fair?) comparison of many max-tree computation algorithms[END_REF]. A formal proof of our algorithm. This paper focuses on how the proposed algorithm works and gives an insight into the reasons why it works; to give a formal proof requires a large amount of materials, the first part of which can be found in[START_REF] Najman | Discrete set-valued continuity and interpolation[END_REF]. About high bit-depths data. That case is not detailled in this paper.

In[START_REF] Najman | Discrete set-valued continuity and interpolation[END_REF], a formal characterization of the discrete topology underlying this novel representation is presented.

[START_REF] Géraud | Mathematical Morphology-From Theory to Applications[END_REF] The authors recommend[START_REF] Mazo | Digital imaging: A unified topological framework[END_REF] and[START_REF] Aubin | Set-Valued Analysis[END_REF] for extra readings about those tools.

Formally the sorting step has the pseudo-polynomial O(k n) complexity, k being the number of different gray values. Though, since we consider low bit-depths data, k shall only be considered as a complexity multiplicative factor.

Acknowledgements. The authors would like to thanks Michel Couprie and Jean Cousty for fruitful discussions. This work received funding from the Agence Nationale de la Recherche, contract ANR-2010-BLAN-0205-03 and through "Programme d'Investissements d'Avenir" (LabEx BEZOUT n • ANR-10-LABX-58).

Algorithm 2: Sorting for tree of shapes computation.

priority push(q, h, U, ) /* modify q */ begin [lower , upper ] ← U (h) if lower > then ← lower else if upper < then ← upper else ← push(q[ ], h) priority pop(q, ) : H /* modify q, and sometimes */ begin if q[ ] is empty then ← level next to such as q[ ] is not empty ← return pop(q[ ]) sort(U ) : Pair(Array[H], Image) begin for all h do deja vu(h)

by . The sorting algorithm is given in Algorithm 2. There are two notable differences with the well-known hierarchical-queue-based propagation. First the d-faces, with d < n, are interval-valued so we have to decide at which (singlevalued) level to enqueue those elements. The solution is straightforward: a face h is enqueued at the value of the interval U (h) that is the closest to (see the procedure priority push). Just also note that we memorize the enqueuing level of faces thanks to the image u (see the procedure sort). Second, when the queue at current level, q[ ], is empty (and when the hierarchical queue q is not yet empty), we shall decide what the next level to be processed is. We have the choice of taking the next level, either less or greater than , such that the queue at that level is not empty (see the procedure priority pop). Practically choosing going up or down the levels does not change the resulting tree since it just means exploring some sub-tree before some other disjoint sub-tree. The result R of the sorting step is the one expected since the image U , in addition with the browsing of level in the hierarchical queue, allows for a propagation that is "continuous" both in domain space and in level space. An interesting property due to the interpolation and the well-composedness of cuts is that the neighborhood N , used for faces in the propagation, corresponds to the c 2n connectivity on the Khalimsky's grid.

Max-tree versus tree of shapes computation

The main body of the tree of shapes computation algorithm is given in Algorithm 3. The major differences between this algorithm and the one dedicated to the max-tree (see the procedure compute tree in Algorithm 1) are the following ones. First the three basic steps (sort, union-find, and canonicalization) are now surrounded by an interpolation and un-interpolation process. Note that the uninterpolation just cleans up both R and parent to keep only elements of D. Second, as emphasized in Section 2.2, the sorting step is of course dedicated to the tree of shapes computation. Last, a temporary image, u , is introduced. It is defined on the same domain as u, namely X, and contains only single-valued elements. This image is the equivalent of the original image u when dealing with the max-tree: it is used to know when an element h is canonical, that is, when u (parent(h)) = u (h) (so that image is thus required by the canonicalization step that runs on X).

Complexity analysis of the algorithm presented here is trivial. The interpolation, canonicalization, and un-interpolation are linear. The modified union-find (once augmented with tree balancing, i.e., union-by-rank) is quasi-linear when values of the input image u have a low quantization (typically 12 bits or less). Last, the time complexity of the sorting step is governed by the use hierarchical queue: it is linear with low quantized data 8 . Eventually we obtain a quasi-linear algorithm. The representation of the tree with the pair (R, parent) allows for any manipulation and processing that one expects from a morphological tree [START_REF] Carlinet | A (fair?) comparison of many max-tree computation algorithms[END_REF].

Related works

The first known algorithm, the "Fast Level Line Transform (fllt)" [START_REF] Monasse | Fast computation of a contrast invariant image representation[END_REF], computes the max-tree and the min-tree of an image and obtains the tree of shapes by merging both trees. The main drawback of the fllt is the need to know that a component has an hole (in order to match it with a component of the other tree). To that aim the Euler characteristic is computed, which can be done locally (while following the border of components) but in 2D only. In [START_REF] Caselles | Constructing the tree of shapes of an image by fusion of the trees of connected components of upper and lower level sets[END_REF][START_REF] Meinhardt-Llopis | Morphological and Statistical Techniques for the Analysis of 3D Images[END_REF] the authors show that this fusion approach is sound in nD with n > 2; yet it cannot be effective in practice due to unacceptable complexity.


Laurent Najman

Jean Cousty

Benjamin Perret

Playing with Kruskal: algorithms for morphological trees in edge-weighted graphs

The goal of this paper is to provide linear or quasi-linear algorithms for producing some of the various trees used in mathemetical morphology, in particular the trees corresponding to hierarchies of watershed cuts and hierarchies of constrained connectivity. A specific binary tree, corresponding to an ordered version of the edges of the minimum spanning tree, is the key structure in this study, and is computed thanks to variations around Kruskal algorithm for minimum spanning tree.

Introduction

In the theoretical companion paper [START_REF] Cousty | Constructive links between some morphological hierarchies on edge-weighted graphs[END_REF] of the present paper, we show how some morphological hierarchies [START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF][START_REF] Nagao | Region extraction and shape analysis in aerial photographs[END_REF][START_REF] Meyer | Morphological scale-space representation with levelings[END_REF][START_REF] Salembier | Binary partition tree as an efficient representation for image processing, segmentation, and information retrieval[END_REF][START_REF] Najman | Geodesic saliency of watershed contours and hierarchical segmentation[END_REF][START_REF] Morris | Graph theory for image analysis: an approach based on the shortest spanning tree[END_REF] defined on an edge-weighted graph G = (V, E, F ) are related, and when it is possible, how they can be computed one from each other. In this paper, we provide efficient quasi-linear or linear algorithms to compute those hierarchies. In particular, this paper contains:

provided that the edges are either already sorted or can be sorted in linear time, a quasi-linear O(|E| × α(|V |)) (where α() is the extremely slowly growing inverse of the single-valued Ackermann function) algorithm that computes a binary partition tree by altitude ordering, a fundamental structure that we post-process in the sequel; -a linear O(|V |) post-processing algorithm that computes the hierarchy of quasi-flat zones [START_REF] Nagao | Region extraction and shape analysis in aerial photographs[END_REF][START_REF] Cousty | Constructive links between some morphological hierarchies on edge-weighted graphs[END_REF] (also know as the α-tree [START_REF] Ouzounis | Pattern spectra from partition pyramids and hierarchies[END_REF]); -a linear O(|V |) post-processing algorithm that computes (hierarchies of) watershed cuts [START_REF] Cousty | Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle[END_REF]; -a linear O(|V |) post-processing algorithm that computes hierarchies by increasing attributes; as detailed here and in [START_REF] Cousty | Constructive links between some morphological hierarchies on edge-weighted graphs[END_REF], such an algorithm can be used to obtain either constrained connectivity hierarchies [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF] or watershed-based hierarchies [START_REF] Meyer | Segmentation, minimum spanning tree and hierarchies[END_REF].

To the best of our knowledge, the only published constrained connectivity algorithm, available in [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF], has an unknown complexity and only computes one level of the hierarchy. An algorithm computing the whole hierarchy, relying on the component tree of the edge-weighted graph G, is roughly sketched in [START_REF] Najman | On the equivalence between hierarchical segmentations and ultrametric watersheds[END_REF], but has a complexity higher than the one proposed in this paper and is less memory efficient. Concerning attribute-based hierarchies, the most efficient algorithm [START_REF] Cousty | Incremental algorithm for hierarchical minimum spanning forests and saliency of watershed cuts[END_REF] has a complexity higher than the one proposed in this paper, and is less efficient.

At the heart of our approach is the minimum spanning tree (MST): this tree T is a connected spanning graph of the graph G such that the weight of T : F (T ) := e∈E(T ) F (e) is the least possible weight for a spanning graph of G. As detailed in section 2.1, we rely on Kruskal algorithm [START_REF] Kruskal | On the shortest spanning subtree of a graph and the traveling salesman problem[END_REF] for computing this MST. However, while producing the MST, we make use of another tree (detailed in section 2.2) that we call the binary partition tree by altitude ordering. Using Tarjan union-find (section 2.3), we propose in section 2.4 an efficient algorithm to compute this binary tree. Post-processing will be studied in section 3.

Binary Partition Tree and Minimum Spanning Tree

Kruskal algorithm

Kruskal's algorithm [START_REF] Kruskal | On the shortest spanning subtree of a graph and the traveling salesman problem[END_REF] is a greedy algorithm that finds a minimum spanning tree for a connected weighted graph. It can be described as follows:

create a forest F (a set of trees), where each vertex in the graph is a separate tree create a set S containing all the edges in the graph while S is nonempty and F is not yet a single tree

• remove an edge with minimum weight from S • if that edge connects two different trees, then add it to the forest, combining two trees into a single tree • otherwise discard that edge.

At the termination of the algorithm, the forest has only one component and forms a minimum spanning tree of the graph.

The efficiency of Kruskal's algorithm relies on a disjoint-set data structure that keeps track of a set of elements partitioned into a number of disjoint (nonoverlapping) subsets. The disjoint set problem consists in maintaining a collection Q of disjoint sets under the operation of union. Each set Q in Q is represented by a unique element of Q, called the canonical element. In the following, q 1 and q 2 denote two distinct elements. Three operations allow to manage the collection:

-MakeSet(q 1 ): add a new element q 1 to the collection Q, provided that the element q 1 does not already belongs to a set in Q. -FindCanonical(q 1 ): return the canonical element of the set in Q which contains q 1 . -Union(q 1 , q 2 ): let Q 1 and Q 2 be the two sets in Q whose canonical elements are q 1 and q 2 respectively (q 1 and q 2 must be different). Both sets are removed from Q, their union

Q 3 = Q 1 ∪ Q 2 is added to Q and a canonical element for Q 3 is selected and returned.
An implementation of Kruskal algorithm is presented in Algorithm 1. In this implementation, we identify any element of V with an integer corresponding to its index in the finite set |V |. We save the edge of the MST in an array MST, hence obtaining a strict order on the edges; this order is necessary for post-processing. When computing an MST, the DoSomething procedure does nothing, i.e., it discards the considered edge. In section 3.3, we show an example when the procedure DoSomething is useful.

Procedure DoSomethingMST({x, y})

// Ignore {x, y}
The main question in implementing Kruskal's algorithm is thus how to represent and implement the collection Q. A good representation is to maintain the collection as a set of trees, i.e each element of Q is a tree. In the sequel of this paper, each set of the collection Q is represented by a rooted tree, where the canonical element of the set is the root of the tree. We are going to play with various tree representations of connected components. Kruskal's algorithm will not change, but different implementations for Union and FindCanonical will lead to two different trees, one of them being useful for connected filtering.

A simple algorithm that yields a binary partition tree by altitude ordering

The first tree we present is the useful one, although the proposed algorithm in this section is not an efficient one (we make a better proposal in section 2.4). As shown in [START_REF] Cousty | Constructive links between some morphological hierarchies on edge-weighted graphs[END_REF] and in the sequel of this paper, post-processing this tree provides the min-tree of MST, the tree of quasi-flat zones and trees of watersheds. The main idea is the following: each time an edge {x, y} is put into the MST, i.e. each time a union is made, we create a new node whose children are the two disjoint sets containing x and y. Intuitively, it is as if we break the edge in two, and add a node between the two points of the edge. The added node becomes the canonical element of the union of these two points (see Fig. corresponding to vertices are from 0 to |V | -1 and nodes corresponding to edges of the MST are from |V | to 2|V | -2. Such a numbering allows the tree itself to be provided thanks to an array parent that stores the parent of a given node. If node i does not have a parent, then parent

[i] := -1. This algorithm is implemented in Q BT .MakeSet, Q BT .FindCanonical and Q BT .Union. Procedure Q BT .MakeSet(q) 1 Q BT .parent[q] := -1; Q BT .size += 1; Function Q BT .FindCanonical(q) 1 while Q BT .parent[q]≥ 0 do q :=Q BT .parent[q]; 2 return q; Function Q BT .Union(c x , c y ) 1 Q BT .parent[c x ]:=Q BT .size; Q BT .parent[c y ]:=Q BT .size; 2 Q BT .MakeSet(Q BT .size); 3 return Q BT .size-1 ;
At the end, we obtain a collection Q BT which is a binary tree. The Q BT tree provides the order in which the edges have been put into the MST, the latest edge added to the tree being the root of the tree, i.e., the highest one. In other words, the edges of the MST are strictly ordered by Q BT , according to their altitude in the tree. We say that Q BT is a binary partition tree by altitude ordering [START_REF] Cousty | Constructive links between some morphological hierarchies on edge-weighted graphs[END_REF]. The complexity of Kruskal implemented with this specific union-find is in O(|V | 2 ), and thus, this process is not very efficient.

Efficient MST implementation with Tarjan Union-Find

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 union and find. The implementation of this algorithm is given in procedure Q T .MakeSet and functions Q T .Union and Q T .FindCanonical. To each element of the collection is associated a parent (as precedently) and a rank 'Rnk'. Both the mapping 'parent' and the mapping 'Rnk' are represented as 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 FindCanonical. It consists, while searching for the root r of the tree which contains q, in considering each element p of the path from q to r (including q), and setting the parent of p 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 Union operation. The rank Rnk(c x ) of a canonical element c x is a measure of the depth of the tree rooted in c x , and is exactly the depth of this tree if the path compression technique is not used jointly with the union by rank technique. If the two canonical elements c x and c y have the same rank, then one of the elements, say c y , is chosen arbitrarily to be the canonical element of the union: c y becomes the parent of c x ; and the rank of c y is incremented by one. 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 Q T .MakeSet(q) 1 Q T .parent[Q T .size] := -1; Q T .Rnk[Q T .size] := 0; Q T .size += 1; Function Q T .FindCanonical(q) 1 r := q; 2 while Q T .parent[r] ≥ 0 do r :=Q T .parent[r]; 3 while Q T .parent[q] ≥ 0 do tmp := q; q :=QT .parent[q]; QT .parent[tmp] := r; 4 return r; Function Q T .Union(c x , c y ) 1 if (Q T .Rnk[c x ] >Q T .Rnk[c y ]) then swap(c x , c y ); 2 if (Q T .Rnk[c x ] == Q T .Rnk[c y ]) then Q T .Rnk[c y ] += 1; 3 Q T .parent[c x ] := c y ; 4 return c y ;
As stated at the beginning of this section, implementing Kruskal's algorithm with Tarjan Union-Find leads to a quasi-linear complexity: provided that the edges are either already sorted or can be sorted in linear time (for example with counting sort or radix sort), the complexity is O(|E| × α(|V |)), where α() is the extremely slowly growing inverse of the single-valued Ackermann function. Unfortunately, the tree built by Tarjan Union-Find is of no use for connected filtering: the path compression technique flattens the tree, and does not preserve the order by which the edges or the vertices are processed in the algorithm. In the next section, we are going to combine the two proposals into an efficient one.

An efficient algorithm that yields a binary partition tree by altitude ordering

The function Q BT .FindCanonical is slow because it takes some time to find the canonical element for a connected component: we have to climb the tree until the root is found. We can use Tarjan Union-Find on Q BT itself, i.e. we can use Q T as a second collection that maintains a compressed representation of Q BT , so that finding the canonical element is now in quasi-constant time. We also have to store the root of any tree in Q BT in an array Root so that such a root can be found in constant time from any tree in Q T . The implementation uses Q EBT .MakeSet, Q T .FindCanonical and Q EBT .Union.

Procedure Q EBT .MakeSet(q) 1 Q EBT .Root[q]:=q; Q BT .MakeSet(q); Q T .MakeSet(q); Function Q EBT .Union(c x , c y ) 1 t u :=Q EBT .Root[c x ]; t v := Q EBT .Root[c y ]; // Union in Q BT (without compression) 2 Q BT .parent[t u ] := Q BT .parent[t v ] := Q BT .size; // If children are needed, add them to the root 3 Q BT .children[Q BT .size].add({t u }); Q BT .children[Q BT .size].add({t v }); 4 c:=Q T .Union(c x ,c y ); // Union in Q T (with compression) 5 Q EBT .Root[c] := Q BT .size; // Update the root of Q EBT 6 Q BT .MakeSet(Q BT .size); 7 return Q BT .size-1 ; Function Q EBT .FindCanonical(q) 1 return Q T .FindCanonical(q);
When Kruskal is finished, the tree Q BT is exactly the same as the one in section 2.2. The only difference is thus the speed of the algorithm: thanks to the use of Tarjan Union-Find, the complexity of this Kruskal algorithm using Q EBT is quasi-linear O(|E| × α(|V |)).

Post-Processing the binary tree

In this section, we are going to detail some linear O(|V |) algorithms that produce, from Q BT , (1) a watershed cut, (2) a hierarchy of quasi-flat zones, and (3) any attribute-based hierarchy (if the attribute is increasing).

As we have seen, each node of Q BT corresponds either to a vertex of the graph or to an edge of the MST. Recall that edges of MST are sorted by a strict order relationship that follows increasing weight-edges: the |V | first nodes of Q BT are the vertices of |V |, and the root of Q BT corresponds to the edge of the MST with the greatest weight. To ease the reading of the algorithms of this section, we provide below two helper functions that clarify how we can pass from the nodes of Q BT to the edges of the MST and how to obtain the weight of the edge of the MST corresponding to a given node of Q BT .

Function getEdge(n)

Data: a (non-leaf) node n of Q BT Result: the edge e of the MST corresponding to the n th node 1 return n -|V |;

Function weightNode(n)

Data: a (non-leaf) node of the tree Result: the weight of the MST edge associated with the n th node of Q BT 1 return F(MST[getEdge(n)]);

Watershed cuts

A watershed cut [START_REF] Cousty | Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle[END_REF] can be easily obtained in our framework. Indeed, those cuts are the "highest separations" of minima of the minimum spanning tree. We propose below a linear algorithm watershed that labels the edges of the MST with a flag stating whether or not an edge is a watershed edge. The main idea is to obtain the latest edge of the MST included in a given minimum in Q BT . This can be done by counting the number of minima of the edge-weighted graph (G, F ) thanks to a simple process than run through the nodes of Q BT by increasing order, and checks whether or not a given node has an altitude lower than its parent and does not contain a minimum. We increment a counter for the ancestors of several minima. An edge is a watershed edge if it merges several catchment basins corresponding to different minima. The complexity of the function watershed is linear in the number of vertices, and thus the whole process that computes a watershed cut is quasi-linear. The set of watershed edges provides a MST of the neighborhood graph of the catchment basins [START_REF] Meyer | Minimum spanning forests for morphological segmentation[END_REF]. In Fig. 2.a, the nodes of Q BT corresponding to watershed edges have a blue point.

By removing from Q BT all nodes that are not watershed ones, we obtain a filtered tree that corresponds to a hierarchy of watershed cuts, more precisely the one corresponding to an ultrametric watershed [START_REF] Najman | On the equivalence between hierarchical segmentations and ultrametric watersheds[END_REF]. In the sequel of this paper, one can use indifferently either the original Q BT or this modified one. In the first case, one is working within the framework of the constrained connectivity, and in the second case, one is working within the framework of watershed cuts.

One of the interests of working with watershed cuts rather than with flat zones is that the hierarchy is smaller, as the super-pixels (catchment basins) provided by watershed cuts are larger than the super-pixels provided by the flat-zones. From a practical point of view, greater speed can thus be achieved thanks to watershed cuts.

Quasi-flat zones hierarchy

In this section, we post-process Q BT to obtain the tree of the quasi-flat zones hierarchy [START_REF] Nagao | Region extraction and shape analysis in aerial photographs[END_REF], which is proved [START_REF] Cousty | Constructive links between some morphological hierarchies on edge-weighted graphs[END_REF] to be Q CT , the min-tree of the minimum spanning tree. The differences between Q BT and Q CT are illustrated in Fig. 1.d and 2.b. Q CT can be computed directly thanks to a dedicated Union-Find procedure that we will describe in an extended version of this paper. Here, due to space constraints, we only propose a short post-processing that transforms Q BT into the desired min-tree Q CT . The implementation uses CanonizeQ BT to post-process Q BT , and needs the children of a node.

Procedure CanonizeQ BT Data: Q BT Result: Q CT , a canonized version of Q BT 1 for all nodes n of Q BT do Q CT .parent[n]:=Q BT .parent[n]; Q CT .size+=1;
2 for each non-leaf and non-root node n of Q BT by decreasing order do The procedure CanonizeQ BT is in O(|V |), and thus the whole process that computes a quasi-flat zones hierarchy is quasi-linear. To the best of our knowledge, this is the most efficient algorithm published for computing this hierarchy. However, for most image-processing tasks, the binary partition tree Q BT can be used instead, and this is what we are going to do in the sequel of this paper.

Attribute-based hierarchies

Attributes It is easy to compute some attributes on each node of Q BT : surface (number of vertices in a node), depth, volume or ordered markers are the most classical attributes [START_REF] Meyer | Segmentation, minimum spanning tree and hierarchies[END_REF]. Another attribute from the constrained connectivity framework [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF] is the range. Any of those previous attributes are increasing: recall that an attribute A is increasing if when there is a parenthood relationship between two nodes n 1 and n 2 , i.e. when the vertices of n 1 are contained in the vertices of n 2 , then A(n 1 ) ≤ A(n 2 ). In the sequel of this section, we denote by attributeComp[n] the attribute of the node n of Q BT .

It is when computing an attribute that the procedure DoSomethingMST can be useful: for example, one can, using DoSomethingMST, register each edge of the graph at the correct place in Q BT , by adding new nodes corresponding to these edges. Then, for example, we can imagine using an attribute such as the surface computed not from the vertices but from the edges.

Hierarchies Attribute-based hierarchies are obtained by filtering the min-tree Q CT of the MST, and computing the full hierarchy amounts to reweighting the MST. Intuitively, it is as if colored water fills up the branches of Q CT , a merging of two components taking place when two different colors meet. The speed of the filling is controled by the attribute, i.e., the branch is completely filled (i.e., the branch is cut) when the amount of water is exactly equal to the attribute of the corresponding node of Q CT . Some nodes exist in Q BT but not in Q CT . Those nodes correspond to an information not present in Q CT : the order of processing, useful when a choice has to be made. A complete illustrative example explaining the process is detailed in Fig. 3. Computing the correct attribute is done thanks to the function getAttribute, that computes from attributeComp the attribute at the time of the merging. It consists in taking the highest attribute of all the children with the same original weight. The final weight of the MST can then be obtained by taking the lowest attribute of the two children of the node corresponding to the edge of the MST in Q BT ; this is done by ComputeMergeAttributeMST, the complexity of the whole post-processing being in O(|V |). Once the re-weighted MST is computed, the hierarchical tree can be obtained by re-applying the algorithms of this paper. As seen in [START_REF] Cousty | Constructive links between some morphological hierarchies on edge-weighted graphs[END_REF], the resulting hierarchy can be either a hierarchy of watershed cuts or a constrainedconnectivity hierarchy. . When a flooding-by-attribute is performed on QCT , node n1 disappears first at value 2, followed by node n2 at value 3. When these the two corresponding branches are filled by water (cut from the tree), n1 and n2 are no longer minima. This is not the case for n3, whose attribute is 6. When the corresponding branch is filled, n3 stills marks a minimum, and stays so until the value 11, corresponding to node n7. Thanks to the processing order embodied in QBT , getAttribute can compute the value at which a node disappears.

Function getAttribute(n)

Procedure ComputeMergeAttributeMST

Data: QBT Result: a reweighted MST G corresponding to the attribute-based hierarchy 1 for any non-leaf node n of QBT do

Conclusion

This paper has presented several elegant yet efficient algorithms for computing several morphological trees. At the heart of the processing is the minimum spanning tree, and in this paper we have proposed some variations on Kruskal algorithm. However, other approaches can be taken, and any other MST algorithm can be used first to produce a tree on which the algorithms of this paper can be applied. This would be needed in some situations, for example if the edges of the original graph do not fit in memory or if some parallel algorithm for minimum spanning tree is first needed. The unification theory provided in [START_REF] Cousty | Constructive links between some morphological hierarchies on edge-weighted graphs[END_REF], together with the algorithms of this paper shed a new light on what has been done in mathematical morphology for a number of years, linking together some previously unrelated parts of the field.

Source code corresponding to this paper is available at http://www.esiee. fr/ ~info/sm/.

Algorithm 1 : 5 if cx = cy then 6 Q

 156 KruskalData: An edge-weighted graph (V, E, F ). Result: A minimum spanning tree MST Result: A collection Q // Collection Q is initialized to ∅ 1 e := 0; 2 for all xi ∈ V do MakeSet(i); 3 for all edges {x, y} by (strict) increasing weight F ({x, y}) do 4 cx := Q.FindCanonical(x); cy := Q.FindCanonical(y); .Union(cx, cy); 7 MST[e] := {x, y};

8 e := e + 1 ; 9 else

 819 DoSomething({x, y});

Fig. 1 .

 1 Fig.1. A simple process for obtaining QBT , a binary tree providing a strict total order relation on the edges of the MST (see text).

Fig. 2 . 12 p

 212 Fig. 2. Two trees derived from QBT by post-processing. (a) A hierarchical tree QBT on which watershed edges have a blue point. Edges M1 (green), M2 (red) and M3 (blue) are the minima of the edge-weighted graph. (b) The tree of the quasi-flat zones hierarchy, a canonized version of QBT (i.e., the min-tree of the MST.)

3 p 6 Q

 36 := Q CT .parent[n]; 4 if (weightNode(p) == weightNode(n)) then 5 for all c ∈ Q BT .children[n] do Q CT .parent[c]:=p; CT .parent[n]:=n; // Delete node n of Q CT // If needed, build the list of children 7 for all nodes n of Q CT do 8 p:=Q CT .parent[n]; if p ≥ 0 and p = n then Q CT .children[p].add(n);

6 for

 6 Data: A node n of QBT Result: The attribute at the time of the merging 1 if (n is the root) or (weightNode(parent[n]) = weightNode(n)) then 2 for all c children of n do getAttribute(c); 3 attribute[n] := attributeComp[n]; all children c of n do

Fig. 3 .

 3 Fig.3. Hierarchical trees QBT and QCT . At the bottom of (a) and (b) is the (edgeweighted) graph of a map with the respective weight of each edge. Attributes (in this example, the surface of a node) can easily be computed on either QBT or QCT . In (a), the numbers (k, l) in parenthesis represents respectively k=attributeComp[n] and l=getAttribute(n). In (a) and (b), the numbers k in ni : k represents k=attributeComp[n]. When a flooding-by-attribute is performed on QCT , node n1 disappears first at value 2, followed by node n2 at value 3. When these the two corresponding branches are filled by water (cut from the tree), n1 and n2 are no longer minima. This is not the case for n3, whose attribute is 6. When the corresponding branch is filled, n3 stills marks a minimum, and stays so until the value 11, corresponding to node n7. Thanks to the processing order embodied in QBT , getAttribute can compute the value at which a node disappears.

2 a1: 4 G

 24 = attribute[children[n].lef t]; 3 a2 := attribute[children[n].right]; [getEdge(n)] := min(a1, a2);

Acknowledgements. This work received funding from the Agence Nationale de la Recherche, contract ANR-2010-BLAN-0205-03 and through "Programme d'Investissements d'Avenir" (LabEx BEZOUT n • ANR-10-LABX-58).


Yongchao Xu
email: yongchao.xu@lrde.epita.fr

Thierry Géraud
email: thierry.geraud@lrde.epita.fr

Laurent Najman
email: l.najman@esiee.fr

Two applications of shape-based morphology: blood vessels segmentation and a generalization of constrained connectivity

, is a generalization of the existing tree-based connected operators. In this paper, two different applications are studied: in the first one, we apply our framework to blood vessels segmentation in retinal images. In the second one, we propose an extension of constrained connectivity. In both cases, quantitative evaluations demonstrate that shape-based filtering, a mere filtering step that we compare to more evolved processings, achieves state-of-the-art results.

Introduction

Mathematical morphology, as originally developed by Matheron and Serra [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF], proposes a set of morphological operators based on structuring elements. Later, Salembier and Serra [START_REF] Salembier | Flat zones filtering, connected operators and filters by reconstruction[END_REF], followed by Breen and Jones [START_REF] Breen | Attribute openings, thinnings, and granulometries[END_REF], proposed morphological operators based on attributes, rather than on elements. Such operators rely on transforming the image into an equivalent representation, generally a tree of components (e.g., level sets) of the image; such trees are equivalent to the original image in the sense that the image can be reconstructed from its associated tree. Filtering then involves the design of a shape attribute that weights how much a node of the tree fits a given shape. Two different approaches for filtering the tree (and hence the image) have been proposed: the more evolved approach consists in pruning the tree by removing whole branches of the tree, and is easy to apply if the attribute is increasing on the tree (i.e., if the attribute is always stronger for the ancestors of a node). This process is depicted in the black path in Fig. 1. However, most shape attributes are not increasing. When the attribute is not increasing, three strategies have been proposed (min, max, Viterbi; see [START_REF] Salembier | Connected operators[END_REF] for more details). They all choose a particular node on which to take the decision, and remove the whole subtree rooted at this node. While it may give interesting results in some cases, it does not take into account the possibility that several relevant objects can have some inclusion relationship, meaning that they are on the same branch of the tree (e.g., a ring object in a tree of shapes, see Fig. 3.a).

Image f

In the simplest approach, one simply removes the nodes of the tree for which the attribute is lower than a given threshold [START_REF] Urbach | Connected shape-size pattern spectra for rotation and scale-invariant classification of gray-scale images[END_REF]. Such a thresholding does not take into account the intrinsic parenthood relationship of the tree. Moreover it is often impossible to retrieve all expected objects with one unique threshold. Fig. 2 shows the evolution of a shape attribute, the circularity, along two branches of the tree of shapes [START_REF] Monasse | Fast computation of a contrast-invariant image representation[END_REF]. The light round shape and the dark one are both meaningful round objects compared to their context. However, their attribute values are very different. In order to obtain the light one, a higher threshold is applied, but some non-desired shapes appear, the ones in the background in Fig. 2.f.

The founding idea of shape-based morphology is to apply connected filters on the space of all the components of the image, such space being structured into a graph by the parenthood relationship (i.e., the neighbors of a node are its children and its parent). This process is illustrated by the black+red path of Fig. 1. This surprising and simple idea has several deep consequences that were first exposed in [START_REF] Xu | Morphological Filtering in Shape Spaces: Applications using Tree-Based Image Representations[END_REF], where it is shown that this framework encompasses the usual attribute filtering operators [START_REF] Xu | Morphological Filtering in Shape Spaces: Applications using Tree-Based Image Representations[END_REF]. Novel connected filters based on non-increasing criterion can also be proposed. When the first tree T is respectively a Min-tree or a Max-tree, such filters are new morphological lower or upper levelings. When the first tree T is the tree of shapes, we introduce a novel family of self-dual connected filters that we call morphological shapings. The proposed framework can also be used to produce extinction-based [START_REF] Vachier | Extinction values: A new measurement of persistence[END_REF] saliency maps [START_REF] Najman | Geodesic saliency of watershed contours and hierarchical segmentation[END_REF][START_REF] Najman | On the equivalence between hierarchical segmentations and ultrametric watersheds[END_REF]. In general, shape-based morphological operators provide better results than threshold-based approaches. As stated above, when we want to process both upper and lower level sets at the same time, we use as T the tree of shapes, and the such created operator is called a morphological shaping. In both Fig. 2 and Fig. 3, the attribute A is the circularity. The result of the shaping on Fig. 2.a, shown in Fig. 2.g, looks indeed better than the one of Fig. 2.f. In Fig. 3, we compare this extinction-based self-dual shaping approach with a variant of the state-of-the-art thresholding approach [START_REF] Urbach | Connected shape-size pattern spectra for rotation and scale-invariant classification of gray-scale images[END_REF]. When the threshold of A is low, some objects do not appear (Fig. 3.b). To get all the expected objects, we have to set a high threshold; however, in this case, too many unwanted objects are present (Fig. 3.c). With our shaping filter, all the expected objects can be found, as depicted in Fig. 3.d. In this paper, we propose to detail how the framework can be used for two different types of applications. The rest of this paper is organized as follows. An application of our proposed shape-based upper levelings to blood vessels segmentation in retinal images is explained in Section 2. In Section 3, we detail how to extend the constraint connectivity framework first introduced in [START_REF] Soille | Constrained connectivity for hierarchical image decomposition and simplification[END_REF]. Finally we discuss and conclude in Section 4.

Blood vessels segmentation in retinal images

Blood vessels segmentation is a very important task in retinal images analysis. Unlike classical linear filters, connected operators perfectly preserve the location and the shape of the contours, which provides a motivation for using them.

Many existing methods work on the green channel of the color retinal image. To improve the visibility of the blood vessels, for each color retinal image f c , a black top-hat transform is applied to the green channel f g . When a mask of eye fundus is available, we combine it with the black top-hat f t . We thus obtain an image f i in which the blood vessels are visible: indeed, the main structures of the blood vessels are present in the Max-tree T of f i (the connected components of the upper level sets of f i , [START_REF] Salembier | Connected operators[END_REF]). .

For each connected component represented by some node N k of the Max-tree T , we compute a shape attribute A characterizing the blood vessels, which are usually long and thin structures. The attribute used here is the elongation A e :

A e (N k ) = |N k |/(π × l 2 max ), (1)
where |.| denotes the cardinality and l max denotes the length of the largest axis of the best fitting ellipse for the connected component represented by N k . Since the blood vessels are long and thin, nodes having a low attribute A e correspond to the blood vessels. The core of this application is the filtering of the Max-tree T . A mere thresholding of the elongation A e is not sufficient, often giving unwanted objects (noise). However, a very low thresholding value t min on A e ensures that thresholded nodes are blood vessels. Those initial extracted nodes are used as seeds in the sequel. We then apply a morphological filtering with a depth criterion: using the Min-tree T T of the node-weighted graph (T , A e ), we only preserve the nodes that have a certain depth d 0 in T T and that furthermore contain the seeds. The connected components represented by the preserved nodes are considered as the segmented blood vessels. The whole process is one of the many variants of shape-based upper levelings [START_REF] Xu | Morphological Filtering in Shape Spaces: Applications using Tree-Based Image Representations[END_REF].

An example of this blood vessels segmentation methods is given in Fig. 4. As compared with the manually segmented blood vessels segmentation (Fig. 4.f), the elongation-based upper leveling (Fig. 4.e) correctly extracts most of the blood vessels.

We have tested this specific shape-based upper leveling on the Digital Retinal Images for Vessel Extraction (DRIVE) database [START_REF]DRIVE: Digital Retinal Images for Vessel Extraction[END_REF], [START_REF] Staal | Ridge based vessel segmentation in color images of the retina[END_REF] and on the STructured Analysis of the Retina (STARE) database [START_REF] Stare | STructured Analysis of the Retina[END_REF], [START_REF] Hoover | Locating blood vessels in retinal images by piece-wise threshold probing of a matched filter response[END_REF]. DRIVE is a database assembled in the Netherlands from a diabetic retinopathy screening program. It For the test cases, two manual segmentations are available: one is used as gold standard; the other one can be used to compare computer generated segmentations with those of an independent human observer. The complete database contains seven pathological cases (four in the test set and three in the training set). The STARE database contains 20 images captured using a TopCon TRV-50 fundus camera at 35 • FOV, and digitized to 700 × 605 pixels, 8 bits per RGB channel. A manual segmentation is available for each image of the database. Masks of the eye fundus, derived from the matched spatial filter [START_REF] Hoover | Locating blood vessels in retinal images by piece-wise threshold probing of a matched filter response[END_REF], are also available. Note that, among the 20 images, 10 images are abnormal. Fig. 5 and6 show four segmentation results respectively from the DRIVE database and the STARE database. Qualitatively, most of the blood vessels are correctly extracted, although some noise points at the end of the vessels are also kept, and some very thin blood vessels are missed.

Quantitative assesment is based on three performance measurements named respectively sensitivity, specificity and accuracy [START_REF] Staal | Ridge based vessel segmentation in color images of the retina[END_REF]. Sensitivity measures the true positive rate (TPR), specificity measures the true negative rate (TNR), and accuracy measures the rate of pixels correctly classified. These measurements are defined below as:

sensitivity = TPR = TP P = TP TP + FN (2)
specificity = TNR = TN N = TN TN + FP (3)
accuracy = TP + TN TP + TN + FP + FN , (4)
where TP stands for true positive, FP for false positive, TN for true negative and FN for false negative.

A benchmark of different approaches (including ours) is provided in Table 1. It shows the good performance of our proposed elongation-based upper leveling for both databases. In the case of the DRIVE database, our result is slightly under the best results given by the method of Mendonça [START_REF] Mendonça | Segmentation of retinal blood vessels by combining the detection of centerlines and morphological reconstruction[END_REF]. Note also that the approaches of Staal [START_REF] Staal | Ridge based vessel segmentation in color images of the retina[END_REF] and Niemeijer [START_REF] Niemeijer | Comparative study of retinal vessel segmentation methods on a new publicly available database[END_REF] are supervised approaches. On the STARE database, our proposed method performs also very well, and is even better than the method of Mendon¸ca [START_REF] Mendonça | Segmentation of retinal blood vessels by combining the detection of centerlines and morphological reconstruction[END_REF]. Both methods give results that are very close to one of the second human observer. Table 2 shows that our method is more robust than others, in the sense that it performs equivalently on both abnormal and normal images.

Last, note that the proposed elongation-based upper leveling is only a "simple" filtering step, whereas other approaches are more complicated. Besides, our process is not complete, since further post-processing can improve the results.

Extending constrained connectivity

From an algorithmic point of view, constrained connectivity [START_REF] Soille | Constrained connectivity for hierarchical image decomposition and simplification[END_REF] is the application of an increasing criterion (e.g., the range) on the Min-tree T of the minimum spanning tree (MST) [START_REF] Najman | On the equivalence between hierarchical segmentations and ultrametric watersheds[END_REF][START_REF] Najman | Playing with kruskal: algorithms for morphological trees in edge-weighted graphs[END_REF]. In the litterature, this tree T has been called the α-tree [START_REF] Ouzounis | Pattern spectra from partition pyramids and hierarchies[END_REF]. In this section, we apply our framework with a non-increasing criterion A derived from a popular work [START_REF] Felzenszwalb | Efficient graph-based image segmentation[END_REF]. Precisely, we use for T a binary Min-tree of the MST [START_REF] Najman | Playing with kruskal: algorithms for morphological trees in edge-weighted graphs[END_REF]. The dissimilarity used for the MST is the maximal distance of the red, blue, and green channels taken independently. The attribute A is derived from [START_REF] Felzenszwalb | Efficient graph-based image segmentation[END_REF], where the authors propose a region merging process that follows the edges of the MST by increasing order of the weights (dissimilarity). When an edge {x, y} is considered, they search for the regions X and Y that respectively contain the points x and y. The regions X and Y are merged if

Diff(X, Y ) < min{ Int(X) + k |X| , Int(Y ) + k |Y | }, (5)
where |.| denotes the cardinality, Diff(X, Y ) is the minimum weight of the edge connecting the two regions X and Y , Int(X) is the largest weight in the MST of the region X, and k is a parameter favouring the merging of small regions (a large k causes a preference for larger components). However, k is not a scale parameter in the sense of the causality principle: as shown in [START_REF] Guimaraes | An efficient hierarchical graph based image segmentation[END_REF] a contour present at a scale k 1 is not always present at a scale k 2 < k 1 . The merging criterion defined by Eq. ( 5) depends on the parameter k at which the regions X and Y are observed. So let us consider the attribute A as the k defined by k

= max Diff(X, Y ) -Int(X) × |X|, Diff(X, Y ) -Int(Y ) × |Y | . That is to say, for each node N k , let N c1 k and N c2 k be the two children of N k in T , then the attribute A for node N k is given by A(N k ) = max Diff(N c1 k , N c2 k ) -Int(N c1 k ) × |N c1 k |, Diff(N c1 k , N c2 k ) -Int(N c2 k ) × |N c2 k | . (6)
The maxima of the attribute A correspond to meaningful regions. We thus compute a Max-tree T T on the node-weigthed graph (T , A), with which we can Method GT Covering Prob. Rand. Index ODS OIS Best ODS OIS FH [START_REF] Felzenszwalb | Efficient graph-based image segmentation[END_REF] 0.43 0.53 0.68 0.76 0.79 Guimarães [START_REF] Guimaraes | An efficient hierarchical graph based image segmentation[END_REF] filter this graph. The contours of the flat zones of one level of filtering of (T T , A) provides a segmentation of the original image. Computing all levels of filtering produces a hierarchy of segmentations. A good representation of such a hierarchy is a saliency map [START_REF] Najman | Geodesic saliency of watershed contours and hierarchical segmentation[END_REF][START_REF] Najman | On the equivalence between hierarchical segmentations and ultrametric watersheds[END_REF], which is equal to the sum of all the contours of every filtered image. An efficient algorithm for computing the saliency map in our framework, based on the notion of extinction value [START_REF] Vachier | Extinction values: A new measurement of persistence[END_REF], will be provided in an extend version of this paper. The saliency map defines a duality between closed, non-self intersecting weighted contours and a hierarchy of regions. Low levels of the hierarchy correspond to weak contours, thus to over-segmentations. High levels of the hierarchy correspond to strong contours, hence to under-segmentations.

Moving between levels gives a continuous trade-off between those two extremes.

A given level can be seen as an observation scale at which we consider the image.

We have tested our extended constrained connectivity framework on the Berkeley Segmentation Dataset BSDS500 [START_REF] Arbelaez | Contour detection and hierarchical image segmentation[END_REF], an extension of the BSDS300 [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF]. The dataset consists of 500 natural images divided into 200 test images, 200 images for training, and 100 validation images, together with human annotations. Each image is segmented by an average of five different subjects. Fig. 7 and8 show some saliency maps computed on some images from the BSDS500 dataset. Two evaluation schemes are also provided by the authors. In the first one, the same fixed threshold level (observation scale) is used for all saliency maps in the dataset; we refer to it as the optimal dataset scale (ODS). In the second one, we evaluate the performance using an image-dependent threshold for each saliency map; we refer to this choice as the optimal image scale (OIS).

Quantitative evaluation is performed using the region-based performance measurements described in [START_REF] Arbelaez | Contour detection and hierarchical image segmentation[END_REF], in terms of Ground-Truth Covering criterion and Probabilistic Rand Index. Here, we compare our results with the graph-based image segmentation (Felz-Hutt) [START_REF] Felzenszwalb | Efficient graph-based image segmentation[END_REF], and with another method named hierarchical graph based image segmentation (Guimarães et al.) [START_REF] Guimaraes | An efficient hierarchical graph based image segmentation[END_REF], also relying on the same criterion popularized by [START_REF] Felzenszwalb | Efficient graph-based image segmentation[END_REF]. The comparison is given in Table 3. Our method ranks first, for both the optimal dataset scale (ODS) and for optimal image scale (OIS).

Discussion and conclusion

This paper has presented two applications of shape-based morphology, a generalization of existing tree-based connected operators. The first application uses one of the many shape-based upper levelings. Althought such a filter is but a "simple" filtering step, it gives results almost as good as the second human observer Fig. 7. Hierarchical segmentation results on the BSDS500. From left to right: Input image, saliency map, and segmentations at the optimal dataset scale (ODS) and at the optimal image scale (OIS).

in the case of blood vessels segmentation. The second application is an extension of the constrainted connectivity framework to non-increasing constraints. A quantitative evaluation based on a criterion given in [START_REF] Felzenszwalb | Efficient graph-based image segmentation[END_REF] shows that our approach compares favorably to previous works.

The potential of shape-based morphology is tremendous. In this paper, we have hardly skimmed the surface of what this theory has to offer to the scientific community. In [START_REF] Xu | Context-based energy estimator : Application to object segmentation on the tree of shapes[END_REF], we have also used this framework to achieve object segmentation on the tree of shapes. The key idea was to propose an efficient context-based energy estimator whose minima correspond to meaningful objects. Some other criterions that can be used are, for example an energy derived from the number of false alarms [START_REF] Cao | Extracting meaningful curves from images[END_REF] or some snake energy [START_REF] Kass | Snakes : Active contour models[END_REF]. A practical problem is that many minima of such energies do not correspond to meaningful components of the input image. In [START_REF] Xu | Context-based energy estimator : Application to object segmentation on the tree of shapes[END_REF], a morphological closing in the space of shapes helps to filter those spurious minima.

Implementations of shape-based filters are easy thanks to the open-source Milena library [START_REF] Levillain | Why and how to design a generic and efficient image processing framework: The case of the Milena library[END_REF]. A demonstration is available online from http://olena. lrde.epita.fr/ICPR2012 (see also http://olena.lrde.epita.fr/ICIP2012). More applications will be studied in a forthcoming paper. Properties of those filters, such as conditions for idempotence, will be also studied.

Acknowledgements. This work received funding from the Agence Nationale de la Recherche, contract ANR-2010-BLAN-0205-03 and through "Programme d'Investissements d'Avenir" (LabEx BEZOUT n • ANR-10-LABX-58). Fig. 8. Additional hierarchical segmentation results on the BSDS500. From top to bottom: Input mage, saliency map, and segmentations at the optimal dataset scale (ODS) and at the optimal image scale (OIS).

Fig. 1 .

 1 Fig. 1. Classical connected operators (black path) and our proposed shape-based morphology (black+red path).

Fig. 2 .

 2 Fig. 2. (a) Evolution of "circularity" on two branches of T ; (b to e): Some shapes; (f) Attribute thresholding; (g) Shaping.

  (a) Input image. (b) Low threshold. (c) Higher threshold. (d) Our shaping.

Fig. 3 .

 3 Fig. 3. Comparison of extinction-based shapings with attribute thresholding.

  (a) Input color image fc. (b) Green channel fg. (c) Inversed black top-hat fi. (d) Inversed filtering result.(e) Segmented blood vessels. (f) Manual segmentation.

Fig. 4 .

 4 Fig. 4. Illustration of the complete process of blood vessels segmentation in retinal image using elongation-based upper leveling (a-e); (f): Manual segmentation.

Fig. 5 .Fig. 6 . 2 nd

 562 Fig. 5. Illustration of four segmented blood vessels from the 20 test retinal images of the DRIVE database. Top: Input color retinal images; Bottom: Corresponding segmented results. (White pixels: true positive; Black pixels: true negative; Blue pixels: false positive; Red pixels: false negative.)

Table 1 .

 1 Benchmark of different blood vessels segmentation approaches on DRIVE and the STARE database.

				9473 (0.0048) 0.8949 0.9390 0.9354 (0.0171)
	mendonça [14] 0.7344 0.9764 0.9452 (0.0062) 0.6996 0.9730 0.9440 (0.0142)
	Staal [25]	0.7193 0.9773 0.9441 (0.0057)	-	-	-
	Niemeijer [19] 0.6793 0.9801 0.9416 (0.0065)	-	-	-
	Our	0.6924 0.9779 0.9413 (0.0078) 0.7149 0.9749 0.9471 (0.0114)
	Zana [31]	0.6696 0.9769 0.9377 (0.0078)	-	-	-
	Al-Diri [1] -	-	-	0.9258 (0.0126)	-	-	-
	Jiang [9]	0.6478 0.9625 0.9222 (0.0070)	-	-	0.9513
	Perez [12]	0.7086 0.9496 0.9181 (0.0240)	-	-	-
	Hoover [8]	-	-	-	0.6751 0.9567 0.9267 (0.0099)

Table 2 .

 2 Benchmark of vessels segmentation methods (STARE images -Normal versus abnormal cases).

	Method	Sensitivity Specificity	Accuracy
		Normal cases
	2 nd human observer	0.9646	0.9236	0.9283 (0.0100)
	Shape-based upper leveling 0.7178	0.9802 0.9493 (0.117)
	mendonça [14]	0.7258	0.9791	0.9492 (0.0122)
	Hoover [8]	0.6766	0.9662	0.9324 (0.0072)
		Abnormal cases
	2 nd human observer	0.0.8252	0.9544	0.9425
	Shape-based upper leveling 0.7120	0.9696 0.9447 (0.0106)
	mendonça [14]	0.6733	0.9669	0.9388 (0.0150)
	Hoover [8]	0.6736		0.9211 (0.0091)

Table 3 .

 3 Preliminary region benchmarks on the BSDS500.

		0.46 0.53 0.60 0.76	0.81
	Ours	0.50 0.57 0.66 0.77	0.82