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

Stochastic permutation ordering watershed

Keywords: Segmentation, stochastic watershed, Hamiltonian path

The stochastic watershed is a morphological approach to segmentation that repeats the application of a seeded watershed from series of uniform random markers. The obtained watershed boundaries are combined to construct a probability density function. We propose an alternative approach called stochastic permutation ordering watershed. Our approach relies on the construction of several permutation orderings of the pixels of an image, starting from a random pixel and using a dedicated Hamiltonian path construction on a graph. From the permutation orderings, several seeded watersheds are applied and averaged. In contrast to the stochastic watershed, our approach enables to take into account any features associated to pixels, such as patches, of prime important to segment textured images. Experimental results show the benefit of the approach.

I. INTRODUCTION

For image segmentation, the watershed [START_REF] Meyer | Morphological segmentation[END_REF] is among the most popular methods that have been proposed so far. This approach belongs to the research field of Mathematical Morphology (MM) [START_REF] Soille | Morphological Image Analysis: Principles and Applications[END_REF], a nonlinear approach to image processing. Starting from region seeds, the watershed performs a flooding on a gradient image to obtain a segmentation. Efficient algorithms [START_REF] Vincent | Watersheds in digital spaces: an efficient algorithm based on immersion simulations[END_REF] have been proposed to perform this flooding process. In an unsupervised segmentation process, region seeds are usually the minima of the gradient. However this is very sensitive to both noise and over-segmentation [START_REF] Meyer | Morphological segmentation[END_REF]. Typical approaches to cope with these problems do either remove local minima in a pre processing step or merge the final regions in a post processing step. Another solution is to consider, instead of the gradient image, an image that reduces the minima inside the objects and enhances the boundaries. With this objective, Angulo and Jeulin have proposed [START_REF] Angulo | Stochastic watershed segmentation[END_REF] the stochastic watershed. The principle is twofold. First, one repeats the application of a seeded watershed from series of uniform random markers using Monte Carlo simulations. Second, the obtained boundaries of the watershed segmentation are combined to construct a probability density function (pdf) for the boundaries in the image. The pdf computed by the stochastic watershed converges to the probability of seeds falling on both sides of each boundary and can be computed using graph algorithms [START_REF] Malmberg | An efficient algorithm for exact evaluation of stochastic watersheds[END_REF]. Several works have been proposed to extend the stochastic watershed. In [START_REF] López-Mir | Probability density function of object contours using regional regularized stochastic watershed[END_REF] a variant based on regional regularization of the image was proposed. In [START_REF] Franchi | Bagging stochastic watershed on natural color image segmentation[END_REF] a multiscale approach including bagging was proposed. In [START_REF] Angulo | Statistical gaussian model of image regions in stochastic watershed segmentation[END_REF], the probabilities of contours are computed from a gaussian model of image regions.

In this paper we propose an alternative approach for the computation of the probability density function for the boundaries in the image. In the original formulation of [START_REF] Angulo | Stochastic watershed segmentation[END_REF], each simulation is performed on a gradient image [START_REF] Zenzo | A note on the gradient of a multi-image[END_REF] that usually only considers the color differences between pixels. Therefore, the classical stochastic watershed cannot use more advanced pixel features such as patches to construct the pdf for the boundaries in the image. However, texture features such as patches have shown to be very efficient in many computer vision tasks [START_REF] Buades | Image denoising methods. a new nonlocal principle[END_REF]. We introduce an approach to build a pdf for the boundaries in the image that can consider any feature describing a pixel. Our approach relies on the construction of several permutation orderings of the pixels of an image, starting from a random pixel and using a dedicated Hamiltonian path construction on a graph. In section II we present how to construct such a stochastic Hamiltonian path and construct a pdf. In Section III, we present some results and comparisons. Last section concludes.

II. STOCHASTIC PERMUTATION ORDERING WATERSHED

A. Notations

We consider that the domain Ω of the image is a graph G = (V, E) where vertices V = {v 1 , . . . , v m } correspond to pixels and edges e ij = (v i , v j ) connect vertices. With this definition, images are represented as graphs signals [START_REF] Lézoray | Manifold-based mathematical morphology for graph signal editing of colored images and meshes[END_REF] that associate vectors to vertices and are defined by the mapping f : G → F ⊂ R n where F is a non-empty set of vectors (we will consider only RGB color vectors, i.e., n = 3). To each vertex v i ∈ G is associated a vector

v i = f (v i ). The set T = {v 1 , • • • , v m }
denotes the vectors associated to all the vertices of the graph. Consequently, one has |V| = m. We will use the notation T [i] = v i to denote the i-th element of the set T .

The graph G used to represent an image can be a classical grid-graph where pixels are connected with 8-adjacency, but it can be interesting to go beyond this classical grid-graph construction and to consider graphs of the higher connectivities. We will consider several possible graphs to represent an image. The graph can use only spatial or spectral information, or both.

• B-adjacency graph (denoted G B ): this graph construction connects each vertex v i to all the vertices contained in a square box of size (2B + 1) × (2B + 1) around v i . 8-adjacency grid graphs do correspond to graphs G 1 . • Global K-Nearest Neighbor graph (denoted G K * ): this graph construction connects each vertex v i to its K nearest neighbors (in terms of spectral distance) within the set of all vertices. The first graphs G B consider only the spatial proximity information to connect vertices whereas graphs G K * consider only the spectral proximity information and are fully nonlocal graphs [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing[END_REF].

B. Complete lattices

The classical construction of algebraic morphological operators relies on complete lattices [START_REF] Ronse | Why mathematical morphology needs complete lattices[END_REF] that impose the need of an ordering relationship between the elements to be processed. Within this theoretical background, morphological operations can be described as mappings between complete lattices. A complete lattice (T , ≤) is a non-empty set equipped with a total ordering relation, such that every non-empty subset P of T has a lower bound ∧P and an upper bound ∨P. If MM is well defined for gray scale functions, there exists no general admitted extension that permits to perform morphological operations on vectors since there is no natural ordering of vectors [START_REF] Barnett | The ordering of multivariate daa[END_REF]. If many approaches have been proposed for color images [START_REF] Aptoula | Multivariate mathematical morphology applied to colour image analysis[END_REF], very few deal with patches [START_REF] Velasco-Forero | On nonlocal mathematical morphology[END_REF], [START_REF] Lézoray | Complete lattice learning for multivariate mathematical morphology[END_REF]. In this paper we will define an ordering relation between the vectors of a set T with the use of h-orderings [START_REF] Goutsias | Morphological operators for image sequences[END_REF]. This corresponds to defining a surjective transform h from T to L where L is a complete lattice equipped with the conditional total ordering [START_REF] Goutsias | Morphological operators for image sequences[END_REF]. We refer to ≤ h as the h-ordering given by:

h : T → L and v → h(v), ∀(v i , v j ) ∈ T × T v i ≤ h v j ⇔ h(v i ) ≤ h(v j ) . (1)
Then, T is no longer required to be a complete lattice, since the ordering of T can be induced upon L by means of h [START_REF] Aptoula | Multivariate mathematical morphology applied to colour image analysis[END_REF]. When h is bijective, this corresponds to defining a space filling curve [START_REF] Chanussot | Bit mixing paradigm for multivalued morphological filters[END_REF] and we will use this property in the sequel.

C. Stochastic Hamiltonian path

To construct a complete lattice with an image-adaptive hordering, we propose to construct a space-filling curve [START_REF] Chanussot | Bit mixing paradigm for multivalued morphological filters[END_REF] on a graph G associated to the image. This corresponds to the construction of an Hamiltonian path on the image: a path that goes through all the vertices of the graph and traverses each vertex only once. We consider a general approach and build an order for all the set T in the form of a global Hamiltonian path.

Given the set T , the construction of an Hamiltonian path amounts to define a sorted permutation P = PT of the vectors of T with P a permutation matrix of size m × m. Let σ be a permutation of the index set

J = {1, • • • , m}. If σ(i) = j, then P ij = 1 and 0 otherwise. The induced permutation is P = PT = {v σ -1 (1) , • • • , v σ -1 (m) } where (σ -1 • σ)
is the identity. Any permutation is not of interest, and spatial and spectral constraints have to be taken into account. To this end we consider the construction of a smooth permutation. The smoothness of an ordered set is expressed by the Total Variation of its elements:

T T V = m-1 i=1 v i -v i+1 (2)
The optimal permutation operator P can be obtained by minimizing the total variation of PT :

P * = arg min P PT T V (3)
This optimization problem provides a permutation such that the corresponding Hamiltonian path is the shortest one and is equivalent to solve the traveling salesman problem. This problem being too computationally demanding for large sets, it can be solved using a randomized version of nearest neighbors heuristics as presented in [START_REF] Ram | Image processing using smooth ordering of its patches[END_REF]. We recall the principle of this algorithm in Algorithm 1. This algorithm starts from an arbitrary vertex and continues by finding the two nearest unexplored neighbor vertices and choosing one of them at random. Both these random choices make the obtained permutation stochastic.

Algorithm 1: Permutation construction as a global Hamiltonian path on a graph.

Input: set of vectors T = {v1, • • • , vm} and graph G; Generate a random probability vector p = [p1, • • • , pm]; Randomly choose an index j ∈ J = {1, • • • , m}; Set σ -1 (1) = j (or equivalently σ(j) = 1); Set P = {v σ -1 (1) } and J = {σ -1 (1)}; for i = 1 to m -1 do Let N (v σ -1 (i) ) = {k|(v σ -1 (i) , v k ) ∈ E} \ I; if |N (v σ -1 (i) )| = 1 then σ -1 (i + 1) = N (v σ -1 (i) ) else if |N (v σ -1 (i) )| ≥ 2 then Find the first vj 1 and second vj 2 nearest neighbors of v σ -1 (i) in N (v σ -1 (i) ) else if |N (v σ -1 (i) )| = 0 then Find the first vj 1 and second vj 2 nearest neighbors of v σ -1 (i) in J \ I else qi = 1 1+exp   v σ -1 (i) -v j 1 -v σ -1 (i) -v j 2 α   ;
If (qi < p σ -1 (i) ) then σ -1 (i + 1) = j2 and σ -1 (i + 1) = j1 otherwise;

P = P ∪ {v σ -1 (i+1) }; I = I ∪ {σ -1 (i + 1)};
Output: ordered set P and associated index set I

After the construction of the permutation, we define the hordering from the constructed permutation as h(v i ) = σ(i) and this defines the complete lattice (T , ≤ h ). Given a graph signal f : G → T , a new representation is obtained in the form of the pair (I, P) with I(v i ) = σ(i). When a graph signal is encoded in this way, the spectral information is not directly carried by the index I, but is stored in a separate piece of data called a palette: the set P of sorted vectors. Figure 1 presents the obtained permutation and index from an original image, as well as the constructed Hamiltonian path, on a 8-adjacency grid graph G 1 . The original graph signal f can be directly recovered since f

(v i ) = P[I(v i )] = T [i] = v i .
It can be interesting to construct graphs that combine the use of both spatial and spectral proximity to infer the set of edges. For graphs G K * , it is even recommended since their construction cannot ensure that the graph is connected. In this special case, Algorithm 1 will not be able to continue the construction of the path in the vertices of the neighborhood if they have been all explored. This case was foreseen in Algorithm 1, and in the case where |N (v σ -1 (i) )| = 0, meaning that no vertex is available to continue the path, the search of the nearest neighbor is performed on the set of all the remaining vertices J \I. These cases can be easily seen as the "jumping" edges in the constructed path shown in Figure 1(b). In the sequel we will consider graphs G 20 ∪ G 21 * . Finally, the spectral proximity used to determine the nearest neighbors within the set of adjacent vertices can consider any distance measure. In particular, instead of using a classical L 2 distance between the vectors v i associated to vertices, one can consider patches

p w i = f (x i + t), ∀t ∈ [-w/2, w/2] 2 T of size w × w.

D. Constructing the probability density function

As Algorithm 1 starts from a random vertex, a different stochastic permutation ordering is obtained for each run of the algorithm. If this is problematic for the construction of morphological operators, we can make the most of it and propose an alternative formulation of the stochastic watershed [START_REF] Angulo | Stochastic watershed segmentation[END_REF]. In this latter approach M segmentations are obtained from N random markers to build a probability density function using a Parzen estimator to combine them altogether. This pdf is then combined with the initial gradient to produce a probabilistic gradient that can be used for watershed segmentation. The stochastic nature of this approach comes from the use of random markers. We propose a novel approach that we will call stochastic permutation ordering watershed. In our case the stochastic aspect will come from the generation of several permutation orderings from different starting vertices. The principle is the following. We built M stochastic permutation orderings using Algorithm 1. This provides M representations 2 presents for an image, several different obtained indexes and permutations. As it can be seen, the orderings of the colors can be very different in the palette (or equivalently the lattice) because of the stochastic nature of the constructed Hamiltonian paths. As a consequence the index palette is also very different from one stochastic permutation ordering to another. Therefore the minima extracted from index images I i will not be located at the same position in all index images and we can obtain different segmentations from these minima, as in the original stochastic watershed. However in our case the germs are not randomly generated but are the minima of the index images generated from stochastic permutation orderings. Figure 2 presents the obtained minima and watershed segmentations for each of the obtained indexes with the Di Zenzo gradient [START_REF] Zenzo | A note on the gradient of a multi-image[END_REF]. For an image f , we denote as W S i (f ) = W S(I i , ∇f ) the watershed obtained from the minima m(I i ) of the ith stochastic ordering I i on the gradient ∇f . Then, as in [START_REF] Angulo | Stochastic watershed segmentation[END_REF], the pdf is computed with the Parzen method:

(I i , P i ) with i ∈ [1, M ]. Figure
pdf(f ) = 1 M M i=0 G(W S i (f )) (4)
with G a Gaussian kernel of variance σ 2 (σ = 2 in our case). Finally this pdf is combined with the initial Di Zenzo Gradient to produce a probabilistic gradient :

∇f = pdf(f ) + ∇f 2 (5)
This enables to reinforce the gradients contours of high probability. First row of Figure 2 shows the obtained probabilistic gradient. This was obtained using 11 combined different stochastic permutations orderings even if only 3 are shown in the Figure . As it can be seen, the probabilistic gradient contains much more information than the classical Di Zenzo gradient. So far, the distance used to generate the different stochastic ordering in Algorithm 1 was based only on the color vectors of the vertices v i = f (v i ). As previously explained, we can used instead a distance on other features such as the patches p w i of size w × w or structure tensors t i of size 2 × 2. This directly enables to have a feature-based stochastic watershed which is not possible with the original formulation of [START_REF] Angulo | Stochastic watershed segmentation[END_REF]. We investigate this in the next Section.

III. EXAMPLES

In this section we investigate the interest of the proposed stochastic permutation ordering watershed with respect to the original stochastic watershed [START_REF] Angulo | Stochastic watershed segmentation[END_REF]. We also investigate the interest of using patches or structure tensors instead of colors to compare vertices in Algorithm 1. We perform our experiments on images of the Bayeux Tapestry 1 . On such images, historian experts need to delineate interactively some characters. Therefore, a precise segmentation is required but by using simple object/background seed labeling by point click to ease the end-users use. The characters are visually easy to identify but the reduced number of colors, the fine embroidery as well as the texture differences in the linen fabric can make the segmentation hard. Figure 3 presents some results on four distinct parts of the tapestry. To obtain the segmentation, a watershed is performed on the probabilistic gradient from the user seeds. We compare our approach with the classical stochastic watershed. To have a fair comparison, for the computation of pdf(f ) we consider the same number of random realizations than ours [START_REF] Lézoray | Manifold-based mathematical morphology for graph signal editing of colored images and meshes[END_REF] and generate 200 random seeds. In contrast to [START_REF] Angulo | Stochastic watershed segmentation[END_REF], we do not add any pre-processing of the image. For each image to segment, Figure 3 presents on two lines, from left to right: the original image f , the user seeds, the probabilistic gradient ∇f and the segmentation result for the classical stochastic watershed, our stochastic permutation ordering watershed with colors v i , 7 × 7 patches p 7 i , 15 × 15 patches p 15 i , and structure tensors t i . As it can be seen, the probabilistic gradient obtained from the classical stochastic watershed (CSW) has few regions in homogeneous areas but some visually important edges do not have a high probability. This is not the case of the probabilistic gradient obtained from our proposed stochastic permutation ordering watershed (SPOW). In the segmentation results, the blue squares do identify the strong segmentation errors. The results obtained with colors vectors for SPOW are always better than for CSW. One can see that for the observed errors this is due to not having edges of high probabilities. The use of 7 × 7 patches can further enhance the segmentation. However the use of very large 15 × 15 patches provides a result close to the CSW proving that if patches can be of interest, the size of the patterns of interest is very important. Finally with structure tensors, the probabilistic gradient is very fine, but this tends to over detect structures in textured areas. As a conclusion, for this kind of images the use of our proposed SPOW with 7 × 7 patches provides the best results and this shows the interest of our proposal.

IV. CONCLUSION

We have proposed an alternative formulation to the stochastic watershed. Starting from a random pixel, an Hamiltonian path is constructed to build a complete lattice. From the latter a new representation of an image in the form of a palette and an index image is obtained. Minima are extracted from the latter and used to perform several watershed that are combined into a probability density function. A probabilistic gradient is computed from the obtained pdf and the Di Zenzo gradient. The advantage of the proposed approach relies in its ability to incorporate any features to represent pixels for the computation of the pdf. The experimental results have shown that this has some benefits with respect to the classical stochastic watershed.

Fig. 1 .

 1 Fig. 1. From left to right: original image, an Hamiltonian path constructed on graph G 1 , the associated index and palette images.

Fig. 2 .

 2 Fig. 2. First row: initial image f , Di Zenzo gradient ∇f , probabilistic gradient ∇f . Second row: stochastic orderings I i . Third row: stochastic palettes P i . Fourth row: minima m(I i ) of the stochastic orderings I i . Fifth row: watershed lines W S i (f ) from the minima m(I i ) on the gradient ∇f .

1

  Panorama of the Bayeux Tapestry, GREYC Image (UMR UNICAEN/ENSICAEN/CNRS), CERTIC (UNICAEN) from the 2017 photographic campaign of La Fabrique de patrimoines en Normandie, 2017-2018.

Fig. 3 .

 3 Fig. 3. Segmentation results, with for each image, on two lines, from left to right: the original image f , the user seeds, the probabilistic gradient ∇f and the segmentation result for the classical stochastic watershed, our stochastic permutation ordering watershed with colors v i , 7 × 7 patches p 7 i , 15 × 15 patches p 15 i , and structure tensors t i .