Pierre Soille

Laurent Najman

On morphological hierarchical representations for image processing and spatial data clustering

Keywords: image representation, segmentation, clustering, ultrametric, hierarchy, graphs, connected components, constrained connectivity, watersheds, min-tree, alpha-tree

Hierarchical data representations in the context of classification and data clustering were put forward during the fifties. Recently, hierarchical image representations have gained renewed interest for segmentation purposes. In this paper, we briefly survey fundamental results on hierarchical clustering and then detail recent paradigms developed for the hierarchical representation of images in the framework of mathematical morphology: constrained connectivity and ultrametric watersheds. Constrained connectivity can be viewed as a way to constrain an initial hierarchy in such a way that a set of desired constraints are satisfied. The framework of ultrametric watersheds provides a generic scheme for computing any hierarchical connected clustering, in particular when such a hierarchy is constrained. The suitability of this framework for solving practical problems is illustrated with applications in remote sensing.

Introduction

Most image processing applications require the selection of an image representation suitable for further analysis. The suitability of a given representation can be evaluated by confronting its properties with those required by the application at hand. In practice, images are often represented by decomposing them into primitive or fundamental elements that can be more easily interpreted. Examples of decomposition (or simply representation) schemes are given hereafter:

-A functional decomposition decomposes the image into a sum of elementary functions. The most famous functional decomposition is the Fourier transform which decomposes the image into a sum of cosine functions with a given frequency, phase, and amplitude. This proves to be a very effective representation for applications requiring to target structures corresponding to well-defined frequencies; -A pyramid decomposition relies on a shrinking operation which applies a low-pass filter to the image and downsamples it by a factor of two and an expand operation which upsamples the image by a factor of two using a predefined interpolation method. Such a scheme is extremely efficient in situations where the analysis can be initiated at a coarse resolution and refined by going through levels of increasing resolution; -A multi-scale representation consists of a one-parameter family of filtered images, the parameter indicating the degree (scale) of filtering. This scheme is appropriate for the analysis of complex images containing structures at various scales; -A skeleton representation consists in representing the image by a thinned version. It is useful for applications where the geometric and topological properties of the image structures need to be measured; -The threshold decomposition decomposes a grey tone image into a stack of binary images corresponding to its successive threshold levels. This decomposition is useful as a basis for some hierarchical representations (see below) and from a theoretical point of view for generalising operations on binary images to grey tone images; -A hierarchical representation of an image can be viewed as an ordered set or tree (acyclic graph) with some elementary components defining its leaves and the full image domain defining its root. Examples of elementary components are the regional minima/maxima/extrema, or the flat zones of the input image. This approach is interesting in all applications where the tree encoding the hierarchy offers a suitable basis for revealing structural information for filtering or segmentation purposes.

Note that these schemes are not mutually exclusive. A case in point is the skeleton representation defined in terms of maximal inscribed disks since it fits the multi-scale representation (with morphological openings with disks of increasing size as structuring elements) as well as the functional decomposition (with spatially localised disks as elementary functions that are unioned to reconstruct the original pattern).

A given representation scheme can be further characterised by considering the properties of the operations it relies on. For example, a representation is linear if it is based on operations invariant to linear transformations of the input image. The multi-scale representation with Gaussian filters of increasing size fulfils this property. Morphological representations are non-linear representations relying on morphological operations. For example, a granulometry is a morphological multi-scale representation originally proposed by Matheron in his seminal study on the analysis of porous media [START_REF] Matheron | Eléments pour une théorie des milieux poreux[END_REF]. The representation does not need to rely exclusively on morphological operations to be considered as morphological. For example, the non-linear scale-space representation with levellings [START_REF] Meyer | Nonlinear scale-space representation with morphological levelings[END_REF] is based on self-dual geodesic reconstruction using Gaussian filters of increasing size as geodesic mask.

This paper deliberately focuses on hierarchical image representations for image segmentation with emphasis on morphological methods. Note that the development of hierarchical representations appeared first in taxonomy in the form of hierarchical clustering methods (see for example [START_REF] Cormack | A review of classification (with discussion)[END_REF] for an old but excellent review on classification including a discussion on hierarchical clustering). In fact, hierarchical image segmentation can be seen as a hierarchical clustering of spatial data. Graph theory is the correct setting for formalising clustering concepts as already recognised in [START_REF] Estabrook | A mathematical model in graph theory for biological applications[END_REF] and [START_REF] Matula | Cluster analysis via graph theoretic techniques[END_REF], see also the enlightening paper [START_REF] Zahn | Graph-theoretical methods for detecting and describing gestalt clusters[END_REF] as well as the detailed survey and connections between graph theory and clustering in [START_REF] Hubert | Some applications of graph theory to clustering[END_REF] (and [START_REF] Hubert | Min and max hierarchical clustering using asymetric similaritly measures[END_REF] for clustering on directed graphs). For this reason, Sec. 2 presents briefly background notions and notations of graph theory used throughout this paper. Then, fundamental concepts of hierarchical clustering methods where the spatial location of the data points is usually not taken into account are reviewed in Sec. 3. Hierarchical image segmentation methods where the spatial location of the observations (i.e., the pixels) plays a central role are presented in a nutshell in Sec. 4. Recent recent paradigms developed for the hierarchical representation of images in the framework of mathematical morphology known as constrained connectivity and ultrametric watersheds are then developed in Sec. 5 while highlighting their links with hierarchical clustering methods. The framework of ultrametric watersheds provides a generic scheme for computing any hierarchical connected clustering, in particular when such a hierarchy is constrained. Before concluding, the problem of transition pixels is set forth in Sec. 6.

Background definitions and notations on graphs

The objects under study (specimens in biology, galaxies in astronomy, or pixels in image processing) are considered as the nodes of a graph. An edge is then drawn between all pairs of objects that need to be compared. The comparison often relies on a dissimilarity measure that assigns a weight to each edge. Following the notations of [START_REF] Diestel | Graph Theory[END_REF], we summarise hereafter graph definitions required in the context of clustering.

A graph is defined as a pair X = (V, E) where V is a finite set and E is composed of unordered pairs of V , i.e., E is a subset of {{p, q} ⊆ V | p = q}. Each element of V is called a vertex or a point (of X), and each element of E is called an edge (of X). If V = ∅, we say that X is non-empty.

As several graphs are considered in this paper, whenever this is necessary, we denote by V (X) and by E(X) the vertex and edge set of a graph X.

Let X be a graph. If u = {p, q} is an edge of X, we say that p and q are adjacent (for X). Let π = p 0 , . . . , p be an ordered sequence of vertices of X, π is a path from p 0 to p in X (or in V ) if for any i ∈ [1, ], p i is adjacent to p i-1 . In this case, we say that p 0 and p are linked for X. We say that X is connected if any two vertices of X are linked for X.

Let X and Y be two graphs. If V (Y ) ⊆ V (X) and E(Y ) ⊆ E(X), we say that Y is a subgraph of X and we write Y ⊆ X. We say that Y is a connected component of X, or simply a component of X, if Y is a connected subgraph of X which is maximal for this property, i.e., for any connected graph Z, Y ⊆ Z ⊆ X implies Z = Y .

Clustering methods generally work on a complete graph (V, V × V ). In this case, the notion of connected component is not an important one, as any subset is obviously connected. On contrary, this notion is fundamental for image segmentation.

Let X be a graph, and let S ⊆ E(X). The graph induced by S is the graph whose edge set is S and whose vertex set is made of all points that belong to an edge in S, i.e., ({p ∈ V (X) | ∃u ∈ S, p ∈ u}, S).

In the sequel of this paper, X = (V, E) denotes a connected graph, and the letter V (resp. E) will always refer to the vertex set (resp. the edge set) of X. We will also assume that E = ∅. Let S ⊂ E. In the following, when no confusion may occur, the graph induced by S is also denoted by S. If S ⊂ E, we denote by S the complementary set of S in E, i.e., S = E \ S.

Typically, in applications to image segmentation, V is the set of picture elements (pixels) and E is any of the usual adjacency relations, e.g., the 4-or 8-adjacency in 2D [START_REF] Kong | Digital topology: Introduction and survey[END_REF]. In all examples, 4-adjacency is used.

We consider in this paper weighted graphs, and either the vertices or the edges of a graph can be weighted. We denote the weight on the vertives of V by f , and the weights on the edges of E by F . For application to image processing, f is generally some information on the pixels (e.g., the grey level of the considered pixel), and F represents a dissimilarity (e.g., F ({p, q}) = |f (p) -f (q)|).

Hierarchical clustering

Clustering can be defined as a method for grouping objects into homogeneous groups (called clusters) on the basis of empirical measures of similarity among those objects. Ideally, the method should generate clusters maximising their internal cohesion and external isolation. Analogously to the categorisation of classification methods proposed in [START_REF] Jones | Some thoughts on classification for retrieval[END_REF], any clustering methodology can be characterised by three main properties. The first concerns the relation between object properties and clusters. It indicates whether the clusters are monothetic or polythetic. A cluster is monothetic if and only if all its members share the same common property or properties. The second property regards the relation between objects and clusters. It indicates whether the clusters are exclusive (i.e., non-overlapping) or overlapping. Non-overlapping clustering methods can be defined as partitional in the sense that they realise a partition of the input objects (a partition of a set is defined as division of this set in disjoint non-empty subsets such that their union is equal to this set). Non-partitional clustering allows for overlap between clusters, see [START_REF] Jardine | A model for taxonomy[END_REF] for an early reference on this topic and [START_REF] Barth à C Lemy | Combinatorial optimization and hierarchical classifications[END_REF] for recent developments. The third property refers to the relation between clusters. It indicates whether the clustering method is hierarchical (also called ordered) or non-hierarchical (unordered).

Because we are chiefly interested in image segmentation applications, we focus on clustering methods that are monothetic, partitional, and hierarchical.

The term hierarchical clustering was first coined in [START_REF] Johnson | Hierarchical clustering schemes[END_REF]. A hierarchical clustering can be viewed as a sequence of nested clusterings such that a cluster at a given level is either identical to a cluster already existing at the previous level or is formed by unioning two or more clusters existing at the previous level. It is convenient to represent this hierarchy in the form of a tree called dendrogram [START_REF] Sokal | Principles of Numerical Taxonomy[END_REF] or taxonomic tree (see [START_REF] Sneath | The application of computers in taxonomy[END_REF] for this latter terminology as well as a procedure which in essence already defined the concept of hierarchical clustering). The first detailed study about the use of trees in the context of hierarchical clustering appeared in [START_REF] Hartigan | Representation of similarity matrices by trees[END_REF].

By construction, a hierarchical clustering is parameterised by a non-negative real number λ indicating the level of a given clustering in the hierarchy. At the bottom level, this number is equal to zero and each object correspond to a cluster so that the finest possible partition is obtained. At the top level only one cluster containing all objects remains. Given any two objects, it is possible to determine the minimum level value for which these two objects belong to the same cluster.

A key property of hierarchical clustering is that the function that measures this minimum level is an ultrametric. An ultrametric is a measurement that satisfies all properties of a metric (distance) plus a condition stronger than the triangle inequality and called ultrametric inequality. It states that the distance between two objects is lower than or equal to the maximum of the distances calculated from (i) the first object to an arbitrary third object and (ii) this third object to the second object. Denoting by d the ultrametric function and p, q, and r respectively the first, second and third objects, the ultrametric inequality corresponds to the following inequality: d(p, q) ≤ max{d(p, r), d(r, q)}.

The ultrametric property of hierarchical clustering was discovered simultaneously in [START_REF] Johnson | Hierarchical clustering schemes[END_REF][START_REF] Jardine | The structure and construction of taxonomic hierarchies[END_REF], see also [START_REF] Benz à C Cri | La taxinomie[END_REF] for a thorough study on ultrametrics in classification. An example of dendrogram is displayed in Fig. 1. The measure of similarity between the input objects requires the selection of a dissimilarity measurement. A dissimilarity measurement between the elements of a set V is a function d from V × V to the set of nonnegative real numbers satisfying the three following conditions: (i) d (p, q) ≥ 0 for all p, q ∈ V (i.e., positiveness), (ii) d (p, p) = 0 for all p ∈ V , and (iii) d (p, q) = d (q, p) for all p, q ∈ V (i.e., symmetry). Starting from an arbitrary dissimilarity measurement, it is possible to construct a hierarchical clustering: if the dissimilarity is increasing with the merging order, an ultrametric distance between any two objects (or clusters) can be defined as the dissimilarity threshold level from which these two objects (or clusters) belong to the same cluster; if if the dissimilarity is not increasing with the merging order, then any increasing function of the merging order can be used.

In practice, the hierarchy is constructed by an iterative procedure merging first the object pair(s) with the smallest dissimilarity value so as to form the first non-trivial cluster(s) (i.e., non reduced to one object). To proceed, the dissimilarity measurement between objects needs to be extended so as to be applicable to clusters. Let C i and C j denote two clusters obtained at a given iteration level. The dissimilarity between between these two clusters is naturally defined as a function of the dissimilarities between the objects belonging to these clusters:

d (C i , C j ) = f {d (p, q)) | p ∈ C i and q ∈ C j }.
Typical choices for the function f are the minimum or maximum. The maximum rule leads to the complete-linkage clustering (sometimes called maximum method) and dates back to [START_REF] Sørensen | A method of establishing groups of equal amplitude in plant sociology based on similarity of species content and its applications to analyses of the vegetation of Danish commons[END_REF]. Complete-linkage is subject to ties in case the current smallest dissimilarity value is shared by two or more clusters. Consequently, one of the possible merge must be chosen and often this can only be achieved by resorting to some arbitrary (order dependent or random) selection. By construction, complete-linkage favours compact clusters. On the other hand, the minimum rule is not subject to ties (and is therefore uniquely defined) and does not favour compact clusters. The resulting clustering is called the singlelinkage clustering1 (sometimes called minimum method). Indeed, only the pair (link) with the smallest dissimilarity value is playing a role.

The single-linkage clustering is closely related to the minimum spanning tree [START_REF] Gower | Minimum spanning trees and single linkage cluster analysis[END_REF], defined as follows. To any edge-weighted graph X, the number F (X) = u=∈E(X) F (u) is the weight of the graph. A spanning tree of a connected graph X is a graph whose vertex set is equal to V (X) and whose edge set is a subset of E(X) such that no cycles are formed. A spanning tree of X with minimum weight is called a minimum spanning tree of X.

Indeed, the hierarchy underlying the single-linkage clustering is at the root of the greedy algorithm of Kruskal [START_REF] Kruskal | On the shortest spanning subtree of a graph and the traveling salesman problem[END_REF] for solving the minimum spanning tree problem 2 . In this algorithm, referred to as 'construction A' in [START_REF] Kruskal | On the shortest spanning subtree of a graph and the traveling salesman problem[END_REF], the edges of the graph are initially sorted by increasing edge weights (in a clustering perspective, the nodes of the graph are the objects and the edge weights are defined by the dissimilarity measurements between the objects). Then, a minimum spanning tree M ST is defined recursively as follows: the next edge is added to M ST if and only if together with M ST it does not form a circuit. That is, there is a one-to-one correspondence between (i) the clusters obtained for a given dissimilarity level and (ii) the subtrees obtained for a distance equal to this level in Kruskal's greedy solution to the minimum spanning tree problem.

While the single-linkage is not subject to ties, it is sensitive to the presence of objects of intermediate characteristics (transitions) that may occur between two clearly defined populations, see [START_REF] Wishart | Mode analysis: a generalization of nearest neighour which reduced chain effect[END_REF] for a detailed discussion as well as Sec. 6. This effect is sometimes called 'chaining-effect' although this latter terminology is somewhat misleading for chaining is the very principle of single-linkage [START_REF] Jardine | The construction of hierarchic and non-hierarchic classifications[END_REF].

Hierarchical image segmentation

After a brief discussion on the definition of image segmentation and hierarchical image segmentation (see Sec. 4.1), methods relying on graph representations are presented (Sec. 4.2) and then those developed in MM (Sec. 4.2).

From image segmentation to hierarchical image segmentation

A segmentation of the definition domain V of an image is usually defined as a partition of V into disjoint connected subsets V i , . . . , V n (called segments) such that there exists a logical predicate P returning true on each segment but false on any union of adjacent segments [START_REF] Horowitz | Picture segmentation by a directed split-and-merge procedure[END_REF][START_REF] Zucker | Region growing: childhood and adolescence[END_REF]. That is, a series of subsets V i of the definition domain V of an image forms a segmentation of this image if and only if the following four conditions are met (i

) ∪ i (V i ) = X, (ii) V i ∩ V j = ∅ for all i = j, (iii) P (V i ) =
true for all i, and (iv) P (V i ∪ V j ) = false if V i and V j are adjacent. The first condition requires that every picture element (pixel) must belong to a segment. The second condition requires that each segment does not overlap any other segment. The third condition determines what kind of properties each segment must satisfy, i.e., what properties the image pixels must satisfy to be in the same segment. The fourth condition ensures that the segments are maximal in the sense that specifies that any merging of any adjacent regions would violate the third condition.

Note that uniqueness of the resulting segmentation given a predicate is not required. If uniqueness is desired, the predicate should rely on an equivalence relation owing to the one-to-one correspondence between the unique partitions of a set and the equivalence relations on it, see for example [31, p. 48]. Interestingly, the relation 'is connected' is an equivalence relation since it is reflexive (a point is connected to itself by a path of length 0), symmetric (if a point p is connected to a point q then q is connected to p since the reversal of a path is path), and transitive (if p is connected to q and q to r then p is connected to r since the concatenation of two paths is a path). Any given connectivity relation partitions the set of pixels of a given input image into equivalent classes called connected components [START_REF] Rosenfeld | Fuzzy digital topology[END_REF]. They are maximal subsets of pixels such that every pair of pixels belonging to such a subset is connected. The resulting partition meets therefore all conditions of a segmentation.

The segments resulting from a segmentation procedure are analogous to the clusters obtained when clustering data. Clustering techniques can be applied to image data for either classification or segmentation purposes. In the former case, the spatial position of the pixels does not necessarily play a role for clusters are searched in a parametric space such as the multivariate histogram. The resulting clusters partition the parametric space into a series of classes and this partition is used as a look-up-table to indicate the class of each pixel of the input image. An example of this approach using morphological clustering is proposed in [START_REF] Soille | Morphological partitioning of multispectral images[END_REF]. Contrary to data clustering applied to non-spatial data, the dissimilarity measurements between the data samples (i.e., the pixels) are not measured between all possible pairs. Indeed, the spatial position of the pixels plays a key role so that measurements are only performed between adjacent pairs of pixels. That is, the full dissimilarity matrix is very sparse: for a image of m×n pixels, there are 2mn -m -n entries in the (mn) 2 × (mn) 2 dissimilarity matrix when considering 4-adjacency relation.

By analogy with hierarchical clustering, hierarchical segmentation can be defined as a family of fine to coarse image partitions (i.e., family of ordered partitions) parameterised by a non-negative real number indicating the level of a given partition in the hierarchy. Hierarchical segmentation is useful to help the detection of objects in an image. In particular, it can be used to simplify the image in such a way that the elementary picture elements are not anymore the pixels but connected sets of pixels. Indeed, in image data, analogues to phonemes and characters correspond to structural primitives that compress the data to a manageable size without eliminating any possible final interpretations [START_REF] Ahuja | On detection and representation of multiscale low-level image structure[END_REF]. It should be emphasised that a hierarchical segmentation does not necessarily deliver segments directly corresponding to the searched objects. This happens for instance when an object is not characterised by some homogeneity/separation criteria but from the consideration of an a priori model of the whole object (e.g. perceptual grouping and Gestalt theory).

There exists a fundamental difference between segmentation and classification. Indeed, contrary to classification, segmentation requires the explicit definition of an adjacency graph or, more generally, a connection [START_REF] Serra | A lattice approach to image segmentation[END_REF][START_REF] Ronse | Partial partitions, partial connections and connective segmentation[END_REF]. Typically, the k-nearest neighbouring graph with k equal to 4 or 8 is used for processing 2-dimensional images. With classification, a decision about the class (i.e., label) of each pixel can be reached without using its spatial context (position) so that it does not necessarily need the definition of an adjacency graph. Nevertheless, any classification can be used to generate a segmentation. Indeed, once an adjacency graph is added to the classified image, the maximal connected regions of pixels belonging to the same class generate a segmentation of the image definition domain. If the considered adjacency graph is the complete graph, a one-to-one correspondence between the classes and the resulting connected components is obtained.

Hereafter, a selection of techniques achieving hierarchical image segmentation is proposed, extending the initial survey proposed in [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF]. We start with generic methods based on graph representations and then proceed with specific methods developed in the context of mathematical morphology. Recent developments related to constrained connectivity and ultrametric watersheds are discussed in Sec. 5.

Methods based on graph representations

Horowitz and Pavlidis [START_REF] Horowitz | Picture segmentation by a directed split-and-merge procedure[END_REF][START_REF] Horowitz | Picture segmentation by a tree traversal algorithm[END_REF] are among the first to suggest a formulation of hierarchical image segmentation in a graph theoretical framework. It is based on the split-and-merge algorithm. Because their implementation relies on a regular pyramid data structure with square blocks, it is not translation invariant and it favours blocky edges owing to the initial regular split of the image. In addition, the grouping stage of split-and-merge algorithms is order dependent, a drawback of all procedures updating the features of a region once new points are added to it.

The idea of applying the single-linkage clustering method to produce hierarchical image segmentation was implemented for the first time by Nagao [START_REF] Nagao | Region extraction and shape analysis in aerial photographs[END_REF][START_REF] Nagao | A Structural Analysis of Complex Aerial Photographs[END_REF] for processing aerial images using grey level differences between adjacent pixels as dissimilarity measurement. For colour images, the resulting dissimilarity vector led to the notion of differential threshold vector in [START_REF] Baraldi | Single linkage region growing algorithms based on the vector degree of match[END_REF]. The application of single-linkage clustering to image data are further developed in [START_REF] Morris | Graph theory for image analysis: an approach based on the shortest spanning tree[END_REF] using a graph theoretic framework. This latter paper also details a minimax SST (Shortest Spanning Tree) segmentation allowing for the initial minimum spanning tree to be partitioned into n subtrees by recursively splitting the subtree with the larger cost into 2 subtrees (see also recursive SST segmentation into n regions). Note that single-linkage clustering based on grey level difference dissimilarity was rediscovered much later in morphological image processing under the term quasi-flat zones [START_REF] Meyer | Morphological scale-space representation with levelings[END_REF][START_REF] Meyer | Nonlinear scale-space representation with morphological levelings[END_REF]. More recently, the more general and appropriate term of α-connected component was proposed in [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF] to refer to any connected component (i.e., maximal set of connected pixel) of pixels such that any pair of pixels of this connected component can be linked by a path such that the dissimilarity value between two successive pixels of the path does not exceed a given dissimilarity threshold value (see details in Sec. 5.1). The ultrametric behind the single-linkage hierarchical image segmentation is analogous to the one defined for single-linkage clustering, see Sec. 3.

The hierarchy of graphs (irregular pyramids) proposed recently in [START_REF] Nacken | Image segmentation by connectivity preserving relinking in hierarchical graph structures[END_REF][START_REF] Kropatsch | Grouping and segmentation in a hierarchy of graphs[END_REF] builds on the graph weighted partitions developed in [START_REF] Felzenszwalb | Image segmentation using local variations[END_REF][START_REF] Felzenszwalb | Efficient graph-based segmentation[END_REF] and inspired by the seminal work of Zahn [START_REF] Zahn | Graph-theoretical methods for detecting and describing gestalt clusters[END_REF] on point data clustering and its extension to graph cut image segmentation in [START_REF] Wu | An optimal graph-theoretic approach to data clustering: theory and its applications to image segmentation[END_REF][START_REF] Shi | Normalized cuts and image segmentation[END_REF]. It relies on weighted graphs where each element of the edge set is given a weight corresponding to the range of the values of its two nodes. The internal contrast of a connected component corresponds to the largest weight of all edges belonging to this connected component (an edge belongs to a connected component if its corresponding nodes belong to it or, alternatively, to a spanning tree of minimum sum of edge weights). The external contrast is defined as the smallest weight of the edges linking a pixel of the considered connected component to another one. The hierarchy is achieved by defining a dissimilarity measure accounting for both the internal and external contrasts. The successive levels of the hierarchy are then obtained by iteratively merging the adjacent connected components of minimum dissimilarity. An upto-date survey (including comparisons) of both regular and irregular pyramidal structures can be found in [START_REF] Marfil | Pyramid segmentation algorithms revisited[END_REF]. A survey on graph pyramids for hierarchical segmentation is proposed in [START_REF] Kropatsch | Multiresolution image segmentations in graph pyramids[END_REF].

The hierarchical image segmentation based on the notion of the cocoons of a graph relies on a complete-linkage hierarchy and its corresponding ultrametric [START_REF] Guigues | The hierarchy of the cocoons of a graph and its application to image segmentation[END_REF]. The same authors introduced the notion of scale-sets [START_REF] Guigues | Scale-sets image analysis[END_REF] where the dissimilarity measurement is replaced by a two-term energy minimization process where the first term accounts for the amount of information required to encode the deviation of the data against the region model (typically taken as the mean of the region) and the second term is proportional to the amount of information required to encode the shape of the model (typically taken as the boundary length of the region).

In [START_REF] Arbel áez | Energy partition and image segmentation[END_REF], the extrema mosaic (influence zones of the image regional extrema) is considered as the base level of the hierarchy. The dissimilarity between the segments is defined as the average gray level difference along the common boundary of these segments. This dissimilarity is increasing with the merging order and is therefore an ultrametric. Generic ultrametric distances obtained by integrating local contour cues along the regions boundaries and combining this information with region attributes are proposed in [START_REF] Arbel áez | Boundary extraction in natural images using ultrametric contour maps[END_REF].

Methods developed in mathematical morphology

Mathametical morphology relies on the notion of lattices, and a theory devoted to segmentation in this context recently appears [START_REF] Serra | A lattice approach to image segmentation[END_REF][START_REF] Ronse | Partial partitions, partial connections and connective segmentation[END_REF]. From a practical point of view, most of the application schemes use either a watershed-based approach or a tree-based approach.

Watershed based The waterfall algorithm [START_REF] Beucher | Segmentation d'images et morphologie mathématique[END_REF][START_REF] Beucher | Watershed, hierarchical segmentation and waterfall algorithm[END_REF][START_REF] Beucher | The morphological approach to segmentation: the watershed transformation[END_REF] can be considered as the first morphological hierarchical image segmentation method. The elementary components of the base level of the tree underlying the waterfall hierarchy are the catchment basins of the gradient of the image. Each basin is then set to the height of the lowest watershed pixel surrounding this basin while the watershed pixels keep their original value. The watersheds of the resulting image delivers basins corresponding to the subsequent level of the hierarchy. The procedure is then iterated until only one basin matching the image domain is obtained. This hierarchy of partitions can be implemented directly on graph data structures as detailed in [START_REF] Vincent | Watersheds in digital spaces: an efficient algorithm based on immersion simulations[END_REF].

Watershed hierarchies using the notion of contour dynamic is proposed in [START_REF] Najman | Geodesic saliency of watershed contours and hierarchical segmentation[END_REF]. The arcs of the watersheds of the gradient of the original image are valued by their contour dynamic. More precisely, the contour dynamic of an arc of a watershed separating two basins is defined as the height difference between the lowest point of this arc and the height of the highest regional minimum associated with these two basins. The contour dynamic is a dissimilarity that satisfies all properties of an ultrametric. The resulting contour dynamic map is a saliency map representing a hierarchy. Indeed, a fine to coarse family of partitions is obtained by thresholding the contour dynamic map for increasing contour dynamic values. By associating other dissimilarity measures to the arcs of the watersheds, other partition hierarchies are obtained.

Note that, if one wants to obtain theoretical results associating definitions and properties [START_REF] Cousty | Incremental algorithm for hierarchical minimum spanning forests and saliency of watershed cuts[END_REF], one has to work on edge-weighted graphs with the watershedcut definition [START_REF] Cousty | Watershed cuts: thinnings, shortest-path forests and topological watersheds[END_REF] that links the watershed with the minimum spanning tree as initially pointed out in [START_REF] Meyer | Minimum spanning forests for morphological segmentation[END_REF].

Tree based Another type of hierarchy is obtained by considering the flat zones of the image as the finest partition and then iteratively merging the most similar flat zones. This resulting tree is called binary partition trees in [START_REF] Salembier | Binary partition tree as an efficient representation for image processing, segmentation, and information retrieval[END_REF]. The tree always represents a hierarchy indexed by the merging order and not always the dissimilarity since the one used in [START_REF] Salembier | Binary partition tree as an efficient representation for image processing, segmentation, and information retrieval[END_REF] is not an ultrametric.

Another tree, known as the component tree [START_REF] Jones | Component trees for image filtering and segmentation[END_REF][START_REF] Jones | Connected filtering and segmentation using component trees[END_REF] of the vertices (called max-tree or min-tree in [START_REF] Salembier | Antiextensive connected operators for image and sequence processing[END_REF] depending on whether its leaves are matching the image maxima or minima) represents the hierarchy of the level sets of the image and are therefore not directly representing a hierarchy of partitions of the image definition domain.

However, when defined not on the vertices but on the edges, we will see below that the component tree is indeed a dendrogram representing a hierarchy of connected partitions.

Reviews on hierarchical methods developed in mathematical morphology based on watersheds are presented in [START_REF] Meyer | An overview of morphological segmentation[END_REF][START_REF] Meyer | Segmentation, minimum spanning tree and hierarchies[END_REF], and on trees in [START_REF] Salembier | Connected operators: A review of region-based morphological image processing techniques[END_REF][START_REF] Salembier | Connected operators based on tree pruning strategies[END_REF]. Recent developments related to constrained connectivity and ultrametric watersheds are developed in the next section.

Constrained connectivity and ultrametric watersheds

Constrained connectivity

Preliminaries Let us first recall the notion of α-connectivity that corresponds to single-linkage clustering applied to image data, see Sec. 4.2. Two pixels p and q of an image f are α-connected if there exists a path going from p to q such that the dissimilarity between any two successive pixels of this path does not exceed the value of the local parameter α. By definition, a pixel is α-connected to itself.

Accordingly, the α-connected component of a pixel p is defined as the set of image pixels that are α-connected to this pixel. We denote this connected component by α-CC(p): α-CC(p) = {p} ∪ q | there exists a path P = p = p 1 , . . . , p n = q , n > 1, such that F ({p i , p i+1 }) ≤ α for all 1 ≤ i < n . In the case of grey level images and when considering the absolute intensity difference as dissimilarity measure, the α-connected components of an image are equivalent to its quasiflat zones [START_REF] Meyer | Morphological scale-space representation with levelings[END_REF][START_REF] Meyer | Nonlinear scale-space representation with morphological levelings[END_REF]. Note that the edges of the connected graph corresponding to a given α-connected component is defined by the pairs of adjacent pixels belonging to this α-connected component such that their associated dissimilarity (weight) does not exceed α.

Definitions and properties

The constrained connectivity paradigm [START_REF] Soille | On genuine connectivity relations based on logical predicates[END_REF][START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF] originated from the need to develop a method preventing the formation of αconnected components whose range values exceed that specified by the local range parameter α (assuming that the dissimilarity between two pixels is the absolute difference of their intensity values, see [START_REF] Soille | Preventing chaining through transitions while favouring it within homogeneous regions[END_REF][START_REF] Gueguen | Frequent and dependent connectivities[END_REF] for other examples of dissimilarity measures). This is simply achieved by looking for the largest αconnected components satisfying a global range constraint referred to as the global range parameter denoted by ω:

(α, ω)-CC(p) = α i -CC(p) α i ≤ α and R α i -CC(p) ≤ ω ,
where the range function R calculates the difference between the maximum and the minimum values of a nonempty set of intensity values. Note that the (α, ω)connected components for α ≥ ω are equivalent to those obtained for α = ω. That is, when α ≥ ω the local range parameter does not play a role. This leads to the concept of (ω)-connected component3 :

(ω)-CC(p) = (α ≥ ω, ω)-CC(p) = α i -CC(p) | R α i -CC(p) ≤ ω .
The corresponding global dissimilarity measurement d Ω between two pixels is defined by the smallest range of the α-connected components containing these two pixels. This dissimilarity measurement satisfies also the ultrametric inequality. Accordingly, we obtain the following equivalent definition of a (ω)-connected component: (ω)-CC(p) = {q | d Ω (p, q) ≤ ω}. In contrast to what happens with the local dissimilarity measurement d A , the range of the values of arbitrary pairs of pixels belonging to the same (ω)-connected component is limited, the maximal value of this range being equal to ω. Therefore, the resulting clustering bears some resemblance to the complete linkage clustering suggested in [START_REF] Sørensen | A method of establishing groups of equal amplitude in plant sociology based on similarity of species content and its applications to analyses of the vegetation of Danish commons[END_REF] but, contrary to the latter procedure, it is unequivocal (see [16, pp. 181-182] for an account on the equivocality of the complete linkage clustering). The generalisation of the concept of constrained connectivity to arbitrary constraints is presented in [START_REF] Soille | On genuine connectivity relations based on logical predicates[END_REF].

Separation value

The separation value ∆ ∧ of an iso-intensity connected component (flat-zone) can be defined in terms of grey tone hit-or-miss transforms [START_REF] Soille | Advances in the analysis of topographic features on discrete images[END_REF] with adaptive composite structuring elements. The adaptive hit-or-miss transform of a pixel with the composite structuring element containing the origin o for the foreground component and its direct neighbours having a strictly lower value N < (o) for the background component outputs the difference between the input pixel value and that of its largest lower neighbour(s) if the set of its lower neighbours is non-empty, 0 otherwise. This adaptive hit-or-miss transform is denoted by HMT (o,N < (o)) :

[HMT (o,N < (o)) ](p) = f (p) -∨{f (q) | q ∈ N < (p)} if N < (p) = ∅, 0 otherwise.
Similarly, the adaptive hit-or-miss transform HMT (N > (o),o) of a pixel outputs the difference between the value of its smallest greater neighbour(s) and that of the pixel itself, if the set of its greater neighbours N > (o) is non-empty, 0 otherwise:

[HMT (N > (o),o) ](p) = ∧{f (q) | q ∈ N > (p)} -f (p) if N > (p) = ∅, 0 otherwise.
The non-zero values of the point-wise minimum between the two hit-or-miss transforms corresponds to the transition pixels in the sense that these pixels have simultaneously lower and greater neighbours (and the point-wise minimum image indicates the minimum height of the transition). The binary mask of transition pixels can therefore be obtained by the following operator denoted by TP:

TP = T >0 [HMT (o,N < (o)) ∧ HMT (N > (o),o) ].
In [START_REF] Soille | Constrained connectivity and transition regions[END_REF], the same mask is obtained by considering the non-zero values of the point-wise minimum of the gradients by erosion and dilation with the elementary neighbourhood (the pixel and its direct neighbours) as structuring element. In this latter case, the point-wise minimum image indicates the maximum height of the transition. The minimum separation value of a pixel of an image is defined as the minimum intensity difference between a pixel and its neighbour(s) having a different value from this pixel if such neighbour(s) exist, 0 otherwise. It is denoted by [∆ ∧ (f )](p) and can be calculated as follows:

[∆ ∧ (f )](p) =    [HMT (o,N < (o)) (f )](p) if [HMT (o,N < (o)) (f )](p) < [HMT (N > (o),o) (f )](p) and [HMT (o,N < (o)) (f )](p) = 0, [HMT (N > (o),o) (f )](p) otherwise.
The minimum separation value of an iso-intensity connected component 0-CC is then defined as the smallest (minimum) separation value of its pixels:

∆ ∧ (0-CC) = ∧{∆ ∧ (q) | q ∈ 0-CC and ∆ ∧ (q) = 0}.
It is equivalent to the smallest α value such that α-CC = 0-CC. Similarly, the operator that sets each pixel of the image to the minimum separation value of the iso-intensity connected component it belongs to is defined as follow:

[∆ ∧ (0-CC(f ))](p) = ∧{∆ ∧ (q) | q ∈ 0-CC(p) and ∆ ∧ (q) = 0}.
It can be viewed as an adaptive operation where the output value at a given pixel depends on the iso-intensity component of this pixel and the neighbouring pixels of this component. By replacing the ∧ operation with the ∨ operation in the minimum separation definitions, we obtain the definitions for maximal separations. Figure 2 illustrates the map of minimal separation of the pixels and iso-intensity connected components of a synthetic image.

1 2 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 2 5 1 5 2 1 2 1 1 1 1 1 1 3 2 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1
Fig. 2. Left: a synthetic 7 × 7 image f with its intensity values [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF]Fig. 2a]. Middle: the map of separation value of its pixels ∆ ∧ (f ). Right: the map of separation value of its flat zones ∆ ∧ (0-CC(f )).

The regional maxima RMAX of ∆ ∧ (0-CC(f )) can be used to flag the flat zones that are the most isolated. Conversely, the regional minima RMIN of ∆ ∧ (0-CC(f )) can be used to flag the flat zones from which an immersion simulation should be iniated to compute the successive levels of the hierarchy of constrained components. By doing so, an algorithm similar to the watershed by flooding simulation [START_REF] Soille | Determining watersheds in digital pictures via flooding simulations[END_REF] can be designed.

Alpha-tree representation Constrained connectivity relies on the definition of α-connectivity. The later boils down to the single-linkage clustering of the image pixels given the underlying dissimilarity measure between adjacent pixel pairs. The corresponding single-linkage dendrogram was described as a spatially rooted tree in [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF]. This spatially rooted tree was introduced as the alpha-tree in [START_REF] Ouzounis | Pattern spectra from partition pyramids and hierarchies[END_REF][START_REF] Ouzounis | Attribute-constrained connectivity and alpha-tree representation[END_REF]. It represents the fine to coarse hierarchy of partitions for an increasing value of the dissimilarity threshold α. The alpha-tree can also be seen as a component tree representing the ordering relations of the α-connected components of the image. The representation in terms of min-tree is developed in Sec. 5.2.

In the case of constrained connectivity, a given (α, ω)-partition corresponds to the highest cut of the alpha-tree such that all the nodes below this cut satisfy the α and ω constraints. Usually this cut is not horizontal. A given (ω)-partition corresponds to the highest cut of the alpha-tree such all the nodes below the cut satisfy the ω constraint. Alternatively, a (ω)-partition can be obtained by performing a horizontal cut in the dendrogram based on the ultrametric d Ω (i.e., the omega-tree). An example of omega-tree is given [START_REF] Soille | Constrained connectivity for the processing of very high resolution satellite images[END_REF]. Note however that the set of all (α, ω)-partitions is itself not ordered given the absence of order between arbitrary pairs of local and global dissimilarity threshold values.

Edge-weighted graph setting and minimum spanning tree By construction, the connected components of the graph G[α] = (V, {{p, q} ∈ E | F ({p, q}) ≤ α}) are equivalent to the α-connected components of f . Since α-connectivity corresponds to single-linkage clustering, there is an underlying minimum spanning tree associated to it (see also section 3 and [START_REF] Morris | Graph theory for image analysis: an approach based on the shortest spanning tree[END_REF] for equivalent image segmentations based on the direct computation of a minimum spanning tree). More precisely, the minimum spanning tree of the edge-weighted graph of an image is a tree spanning its pixels and such that the sum of the weights associated with the edges of the tree is minimal. Denoting by E min the edge set of a minimum spanning tree of the edge-weighted graph of an image, the connected components of the graph (V, {{p, q} ∈ E min | F ({p, q}) ≤ α}) are equivalent to those of G[α] (equivalent in the sense that given any node, the set of nodes of the connected component of (V, {{p, q} ∈ E min | F ({p, q}) ≤ α}) containing this node is identical to the set of nodes of the connected component of G[α] containing this very node). Since the minimum spanning tree representation contains less edges than the initial edge-weighted graph, it is less memory demanding for further computations such as global range computations. However, not all computations can be done on the minimum spannning tree (for example, connectivity constraints relying on the computation of a connectivity index [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF] cannot be derived from it).

Ultrametric watersheds: from hierarchical segmentations to saliency maps

We have several different ways to deal with hierarchies: dendrograms and minimum spanning trees. In the case where a hierarchy is made of connected regions, then we can also use its connected component tree, e.g., min-tree, max-tree or alpha-tree. None of these three tools allows for an easy visualisation of a given hierarchy as an image. We now introduce ultrametric watershed [START_REF] Najman | Ultrametric watersheds[END_REF][START_REF] Najman | On the equivalence between hierarchical segmentations and ultrametric watersheds[END_REF] as a tool that helps visualising a hierarchy: we stack the contours of the regions of the hierarchy; thus, the more a contour of a region is present in the hierarchy, the more visible it is. Ultrametric watershed is the formalisation and the caracterisation of a notion introduced under the name of saliency map [START_REF] Najman | Geodesic saliency of watershed contours and hierarchical segmentation[END_REF].

Ultrametric watersheds

The formal definition of ultrametric watershed relies on the topological watershed framework [START_REF] Bertrand | On topological watersheds[END_REF].

Let X be a graph. An edge u ∈ E(X) is said to be W-simple (for X) if X has the same number of connected components as X + u = (V (X), E(X) ∪ {u}).

An edge u such that F (u) = λ is said to be W-destructible (for F ) with lowest value λ 0 if there exists λ 0 such that, for all λ 1 , λ 0 < λ 1 ≤ λ, u is W-simple for G[λ 1 ] and if u is not W-simple for G[λ 0 ].

A topological watershed (on G) is a map that contains no W-destructible edges.

An ultrametric watershed is a topological watershed F such that F (v) = 0 for any v belonging to a minimum of F .

There exists a bijection between ultrametric distances and hierarchies of partitions [START_REF] Johnson | Hierarchical clustering schemes[END_REF]; in other word, to any hierarchy of partitions is associated an ultrametric, and conversely, any ultrametric yields a hierarchy of partitions, see also Sec. 3. Similarly, there exists a bijection between the set of hierarchies of connected partitions and the set of ultrametric watersheds [START_REF] Najman | Ultrametric watersheds[END_REF][START_REF] Najman | On the equivalence between hierarchical segmentations and ultrametric watersheds[END_REF]. In [START_REF] Cousty | Incremental algorithms for hierarchical minimum spanning forests and saliency of watershed cuts[END_REF], it is proposed a generic algorithm for computing hierarchies and their associated ultrametric watershed.

Usage: gradient and dissimilarity Constrained connectivity is a hierarchy of flat zones of f , in the sense where the 0-connected components of f are the zones of f where the intensity of f does not change. In a continuous world, such zones would be the ones where the gradient is null, i.e. ∇f = 0. However, the space we are working with is discrete, and a flat zone of f can consist in a single point. In general, it is not possible to compute a gradient on the points or on the edges such that this gradient is null on the flat zones. To compute a gradient on the edges such that the gradient is null on the flat zones, we need to "double" the graph, for example we can do that by doubling the number of points of V and adding one edge between each new point and the old one.

More precisely, if we denote the points of V by V = {p 0 , . . . , p n }, we set V = {p 0 , . . . , p n } (with V ∩ V = ∅), and

E = {{p i , p i } | 0 ≤ i ≤ n}. We then set V 1 = V ∪ V and E 1 = E ∪ E .
By construction, as G = (V, E) is a connected graph, the graph G 1 = (V 1 , E 1 ) is a connected graph. We also extend f to V , by setting, for any p ∈ V , f (p ) = f (p), where {p, p } ∈ E .

We set, as in section 5.1, F ({p, q}) = |f (p) -f (q)|. The map F can be seen as the "natural gradient" of f [START_REF] Mattiussi | The Finite Volume, Finite Difference, and Finite Elements Methods as Numerical Methods for Physical Field Problems[END_REF]. We can then apply the same scheme on this F as in section 5.1 to find the hierarchy of α-connected components.

We denote by L(G 1 ) the edge graph (also called line graph) of G 1 . That is, each vertex of L(G 1 ) represents an edge of G 1 and two vertices of L(G 1 ) are adjacent if and only if their corresponding edges in G 1 share a common endpoint in G 1 . While the edges of L(G 1 ) are not weighted, the weights of its nodes are given by the weights of the corresponding edges of G 1 . It follows that the minima of L(G 1 ) are equivalent to the 0-connected components of G 1 . More generally, the alpha-tree of G 1 is contained in the min-tree of L(G 1 ). Interestingly, the mintree of L(G 1 ) can be computed efficiently thanks to the quasi-linear algorithm described in [START_REF] Najman | Building the component tree in quasi-linear time[END_REF]. Hence, the morphological framework of attribute filtering [START_REF] Breen | Attribute openings, thinnings, and granulometries[END_REF] can be applied to this min-tree [START_REF] Jones | Component trees for image filtering and segmentation[END_REF][START_REF] Salembier | Antiextensive connected operators for image and sequence processing[END_REF][START_REF] Jones | Connected filtering and segmentation using component trees[END_REF], similarly to the segmentation of an image into k regions proposed in [START_REF] Cousty | Watershed cuts: minimum spanning forests and the drop of water principle[END_REF]. This is in particular useful when the filtering is performed before computing a watershed and this is illustrated in the next paragraph for the computation of a hierarchy based on constrained connectivity.

Finding the (α, ω)-CCs can be done by filtering the ultrametric watershed W of F with R that acts as a flooding on the topological/ultrametric watershed W of F , and then finding a (topological) watershed of the filtered image. Repeating these steps for a sequence of ordered (α, ω) vectors, we build a constrained connectivity hierarchy. In effect, we are viewing a hierarchy as an image (edgeweighted graph) and transforming it into another hierarchy/image. Thus, classical tools from mathematical morphology can be applied to constrain any hierarchy. Similar examples exist in the literature, for example [START_REF] Guigues | Scale-sets image analysis[END_REF], where the authors compute what they called a non-horizontal cut in the hierarchy, in other words, they compute a flooding on a watershed. In their framework, the flooding is controlled by an energy.

The advantages of using an ultrametric watershed are numerous. Let us mention the two following ones:

1. an ultrametric watershed is visible. A dendrogram or a component tree can be drawn, but less information is available from such a drawing, and visualising a MST is not really useful; 2. an ultrametric watershed allows the use any information in the contours between regions; such information is not available on the component tree, and is only partially available with a MST (which contains only the pass between regions).

Let us note that those concepts are theoretically equivalent: even their respective computational time is in practice nearly identical; thus we can choose the one the most adapted to the desired usage.

Visualising the hierarchy of constrained connectivity as an ultrametric watershed allows ones to assess some of its qualities. One can notice in Fig. 3.c a large number of transition regions (small undesirable regions that persist in the hierarchy), which is the topic of the next section.

Transition pixels

Constrained connectivity prevents the formation of connected components that would otherwise be created in case samples of intermediate value (transition pixels) between two populations (homogeneous image structures) are present. Indeed, these components would violate the global range or other appropriate constraint. However, sometimes the formation of two distinct connected components cannot occur at all. In the extreme case represented in Fig. 4. either each pixel is a connected component (flat zone) or there is a unique connected component. One way to address this problem is to propose a definition of transition pixels and perform some pre-processing to suppress them. This approach is advocated in [START_REF] Soille | Constrained connectivity and transition regions[END_REF][START_REF] Soille | Constrained connectivity for the processing of very high resolution satellite images[END_REF]. For example, assuming that local extrema correspond to non-transition pixels, they are extracted on then considered as seeds whose values are propagated in the input image using a seeded region growing algorithm [START_REF] Adams | Seeded region growing[END_REF].

Note that this approach is linked with contrast enhancement techniques since it Fig. 4. A synthetic sample image with its intensity values and its two possible partitions into constrained connected components whatever the considered constraints in case standard α-connectivity is used in the definitions. The two homogeneous regions show intensity variations of 1 level while the ramp between the two regions also proceeds by steps of 1 intensity level. In the image at the right, adjacent pixels are linked by an edge if and only if their range does not exceed 1.

aims at increasing the external isolation of the obtained connected components.

A number of classical morphological schemes (e.g., area filtering of the ultrametric watershed) can be used to remove those transition zones (see Fig. Another approach is to substitute the α-connectivity with a more restrictive connectivity. Indeed, the local range parameter α defined in [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF] as the intensity difference between adjacent pixels can be viewed as a special case of dissimilarity measurement. Although this measurement is the most natural, other dissimilarity measurements may be considered. For example, the following alternative definition of alpha-connectivity may be considered to tackle the problem of transition regions. Let the α-degree of a pixel (node) be defined as the number of its adjacent pixels that are within a range equal to α: α -deg(p) = #{q | {p, q} ∈ E and |f (q) -f (p)| ≤ α}.

Then two pixels p and q are said to be α n -connected if and only if there exists an α-path connecting them such that every pixel of the path has a α-degree greater of equal to n. We obtain therefore the following definition for the α n -connected component of a pixel p: If necessary, other constraints can be considered. Note that α-connectivity is a special case of α n -connectivity obtained for n = 1. In addition, the following nesting property holds: the non singleton 1 3 -connected components match the core of the two homogeneous regions. Singleton connected components correspond to pixels whose degree is smaller than 3. Non-singleton connected components can be used as seeds for coarsening the obtained partition. Special care is needed to produce connected components matching one-pixel thick non-transition regions. Alternative approaches to tackle the problem of transition regions are also presented in [START_REF] Soille | Preventing chaining through transitions while favouring it within homogeneous regions[END_REF] using a dissimilarity value taking into account the values of the gradient by erosion and dilation at the considered adjacent pixels and in [START_REF] Gueguen | Frequent and dependent connectivities[END_REF] using image statistics.

Conclusion and perspectives

In this paper, we have presented several equivalent tools dealing with hierarchies of connected partitions. Such a review invites us to look more closely at links between what have been done in different research domains as, for example, between clustering and lattice theory [START_REF] Hubert | Some extension of Johnson's hierarchical clustering[END_REF]. A first step in that direction is [START_REF] Cousty | Raising in watershed lattices[END_REF], and there is a need for in-depth study of operators acting on lattices of graphs [START_REF] Cousty | Some morphological operators in graph spaces[END_REF] (or the one of complexes [START_REF] Dias | Some Morphological Operators on Simplicial Complex Spaces[END_REF]). The question of transition pixels is not only a theoretical one, regarding its significance for applications. Finally, we want to stress the importance of having frame work allowing a generic implementation of existing algorithms, not limited to the pixel framework, but also able to deal transparently with edges, or, more generally, with graphs and complexes [START_REF] Levillain | Writing reusable digital topology algorithms in a generic image processing framework[END_REF]. Finally, when dealing with very large images such as those encountered in remote sensing or biomedical imaging, the computation of the min-tree of the edge graph of an image may be prohibitive in terms of memory needs (without mentioning the additional cost of doubling the graph to make sure that each flat zone of the original image is matched by a minimum of the edge graph). In this situation, the direct computation of the alpha-tree of the image may be a valid alternative. An efficient implementation based on the union-find as originally presented for the computation of component trees [START_REF] Najman | Building the component tree in quasi-linear time[END_REF] is presented in [START_REF] Ouzounis | Attribute-constrained connectivity and alpha-tree representation[END_REF].

Fig. 1 .

 1 Fig. 1. An example of dendrogram starting from 6 objects at the bottom of the hierarchy (level λ = 0). At the top of the hierarchy, there remains only one cluster containing all objects.

Fig. 3 .

 3 Fig. 3. Constrained connectivity and ultrametric watersheds. (a) Original image (extract from the panchromatic channel of a Quickbird Imagery c DigitalGlobe Inc., 2007, distributed by Eurimage). (b) Ultrametric watershed W 1 for the α-connectivity (the grey level of a contour corresponds to the α value above which the contour disappears in the α-hierarchy. (c) Ultrametric watershed W 2 for the constrained connectivity (the grey level of a contour corresponds to the α = ω value above which the contour disappears in the (α, ω = α)-hierarchy). (d) Ultrametric watersheds corresponding to one of the possible hierarchies of area-filterings on W 2 .

  3.d for an example).

  α n -CC(p) = {p} ∪ {q | there exists a path p = p 1 , . . . , p n = q , n > 1, such that |f (p i ) -f (p i+1 | ≤ α and α -deg(p i ) ≥ n}.

Fig. 5 .

 5 Fig.5. A synthetic sample image with its intensity values, the corresponding 1-deg map, and 13-connected components.

The concept of single-linkage and its use for classification purposes were apparently suggested for the first time in[START_REF] Florek | Sur la liaison et la division des points d'un ensemble fini[END_REF] while the terminology single-linkage seems to be due to Sneath, see[16, p. 180] where it is also called Sneath's method.

The first explicit formulation of the minimum spanning tree problem is attributed to[START_REF] Boruvka | O jistém problému minimálním [On a certain minimal problem[END_REF], see detailed account on the history of the problem in[START_REF] Graham | On the history of the minimum spanning tree problem[END_REF].

The parenthesis is not dropped to avoid confusion with α-connected components when the Greek letters are replaced by a numerical value indicating the actual value of the corresponding range parameter.


Roland Levillain
email: roland.levillain@lrde.epita.fr

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

Laurent Najman
email: l.najman@esiee.fr

Writing Reusable Digital Topology Algorithms in a Generic Image Processing Framework

Digital Topology software should reflect the generality of the underlying mathematics: mapping the latter to the former requires genericity. By designing generic solutions, one can effectively reuse digital topology data structures and algorithms. We propose an image processing framework focused on the Generic Programming paradigm in which an algorithm on the paper can be turned into a single code, written once and usable with various input types. This approach enables users to design and implement new methods at a lower cost, try cross-domain experiments and help generalize results.

Introduction

Like Mathematical Morphology (MM), Digital Topology (DT) has many applications in image analysis and processing. Both present sound mathematical foundations to handle many types of discrete images. In fact most methods from Mathematical Morphology or Digital Topology are not tied to a specific context (image type, neighborhood, topology): they are most often described in abstract and general terms. Thus they are not limiting their field of application. However, software packages for MM and DT rarely take (enough) advantage of this generality: an algorithm is sometimes reimplemented for each image and/or each value type, or worse, written for a unique input type. Such implementations are not reusable because of their lack of genericity. These limitations often come from the implementation framework, which prohibits a generic design of algorithms. A recent and notable exception is the DGtal project, which proposes Digital Geometry (DG) software tools and algorithms built in a generic C ++ framework [START_REF]DGtal: Digital geometry tools and algorithms[END_REF].

Thanks to the Generic Programming (GP) paradigm, provided in particular by the C ++ language, one can design and implement generic frameworks. This paradigm is especially well-suited to the field of scientific applications where the efficiency, widespread availability and standardization of C ++ are real assets. To this end, we have designed a paradigm dedicated to generic and efficient scientific software [START_REF] Géraud | Semantics-driven genericity: A sequel to the static C++ object-oriented programming paradigm (SCOOP 2)[END_REF] and applied the idea of generic algorithms to MM in Image Processing (IP) [START_REF] Levillain | Milena: Write generic morphological algorithms once, run on many kinds of images[END_REF], as suggested by d'Ornellas and van den Boomgaard [START_REF] Ornellas | The state of art and future development of morphological software towards generic algorithms[END_REF]. The result of our experiments is a generic library, Milena, part of the Olena image processing platform [START_REF]The Olena image processing platform[END_REF].

Lamy suggests to implement digital topology in IP libraries [START_REF] Lamy | Integrating digital topology in image-processing libraries[END_REF]. The proposed solution, applied to the ITK library [START_REF] Ibáñez | the Insight Software Consortium: The ITK Software Guide[END_REF][START_REF]National Library of Medicine: Insight segmentation and registration toolkit (ITK)[END_REF] "works for any image dimension". In this paper, we present a framework for the generic implementation of DT methods within the Milena library, working for any image type supporting the required notions (value types, geometric and topological properties, etc.). Such a generic framework requires the definition of concepts from the domain (in particular, of an image) to organize data structures and algorithms, as explained in Sect. 2. Given these concepts it is possible to write generic algorithms, like a homotopic thinning operator making use of various definitions of the notion of simple point. We present a generic definition of such an operator in Sect. 3 and show some illustrations in Sect. 4. Section 5 concludes on the extensibility of this work along different axes: existing algorithms, new data structures and efficiency.

Genericity in Image Processing

In order to design a generic framework for image processing, we have previously proposed the following definition of an image [START_REF] Levillain | Milena: Write generic morphological algorithms once, run on many kinds of images[END_REF].

Definition. An image I is a function from a domain D to a set of values V ; the elements of D are called the sites of I, while the elements of V are its values.

For the sake of generality, we use the term site instead of point; e.g. a site could represent a triangle of a surface mesh used as the domain of an image. Classical site sets used as image domains encompass hyperrectangles (boxes) on regular n-dimensional grids, graphs and complexes (see Sect. 3).

In the GP paradigm, these essential notions (image, site set, site, value) must be translated into interfaces called concepts in Milena (Image, Site Set, etc.) [START_REF] Levillain | Why and how to design a generic and efficient image processing framework: The case of the Milena library[END_REF]. These interfaces contain the list of services provided by each type belonging to the concept, as well as its associated types. For instance, a type satisfying the Image concept must provide a domain() routine (to retrieve D), as well as a domain t type (i.e. the type of D) satisfying the Site Set concept. Concepts act as contracts between providers (types satisfying the concept) and users (algorithms expressing requirements on their inputs and outputs through concepts). For instance, the breadth first thinning routine from Algorithm 1.3 expects the type I (of the input image) to fulfill the requirements of the Image concept. Likewise nbh must be a Neighborhood; and is simple and constraint must be functions taking a value of arbitrary type and returning a Boolean value (Function v2b concept).

Generic Implementation of Digital Topology

Let us consider the example of homotopic skeletonization by thinning. Such an operation can be obtained by the removal of simple points (or simple sites in the Milena parlance) using Algorithm 1.1 [START_REF] Bertrand | Transformations topologiques discrètes[END_REF]. A point of an object is said to be simple if its deletion does not change the topology of the object. This algorithm takes an object X and a constraint K (a set of points that must not be removed) and iteratively deletes simple points of X\K until stability is reached. Algorithm 1.1 is an example of an algorithm with a general definition that could be applied to many input types in theory. But in practice, software tools often allow a limited set of such input types (sometimes just a single one), because some operations (like "is simple") are tied to the definition of the algorithm [START_REF] Levillain | Milena: Write generic morphological algorithms once, run on many kinds of images[END_REF].

Algorithm 1.2 shows a more general version of Algorithm 1.1, where implementation-specific elements have been replaced by mutable parts: a predicate stating whether a point p is simple with respect to a set X (is simple); a routine "detaching" a (simple) point p from a set X (detach); and a predicate declaring whether a condition (or a set of conditions) on p is satisfied before considering it for removal (constraint). The algorithm takes these three functions as arguments in addition to the input X. Algorithm 1.2 is a good candidate for a generic C ++ implementation of the breadth-first thinning strategy and has been implemented as Algorithm 1.3 in Milena1 . This algorithms implements the breadth-first traversal by using a FIFO (First In, First Out) queue. The set X is represented by a binary image (V = {true, false}), that must be compatible with operations performed within the algorithm. Inputs is simple, detach and constraint2 have been turned into function objects (also called functors). The breadth first thinning routine creates and returns an image with type mln concrete(I); it is an image type equivalent to I that allows to store data for every sites independently (which is not the case for some image types).

Simple Point Characterization Implementation

There are local characterizations of simple points in 2D, 3D and 4D, which can lead to look-up table (LUT) based implementations [START_REF] Couprie | New characterizations of simple points in 2D, 3D, and 4D discrete spaces[END_REF]. However, since the number of configurations of simple and non-simple points in Z d is 2 3 d -1 , this approach can only be used in practice in 2D (256 configurations, requiring a LUT of 32 bytes) and possibly in 3D (67,108,864 configurations, requiring a LUT of 8 megabytes). The 4D case exhibits 2 80 configurations, which is intractable using a LUT, as it would need 128 zettabytes (128 billions of terabytes) of memory. Couprie and Bertrand have proposed a more general framework for checking for simple points using cell complexes [START_REF] Couprie | New characterizations of simple points in 2D, 3D, and 4D discrete spaces[END_REF] and the collapse operation. Intuitively, complexes can be seen as a generalization of graphs. An informal definition of a simplicial complex (or simplicial d-complex) is "a set of simplices" (plural of Algorithm 1.1. Breadth-First Thinning.

Data : E (a set of points/sites), X ⊆ E (initial set of points), K ⊆ X (a set of points (constraint) that cannot be removed)

Result : X P ← { p ∈ X | p is simple for X } while P = ∅ do S ← ∅ for each p ∈ P do if p ∈ K and p is simple for X then X ← X -{p} for each n ∈ N (p) ∩ X do S ← S ∪ {n} P ← ∅ for each p ∈ S do if p is simple for X then P ← P ∪ {p} Algorithm 1.2.
A generic version of Algorithm 1.1.

Data : E, X ⊆ E, N (neighborhood), is simple (a function saying whether a point is simple), detach (a routine detaching a point from X), constraint (a function representing a constraint) simplex), where a simplex or n-simplex is the simplest manifold that can be created using n points (with 0 ≤ n ≤ d). A 0-simplex is a point, a 1-simplex a line segment, a 2-simplex a triangle, a 3-simplex a tetrahedron. A graph is indeed a 1-complex. Figure 1(a) shows an example of a simplicial complex. Likewise, a cubical complex or cubical d-complex can be thought as a set of n-faces (with 0 ≤ n ≤ d) in Z d , like points (0-faces), edges (1-faces), squares (2-faces), cubes (3-faces) or hypercubes (4-faces). Figure 1(b) depicts a cubical complex sample.

Result : X P ← { p ∈ X | is simple(p, X) } while P = ∅ do S ← ∅ for each p ∈ P do if constraint (p) and is simple(p, X) then X ← detach(X, p) for each n ∈ N (p) ∩ X do S ← S ∪ {n} P ← ∅ for each p ∈ S do if is simple (p,
Complexes support a topology-preserving transformation called collapse. An elementary collapse removes a free pair of faces of a complex, like the square face f 1 and its top edge e 1 , or the edge e 2 and its top vertex v, in Fig. 1(b). The pair (f 2 , e 3 ) cannot be removed, since e 3 also belongs to f 3 . Successive elementary collapses form a collapse sequence that can be used to remove simple points. Collapse-based implementations of simple-point deletion can always be used in 2D, 3D and 4D, though they are less efficient than their LUT-based counterparts. On the other hand, they provide some genericity as the collapse operation can have a single generic implementation on complexes regardless of their structure.

Illustrations

Using this generic approach, Algorithm 1.3 can be used to compute skeletons of various input images.

Skeleton of a 2D Binary Image

Our first illustration uses a classical 2D binary image built on a square grid (Fig. 2 I and N are introduced as aliases of the image and neighborhood types for convenience. The breadth first thinning algorithm is called with five arguments, as expected. The first two ones are the input image the (4-connectivity) neighborhood used in the algorithm. The last three ones are the functors governing the behavior of the thinning operator. The call is simple point2d<I, N>(c4(), c8()) creates a simple point predicate based on the computation of the 2D connectivity numbers [START_REF] Bertrand | Transformations topologiques discrètes[END_REF] associated with the 4-connectivity for the foreground and the 8-connectivity for the background. To compute these numbers efficiently, is simple point2d uses a LUT containing all the possible configurations in the 8-connectivity neighborhood of a pixel. detach point<I> is a simple functor removing a pixel by giving it the value "false". Finally, no constraint is an empty functor representing a lack of constraint. We also present a variation of the previous example where the fifth argument passed to the function is an actual constraint, preserving all end points of the initial image (see Fig. 2(c)). This result is obtained by invoking the generic functor is not end point in the following lines. This call creates a predicate characterizing end points by counting their number of neighbors.

I output_with_end_points =

breadth_first_thinning(input, c4(), is_simple_point2d<I, N>(c4(), c8()), detach_point<I>(), is_not_end_point<I, N>(c4(), input));

Skeleton of a 3D Binary Image

This second example in 3D is similar to the previous one in 2D. The domain of the image is a box on a cubical grid; the 26-and the 6-connectivity are respectively used for the foreground and the background. The output of Fig. 3(b) is obtained from the 3D volume shown in Fig. 3 The image of curvature 4(b) is computed on the edges of the mesh, and simplified using an area opening filter. All curvature regional minima are then removed from the mesh, and the skeleton 4(c) is obtained with Algorithm 1.3 using the collapse operation.

The only real difference with the previous example is the use of the functor is simple point3d. The default implementation of this predicate uses an onthe-fly computation of 3D connectivity numbers. We have also implemented a version based on a precomputed LUT which showed significant speed-up improvements.

Please note that the predicates is simple point2d and is simple point3d are specifically defined for a given topology in order to preserve performances.

Thick Skeleton of a 3D Mesh Surface

In this third example, we manipulate discrete mesh surfaces composed of triangles. The input of the thinning operator is a surface containing "holes", obtained from the mesh shown in Fig. 4(a) by removing triangles located in regional minima of the surface's curvature (darkest areas of Fig. 4(b)). The result presented in Fig. 4(c) is obtained with the following lines. Types are not shown to make this code more readable.

output = breadth_first_thinning(input, nbh , is_simple_triangle, detach_triangle, no_constraint());
In the previous code, input is a triangle-mesh surface represented by an image built on a simplicial 2-complex and nbh represents an adjacency relationship between triangles sharing a common edge. Function objects is simple triangle and detach triangle are operations compatible with input's type; they are generic routines based on the collapse operation mentioned in Sect. 3, working with any complex-based binary image. The input image is constructed so that the sites browsed by the for all loops in Algorithm 1.3 are only 2-faces (triangles), while preserving access to values at 1-faces and 0-faces. Thus, even though they receive 2-faces as input parameters, is simple triangle and detach triangle are able to inspect the adjacent 1-faces and 0-faces and determine whether and how a triangle can be completely detached from the surface through a collapse sequence.

The resulting skeleton is said to be thick, since it is composed of triangles connected by a common edge. The corresponding complex is said to be pure, as it does not contain isolated 1-faces or 0-faces (that are not part of a 2-face).

Thin Skeleton of a 3D Mesh Surface

To obtain a thin skeleton, we can use a strategy based on successive n-collapse operations, with n decreasing [START_REF] Cousty | Collapses and watersheds in pseudomanifolds[END_REF]. From the input of the previous example, we can obtain a ultimate 2-collapse by removing all simple pairs composed of a 2face and a 1-face (a triangle and an adjacent edge). The following lines compute such an ultimate 2-collapse. The iteration on input's domain is still limited to triangles (2-faces). collapse2 = breadth_first_thinning(input, nbh, is_triangle_in_simple_pair, detach_triangle_in_simple_pair, no_constraint());

Functor is triangle in simple pair checks whether a given triangle is part of a simple pair, and if so detach triangle in simple pair is used to remove the pair. Thinning the initial surface with this "simple site" definition produces a mesh free of 2-faces (triangles), as shown in Fig. 5(a). From this first skeleton, we can compute an ultimate 1-collapse, by removing all simple pairs composed of an edge (1-face) and a vertex (0-face). This skeleton is produced with the following code, where input2 is an image created from collapse2, and for which the domain of has been set to the edges of the complex, (instead of the triangles). collapse1 = breadth_first_thinning(input2, nbh, is_edge_in_simple_pair, detach_edge_in_simple_pair, no_constraint());

Here is edge in simple pair and detach edge in simple pair respectively test and remove an edge along with a vertex that form a simple pair. The result is a simplified skeleton, with no isolated branches, as the lack of constraint (no constraint) does not preserve them. The output of the ultimate 1-collapse on the bunny mesh is depicted in Fig. 5(b). It contains the crest lines that form the boundaries of catchment basins, such as in the watershed transform, and, in addition, the crest lines that make the previous ones connect one to another. Note that in both cases, the neighborhood object nbh is the same, as it represents the adjacency of two n-faces connected by a common adjacent (n-1)face. In the case of the 2-collapse, the neighborhood of a site (triangle) is the set of adjacent triangles connected by an edge, while in the case of the 1-collapse, the neighborhood of a site (edge) is the set of adjacent edges connected by a vertex.

Execution Times

Table 1 shows the execution times of the previous illustrations, computed on a PC running Debian GNU/Linux 6.0.4, featuring an Intel Pentium 4 CPU running at 3.4 GHz with 2 GB RAM at 400 MHz, using the C ++ compiler g++ (GCC) version 4.4.5, invoked with optimization option '-03'. The first three test cases use a simple point criterion based on connectivity numbers, while the last three use a collapse-based definition.

Conclusion

We have presented building blocks to implement reusable Digital Topology algorithms in an Image Processing framework, Milena. Given a set of theoretical constraints on its inputs, an algorithm can be written once and reused with many compatible image types. This design has previously been proposed for Mathematical Morphology, and can be applied to virtually any image processing field. Milena is Free Software released under the GNU General Public License, and can be freely downloaded from http://olena.lrde.epita.fr/.

Input

Input size Constraint Output Time 2D image (Fig. 2 .18 s (Fig. 5(a))

(1-collapse)

A strength of generic designs is their ability to extend and scale easily and efficiently. First, generic algorithms are extensible because of their parameterization. For instance, the behavior of Algorithm 1.3 can be changed by acting on the simple point definition or the set of constraints. The scope of this algorithm, initially designed to produce homotopic thinnings of binary skeleton, can even be extended further to handle gray-level images and produce gray-level thinnings. From a theoretical point of view, gray-level images can be processed by decomposing them into different sections. The equivalent of detaching a simple point in a binary image is the lowering of a destructible point in a gray-level context [START_REF] Couprie | Topological operators for grayscale image processing[END_REF].

We have been able to produce gray-level skeletons with Algorithm 1.3 by simply replacing the is simple and detach operations by is destructible and lower functors (see Fig. 6). In the case of a 2D regular images on a square grid, this operation is straightforward as a destructible point can also be characterized locally using new definitions of connectivity numbers.

Generic algorithms can thereafter be turned into patterns or canvases [START_REF] Ornellas | Algorithmic Patterns for Morphological Image Processing[END_REF] allowing the implementation of many algorithms sharing a common core. For example Milena implements morphological algorithms like dilation and erosion, reconstructions, etc. depending on the browsing strategy. Digital Topology could also benefit from a canvas-based approach. The framework can also be extended with respect to data structures. Milena provides site sets based on boxes, graphs and complexes, but more can be added to the library (e.g. combinatorial maps, orders, etc.) and benefit from existing algorithms and tools.

Finally, our approach can take advantage of properties of input types (regularity of the site set, isotropic adjacency relationship, etc.) and allow users to write specialized versions of their algorithms for such subsets of data types, leading to faster or less memory-consuming implementations [START_REF] Levillain | Une approche générique du logiciel pour le traitement d'images préservant les performances[END_REF]. tem@tic Paris-Région" Cluster (France). This project is partially funded by the French Government, its economic development agencies, and by the Paris-Région institutions.

  (a) A simplicial 3-complex, composed of 0-faces (vertices), 1-faces (edges), 2faces (triangles) and a 3-face (tetrahedron). A cubical 2-complex, composed of 0-faces (vertices), 1-faces (edges) and 2-faces (squares).

Fig. 1 .

 1 Fig. 1. Examples of cell complexes.

  (a) 2D binary image. (b) Skeleton of (a) with no constraint (c) Skeleton of (a) where end points of the initial image have been preserved.

Fig. 2 .

 2 Fig. 2. Computation of skeletons from a 2D binary regular image.

  (a)). The following lines produces the result shown on Fig. 2(b).

  (a) 3D binary image. (b) Skeleton (dark voxels) of (a) with no constraint superimposed on the initial image (light voxels).

Fig. 3 .

 3 Fig. 3. Computation of the skeleton of 3D binary regular image.

Fig. 4 .

 4 Fig. 4. Computation of a skeleton using breadth-first thinning. The triangle mesh surface 4(a) (35,286 vertices and 70,568 triangles) is seen as a simplicial 2-complex.The image of curvature 4(b) is computed on the edges of the mesh, and simplified using an area opening filter. All curvature regional minima are then removed from the mesh, and the skeleton 4(c) is obtained with Algorithm 1.3 using the collapse operation.

Fig. 5 .

 5 Fig. 5. Thin skeleton obtain by 2-and 1-collapse.

  (a) 2D gray-level image. (b) Gray-level skeleton.

Fig. 6 .

 6 Fig. 6. Computation of a gray-level skeleton.

  (a) with the following lines.

	typedef image3d<bool> I;
	typedef neighb3d N;
	I output =
	breadth_first_thinning(input,
	c26(),
	is_simple_point3d<I, N>(c26(), c6()),
	detach_point<I>(),
	no_constraint());

Table 1 .

 1 Execution times of Algorithm 1.3 for various inputs. Figures correspond to the time spent in the breadth first thinning routine only.

In Algorithm 1.3, mln ch value(I, V) and mln concrete(I) are helper macros. The former returns the image type associated to I where the value type has been set to V. The latter returns an image type corresponding to I with actual data storage capabilities. In many cases, mln concrete(I) is simply equal to I.

Note that the notion of "constraint" is not the same in Algorithm 1.1 and Algorithm 1.3: in the former, it is the set of points to preserve, while in the latter is it a predicate that a candidate point must pass to be removed.

Acknowledgments The authors thank Jacques-Olivier Lachaud, who reviewed this paper, for his valuable comments, as well the initial reviewers from the WADGMM workshop.

This work has been conducted in the context of the SCRIBO project (http: //www.scribo.ws/) of the Free Software Thematic Group, part of the "Sys-