Image Anal Stereol 2011;30:101-109 Original Research Paper FAST COMPUTATION OF ALL PAIRS OF GEODESIC DISTANCES Guillaume Noyel, Jesus Angulo and Dominique Jeulin MINES ParisTech, CMM - Centre de Morphologie Mathematique, Mathematiques et Systemes, 35 rue Saint Honore - 77305 Fontainebleau cedex, France e-mail: guillaume.noyel@ensmp.fr, jesus.angulo@ensmp.fr, dominique.jeulin@ensmp.fr (Accepted May 15, 2011) ABSTRACT Computing an array of all pairs of geodesic distances between the pixels of an image is time consuming. In the sequel, we introduce new methods exploiting the redundancy of geodesic propagations and compare them to an existing one. We show that our method in which the source point of geodesic propagations is chosen according to its minimum number of distances to the other points, improves the previous method up to 32 % and the naive method up to 50 % in terms of reduction of the number of operations. Keywords: all pairs of geodesic distances, fast marching, geodesic propagation. INTRODUCTION An array of all pairs of geodesic distances, between nodes of a graph, is very useful for several applications such as clustering by kernel methods or graph-based segmentation or data analysis. However, computing this array of distances is time consuming. That is the reason why two new methods, that fill in a fast way the distances array, are presented and compared in this paper. Our methods are available on general graphs. In this paper we present them on images of which the pixels are considered as the nodes of the graph and the links between the neighbors corresponds to the edges of the graph. In image processing, it is interesting to use the geodesic distance between pixels in place of another distance. The array of all pairs of geodesies distances allows to segment an image in geodesically connected regions (i.e., geodesic balls, Noyel e/^aZ, 2007b). All pairs of geodesic distances are also usefial on defining adaptive neighborhoods of filters used for edge-preserving smoothing (Bertelli and Manjunath, 2007; Lerallut et al, 2007; Grazzini and Soille, 2009). Nonlinear dimensionality reduction techniques are mostly based on multidimensional scaling on a Gram matrix of distances between the pairs of variables. Particularly interesting for estimating the intrinsic geometry of a data manifold is the Isometric feature mapping Isomap (Tenenbaum, 1997; Tenenbaum e? a/., 2000). After defining a neighborhood graph of variables, Isomap calculates the shortest path between every pairs of vertices, which is then low-dimensional embedded via multidimensional scaling. The application of Isomap to hyperspectral image analysis requires the computation of all pairs of geodesic distances for a graph of all the pixels of an image (Mohan et al, 2007). However, the computation of all pairs of geodesic distance in an image is time consuming. In fact, the naive approach consists in repeating N times the algorithm to compute the geodesic distance from each of the N pixels to all others. Computing the geodesic distance from one pixel to all others is called a geodesic propagation. The pixel at the origin of a propagation is called the source point and the array containing all pairs of geodesic distances is named the distances array. Therefore, the naive approach is of complexity 0{Nx M), with M the complexity of a geodesic propagation. Several algorithms for geodesic propagations are available: the most famous is the Fast Marching Algorithm introduced by Sethian (1996; 1999) and of complexity 0{Nlog{N)). This algorithm consists in computing geodesic distances in a continuous domain, using a first order approximation, to obtain the distance in the discrete domain. Another algorithm was developed by Soille (1991) for binary images and for grey level images (Soille, 1992). In order to compare our results, we will use the Soille's algorithm of "geodesic time fianction" (Soille, 1994; 2003). Recent implementations of Soille's algorithm for binary images are in 0{Nlog{N)) (Coeuijolly ez-aZ, 2004). Bertelli e?a/. (2006) have introduced a method to exploit the redundancy when several geodesic propagations are computed. In fact, when we perform a geodesic propagation from one pixel, all the geodesic paths from this pixel are stored in a tree. Using this tree, we know all the geodesic distances between any pairs of points along the geodesic path connecting two points. This redundancy is also exploited in earlier algorithm for computing the propagation function well known in mathematical morphology (Lantuejoul and Maisonneuve, 1984). In order to choose the source points of the geodesic propagations, Bertelli et al. (2006) have proposed to select them randomly in a spiral like order starting from the edges of the image and going to the center. In the sequel, we test several deterministic approaches to select the source points and we show that a method based on the filling rate of the distances array can reduce the number of operations up to 32 % compared to Bertelh et al. (2006) method. After discussing some prerequisites about the definition of a graph on an image, the geodesic distances and the exploitation of redundancy between geodesic propagations, we introduce several methods of computation of all pairs of distances and we compare them. PREREQUISITES An image / is a discrete fianction defined on the finite domain E C N^, with N the set of positive integers. The values of a gray level image belongs to C M. For a color image {i.e., with 3 channels) the values are in ^^ = ^ X ^ X ^ and for a multivariate image of L channels the values are in In what follows, we consider ,9' C M+. The whole results presented in the current paper are directly extendable to color or multivariate images (Noyel et al. 2007a, 2007b), (Noyel, 2008). An image is represented on a grid on which the neighborhood relations can be defined. Therefore, an image is seen as a non oriented graph G= {Vq,Eq} in which the vertices Vq correspond to the coordinates of pixels, Vg G Z X Z, and the edges, Eq g I?, give the neighborhood relations between the pixels. Then, the notion of neighborhood of a pixel p in the grid is introduced as the set of pixels which are directly connected to it: \/p,qeVc, p and q are neighbors ^{p,q)GEc, (1) where the ordered pair {p, q) is the edge which joins the points p and q. We assume that a pixel is not its own neighbor and the neighboring relations are symmetrical. The neighborhood of pixel p, Nq{p), defined a subset of Vq of any size, such as: Vp, q G Kg, Ng{p) = {g g Kg, (p, q) G Eg} . (2) Usually in image processing, the following neighborhoods are defined: 4-neighborhood, 8-neighborhood or 6-neighborhood. In the sequel, we use the 8-neighborhood. For our study, the choice of the neighborhood has no influence since we compare several methods using for each one the same neighborhood. A path between two points x and / is a chain of points (xo, ,..., Xj,..., X/) G £ such as xq = xand X] = y, and for all 1, are neighbours. Therefore, a path (xo, ,..., xj) can be seen as a subgraph in which the nodes corresponds to the points and the edges are the connections between neighbouring points. The geodesic distance or geodesic time, between two points of a graph, xq and xj, is defined as the minimum distance, or time, between these two points, The geodesic path ^geo is one of the paths linking these two points with the minimum distance, ^geo{xo,x]) = {xo,... ,x]): ^gso{xo,xi) = {xo,...,xi) (3) such as c/geo(xo,Z/) = inf{tj^{xo,xi)} If the edges of the graph are weighted by the distance between the nodes, t{xi, x/), the geodesic path is one of the sequences of nodes with the minimum weight. To generate a geodesic distance, several measures of dissimilarities can be considered between two neighbors pixels of position xj and x/ and of positive grey values f{xi) and f{x/): - Pseudo-metric 11: dL,{xi,x,,) = \f{xi)-f{x^)\. (4) - Pseudo-metric sum of grey levels: d+{xi,x/) = f{xi) + f{x/) . (5) - Pseudo-metric mean of grey levels (similar to the previous pseudo-metric): d^ixi,x/) = f{xi) + f{xi) (6) The corresponding distances along a path ^ = (xo,... ,xi) are the sum of the pseudo-metric along this path t t,g, {xo ,xi) = ^ d{xi^ 1, xi) . (7) i=l The geodesic distance is defined as the distance along the geodesic path. It is also the minimum of distances over all paths connecting two points xq and xf. dg<,o{xo,Xi) = ^{d{Xi^i,Xi)\Xi^i,Xi G ^geo} i=l = mf{t,^{xo,xi)} . Sr In this paper, we only use the pseudo-metric sum of grey levels c/+ which presents the advantage (compared to d^) not to be null when two pixels have the same strictly positive value (if f{xj) = f{xj) > 0 then dj^{xi,x/) = 2f{xj) > 0). The distance associated to the pseudo-metric sum of grey levels is defined as: 1 i=l i=l (9) = f{x(^) + f{xi)+2 i=i In order to reduce the number of geodesic propagations, Bertelli et al. (2006) have used the following properties. In order to compute the distances between the points along the geodesic path, Bertelli et al. (2006) proposed to build a geodesic tree which has three kinds of nodes: 1. the root, which is the source point for a geodesic propagation. Its distance is null and it has no parents; 2. the nodes, which are points having both parents and children; 3. the leaves, which are points without children. Starting from the leaves to the nodes (or the opposite), the distances between points belonging to the same geodesic path are easily computed. The following additional property is very usefial to compute all pairs of geodesic distances. Property 1 Given a geodesic path (po, Pi, • • •, Pn), the geodesic distance, dg^oiPhPj)' between two points pj et pj along the path, i < J, is equal to the difference dgUPo,Pj)-dgUPo,Pi)- In Fig. 1, if C and D belong to a geodesic path connecting A and B, the geodesic distance between C and D is equal to: 8 =20-12 ^ ' with = 8, = 20, and c/geo(^, Q = 12. (8|>-—--------jl^----------------- 1)— Property 2 The longer the geodesic paths are, the higher numbers of pairs of geodesic distances are computed. In order to have the benefit of the property 2, Bertelli et al. (2006) have chosen as sources, of the geodesic propagations, random points in a spiral-like order: starting from the edges of the image and going to the center. Indeed, the points on the edges of the image tends to have longer geodesic paths than points located at the center. Consequently, we want to test their remark by comparing their method to some others. In order to make this comparison, we measure the filling rates of the distances array D. For an image containing N pixels, the distances array is a square matrix of size NxN=J\P' elements. By symmetry, the number of geodesic paths to compute is equal to: A = N^-N (11) Fig. 1. Discrete geodesic path on a 8-neighborhood graph. The points C and D belong to a geodesic path between A et B. Therefore the geodesic distance between C and D is also known. Using the property 1, when the geodesic distance between two points of the image is computed, the distances between all pairs of points along the associated geodesic path are known. Consequently, the distances array is filled faster using this redundancy. The number of computed paths a is determined by counting the unfilled elements of the distances array D. In practice, it is usefial to use a boolean matrix Z^mrk, of size Nx N, with elements equal to 1 if the distance between the pixel located by the line number and the pixel located by the column number is computed, and 0 otherwise. By convention in the algorithm, we impose to each element of the diagonal of the boolean matrix to be equal to one. Vi : L'mrk(^W) = 1, because the distance from one pixel to itself is equal to zero. Due to the symmetric properties of array D, the number of computed paths is equal to: ^=2 / N N \ X X 1) - trace{D^,^) / N N \ l\k=U=l -N (12) Consequently, the filling rate of the distances array is defined by: The proportion of paths to compute for a given pixel Xi, is named the filling rate of the point Xj, and is defined by: = A^-l A^-l (14) When all the distances from one point to the others are computed, this point is said to be "filled", i.e., T(X,) = 1. INTRODUCTION OF NEW METHODS FOR FAST COMPUTATION OF ALL RMRS OF GEODESIC DISTANCES In order to use homogeneous measures for all methods, the Soille's algorithm (Soille, 2003), called "geodesic time fianction" is used, with a discrete neighborhood of size 3x3 pixels. In fact, we do not need an Euclidean geodesic algorithm to make this comparison study. The Euclidean version is described in (Soille, 1992) and improved in (Coeuijolly e? a/., 2004). SPIRAL METHOD Bertelli et al. (2006) affirms that the source points of the geodesic propagations, with the longest paths, are in general situated on the borders of the image. As these points are usefial to reduce the number of operations, the source points are chosen in a random way on concentric spiral turns of image pixels. A concentric spiral turn is a frame, of one pixel width, in which the top left comer is at position (1,1) or (2,2) or (3,3) or etc. The figure 3 gives an example. While not all the pixels of the spiral turns have been selected, we draw one pixel, among them, in a uniform random way ; otherwise we switch to the next spiral turn. Fig. 3. The spiral turns of an image 9x9 pixels. In the current section, we initially present Bertelli et al. (2006) method, named the "spiral method", and then we introduce two new methods before making comparisons : 1) a geodesic extrema method and 2) a method based on the filling rate of all distances pairs array. The empirical comparisons are made on three different images of size 25 x 25 pixels: "bumps", "hairpin bend", "random" (Fig. 2). Image "bumps" Image "hairpin bend" Image "random Fig. 2. Images of size 25 x 25 "hairpin bend " and "random between 1 and 255. pixels "bumps", whose grey levels are Table 1. Algorithm: Spiral method 1: Given L'the distances array of size Nx N 2: while D is not filled do 3: Select the most exterior spiral turn not yet filled 4: Determine 5 the list of pixels of the spiral turn not yet filled 5: while 5 is not empty do 6: Select a source point 5 randomly in 5 7: Compute the geodesic tree from 5 8: Fill the distances array D 9: Remove the points of 5 which are filled 10: end while 11: end while For each image, the filling rate is plotted versus the number of propagations (fig. 4). TTie number of propagations which are necessary to fill the distances array by the spiral method and the naive method are also given in this figure. The relative difference of the number of propagations of the spiral method compared to the number of propagations of the naive method is written A^(naive). We notice that the spiral method reduces the number of propagations by a factor ranging between 13.4 % and 25.6 %, as compared to the naive one. Consequently, it is very useful to exploit the redundancy in the propagations by building a geodesic tree. Filling rate ^^— ? 465 / / ! / Naive method ^ Spiral method f. 54J/625 1 ! Naive method ^ Spiral method 8 is 0 100 200 300 400 500 600 Number of propagations (a) "bumps" A^(naive) =25.6% Filling rate 0 100 200 300 400 500 600 Number of propagations (b) "hairpin bend" A^(naive) = 13.4 % —yf ^^ 52^625 / 1 : / ^ Naive method - Spiral method 0 100 200 300 400 500 600 Number of propagations (c) "random" A^(naive) = 16.2 % Fig. 4. Comparison of the ßlling rates T of the distances array, between the naive method (in green) and the spiral method (in blue) for the images "bumps ", "hairpin bend " and "random ". The relative differences Ar{naive) of the number of propagations of the spiral method compared to the number of propagations of the naive method are given on the bottom line. SPIRAL METHOD WITH REPULSION Two neighbours points have a high probability to have similar geodesic trees. Consequently, a first improvement of the spiral method is to introduced a repulsion distance between the points drawn randomly in a spiral like-order (algorithm of table 2). Several tests have shown us that a repulsion distance of three pixels on both sides of a source point gives the best filling rates. These tests are empirical tests. In fact, several repulsion distances were tried and it has been noticed that the value of three pixels gives the best results in order to fill the array of distances. This value of three pixels is certainly related to the image size, because, generally, the farther the source points of the geodesic propagations are the faster the array of all pairs of geodesic distances is filled. Table 2. Algorithm: Spiral method with repulsion 1: Given L'the distances array of size Nx N 2: Given h the repulsion distance: h pixels 3: while is not filled do 4: Select the most exterior spiral turn not yet filled 5: Determine 5 the list of pixels of the spiral turn not yet filled 6: while 5is not empty do 7: Select a source point 5 randomly in 5 8: Remove in the list 5 the two left points of 5 and the two right points of 5 if they are still in 5 9: Compute the geodesic tree from 5 10: Fill the array D 11: Remove the points of 5 which are filled 12: end while 13: end while Fig. 5 shows that the relative differences in the number of propagations necessary to fill the distances array are larger than 15 % for images "bumps" and "hairpin bend" and than 5.3 % for the image "random". Therefore, the spiral method with the repulsion distance is faster than the spiral method to fill the distances array. GEODESIC EXTREMA METHOD As the longest geodesic paths are those which fill the most the distance array, we look for the geodesic extrema of the image. To compute the geodesic extrema of the image, we use two geodesic propagations: - a first propagation starts from the edges of the image. Then we select one of the points c with the longest distance from the edges. This point is called the geodesic centroid ; - a second propagation from the geodesic centroid gives the farthest points from the centroid, i.e., the geodesic extrema of the image. On the figure 6, which shows the results from these two propagations, we notice that the geodesic extrema are mainly on the edges of the image, which ascertains the motivation to use a spiral method. g o Filling rate 1-O Fllllng rate «5 460 541 ^oo £ 1» 1° r eO" S f? = r 1 ^ Spiral method f repulsion 3 » Spiral method o- 1 » Spiral method ^ repulsion 3 Spiral method 0 100 200 300 400 500 600 Number of propagations 0 100 200 300 400 500 600 Number of propagations (a) "bumps" (b) "hairpin bend" Ar(spiral) = 17.4% Ar(spiral) = 15.0% Ar(naive) = 38.6% Ar(naive) = 26.4% g s Filling rate ^^ 4% 524 t^ 1° s 1 ^ Spiral method f repulsion 3 » Spiral method 0 100 200 300 400 600 600 Number of propagations (c) "random" A^(spiral) = 5.3% A^(naive) = 20.6% Fig. 5. Comparison of the ßlling rates T of the distances array, between the spiral method with repulsion (in red) and the spiral method (in blue) for the images "bumps ", "hairpin bend " and "random ". The relative differences Aj-{spiral) (resp. Aj-{naive)) of the number of propagations of the spiral method with repulsion compared to the number of propagations of the spiral (resp. naive) method are given on the bottom lines. Then the pixels are sorted by geodesic distance from the centroid into the list of geodesic extrema Ext. This list is used to choose the source points of the geodesic propagations. If several points have the same geodesic distance from the centroid, then the less filled is selected (algorithm of table 3). The filling rates of the spiral method and the geodesic extrema method versus the number of propagations are plotted for each image (fig. 7). We notice that the filling rates of the geodesic extrema method are at the beginning inferior or similar to these of the spiral method. However, at the end the filling rates of the geodesic extrema method are better than these of the spiral method. In fact, we are looking for an approach filling totally the distances array in the fastest way. As the number of propagations of the geodesic extrema method necessary to fill the distances array are lower than in the spiral method, the geodesic extrema method fills faster the distances array than the spiral method. "bumps" "hairpin bend" "random" Geodesic distance from the edges of the image Geodesic distance from the centroid (red point) Fig. 6. Geodesic distance from the edges of the image (second line) and from the centroid in red (third line) for several images. Table 3. Algorithm: Geodesic extrema method Given L'the distances array of size Nx N Compute the geodesic centroid c Compute the decreasing list of geodesic extrema Ext whose first element Ext[l\ is the greatest geodesic extrema not yet filled 4: Initialise to zero, the list, of size N, of the filling rate of points T. 5: while D is not filled do 6: B^ {xExt\d^^o{x, c) = d^^o{Ext[\],c)} 7: s*— x 8: Compute the geodesic tree from 5 9: Fill the distances array D 10: Remove the points of the list Ext which are filled 11: Update the list of filling rates T 12: end while Filling rate I fsi r 437 465 / / f Extrema method ^ Spiral method Filling • Extrema method ^ Spiral method 0 100 200 300 400 500 600 Number of propagations (a) "bumps" A^(spiral) = 6.0% A^(naive) =30.1% Filling rate 0 100 200 300 400 500 600 Number of propagations (b) "hairpin bend" A^(spiral) = 5.2% A^(naive) = 17.9% Ž ^---- / ^^ 506524 / 't - Extrema method - Spiral method 0 100 200 300 400 500 600 Number of propagations (c) "random" A^(spiral) =3.4% A^(naive) = 19.0% Fig. 7. Comparison of the ßlling rates T of the distances array, between the geodesic extrema method (in red) and the spiral method (in blue) for the images "bumps ", "hairpin bend " and "random ". The relative differences Ar{spiraf) (resp. Ar{naive)) of the number of propagations of the the geodesic extrema method compared to the number of propagations of the spiral (resp. naive) method are on the bottom lines. In order to get an exact comparison it is necessary to generate two propagations, corresponding to the determination of the geodesic extrema, to the number of propagations necessary to fill the distances array. Even, with this modification, the extrema method fills the distances array faster than the spiral method (the relative differences A^(spiral) are between 3.4 % and 6 %). METHOD BASED ON THE FILLING RATE OF THE DISTANCES ARRAY In place of selecting the source points from their distance from the geodesic centroid propagation, we select first the less filled points. To this aim, after each geodesic propagation the filling rate is computed for each point. Then the less filled point is selected as a source of the propagation. If several points are among the less filled, then the greatest geodesic extrema is chosen among these points (algorithm of table 4). Table 4. Algorithm: Method based on the filling rate of the distances array. 1: Given L'the distances array of size Nx N 2: Compute the decreasing list of geodesic extrema Ext 3: Initialise to zero, the list, of size N, of the filling rate of points T. 4: while is not filled do 5: {argmin^££T[;s^} 6: if Card{B} > 1 then 7: 5 ^ 8: else 9: 5 ^ argmin^g T [x] 10: end if 11: Compute the geodesic tree from 5 12: Fill the distances array D 13: Remove the points of the list Ext which are filled 14: Update the list of filling rates T 15: end while As for the geodesic extrema method, we compare the filling rates of the method based on the filling rate of the distances array and the spiral method versus the number of propagations (fig. 8). We notice that the method based on the filling rate of D reduces of 32.3 % the number of propagations of the spiral method on the image "bumps" and 18.7 % on the image "hairpin bend". Consequently this method fills faster the distances array than the spiral method. Even for the image "random" the method based on the filling rate of D still improves the spiral method of 2.7 %. However it is not very common to compute all pairs of geodesic distances on a strong unstructured image such as the "random" one. By comparison to the naive approach, the method based on the filling rate of D reduces the number of propagations by a 49.6 % rate (resp. 29.6 %) on the image "bumps" (resp. "hairpin bend"). As for the previous method, in order to get an exact comparison, it is necessary to add two propagations, corresponding to the determination of the geodesic extrema, to the number of propagations necessary to fill the distances array. Even, with this modification, the number of propagations necessary to fill the distances array is still lower for the method based on the filling rate of D. Filling rate --- 465 if ij 1 ^ Fiiiing rate metiiod » Spirai method Filling rate 0 100 200 300 400 500 600 Number of propagations (a) "bumps" Dr (spiral) = 32.3% Ar(naive) = 49.6% Filling rata 100 200 300 400 500 600 Number of propagations (b) "hairpin bend" Ar(spiral) = 18.7% Ar (naive) = 29.6% 1° ü / yT^ 510 524 f • Fiiiing rate method » Spirai method 0 100 200 300 400 500 600 Number of propagations (c) "random" Ar(spiral) = 2.7% Ar(naive) = 18.4% Fig. 8. Comparison of the filling rates t of the distances ar^ray, betw^een the method based on the filling rate of the distances ar^ray (in red) and the spiral method (in blue) for the imagoes "bumps ", "hairpin bend " and "random ". The relative differences Ar{spiral) (resp. Ar(naive)) of the number ofpropag^ations of the method based on the filling rate compared to the number of propagations of the spiral (resp. naive) method a^re given on the bottom lines. on the filling rate of D reduces the number of operations between 19 % and 32 %, as compared to the spiral approach and between 30 % and 50 %, as compared to the naive method, on standard images. Even on "random" image the improvements is of 3 % (resp. 18 %) compared to the spiral (resp. naive) method. Consequently, the filling rate of the distances array, combined with the geodesic extrema when several points have the minimum filling rate, seems to be the best criterion to fill efficiently the distances array. CONCLUSION From a comparison between different approaches, it turns out that a method based on the optimization of the filling rate of the distances array is the most efficient to compute the geodesic distances between all pairs of pixels in an image. Besides, in the current paper, we have shown our results on grey images. They can be directly extended to hyperspectral images using appropriate pseudo-metrics (Noyel et al., 2007a). This can be a useful step for a subsequent clustering by kernel methods or data reduction approaches on multivariate images. Filling rate Filling rate DISCUSSION After having presented and tested several methods to fill the distances array, we have compared them for the three test images "bumps", "hairpin bend" and "random" on the figure 9 and in the table 5. For the images "bumps" and "hairpin bend", the method based on the filling rate of the distances array is faster than the other algorithms. For the "random image" (an extreme case presenting no texture) the methods introduced here give similar results, since the relative difference between the maximum number of propagations is less than 2.7 %. Therefore, we conclude that the method based on the filling rate of the distances array is the best one to calculate the array of all pairs of distances. According to their performances, the others are ranked in the following order: 1) the spiral method with repulsion distance, 2) the geodesic extrema method and 3) the spiral method ofBertelli eta^^. (2006). Moreover, we have shown that the method based 184 465 315 437 spiral geodesic extreme fiiiing rate spiral rspuision 3 0 100 300 500 Number of propagations (a) "bumps" Fiiiing r spirai geodesic extrema fiiiing rate spirai repulsion 3 0 100 300 500 Number of propagations (c) "random" 0 100 300 500 Number of propagations (b) "hairpin bend" Fig. 9. Comparison of ^e filling rates, of the array of distances, for the methods: spiral, geodesic extrema, the method based on the filling and the spiral method w^ith a repulsion distance of 3 pixels for the images "bumps ", "hairpin bend " and "random ". Table 5. Relative difference values of filling rates (a) between the different methods and the native approach, or (b) betw^een the different methods a^nd the spiral approach. For each image is given in bold the best relative difference value. (a) Spiral spiral geodesic filling Ar (naive) method method extrema rate with method method repulsion "Bumps" 25.6 % 38.6 % 30.1 % 49.6 % "Hairpin bend" 13.4 % 26.4 % 17.9 % 29.6 % "Random" 16.2 % 20.6 % 19.0 % 18.4 % (b) spiral geodesic filling Ar (spiral) method with extrema rate repulsion method method "Bumps" 17.4 % 6.0 % 32.3 % "Hairpin bend" 15.0 % 5.2 % 18.7 % "Random" 5.3 % 3.4 % 2.7 % The main motivation for our developments is the computation of all pairs of geodesic distances for the pixels of an image, which is usually a graph of thousands of vertices arranged spatially. Nevertheless, our approach is valid on more general graphs, than those associated to bitmap images, after determining the "boundary vertices" of the graph, since then, the computation of the geodesic centre (and the geodesic extremities) of a graph can be obtained by a first propagation from the "boundary vertices". The "boundary vertices" can be defined, for instance, as the vertices having less neighbouring vertices than the average number of connectivity in the graph. REFERENCES BertelliL, SumengenB,ManjunathBS (2006). Redundancy in All Pairs Fast Marching Method. In: Proc IEEE Int Conf Image Process ICIP'06. 3033-6. Bertelli L, Manjunath BS (2007). Edge Preserving Filters using Geodesic Distances on Weighted Orthogonal Domains. In: Proc IEEE Int Conf Image Process ICIP'07. I:321-4. Coeurjolly D, Miguet S, Tougne L (2004). 2D and 3D visibility in discrete geometry: an application to discrete geodesic paths. Pattern Recogn Lett 25:561-70. Grazzini J, Soille P (2009). Edge-preserving smoothing using a similarity measure in adaptive geodesic neighbourhoods. Pattern Recogn 42(10):2306-16. Lantuejoul C, Maisonneuve F (1984). Geodesic methods in image analysis. Pattern Recogn 17:177-87. Lerallut R, Decenciere E, Meyer F (2007). Image filtering using morphological amoebas. Image Vision Comput 25(4):395-404. Mohan A, Sapiro G, Bosch E (2007). Spatially Coherent Nonlinear Dimensionality Reduction and Segmentation of Hyperspectral Images. IEEE Geosci Remote Sens 4(2):206-10. Noyel G, Angulo J, Jeulin D (2007a). Morphological segmentation of hyperspectral images. Image Anal Stereol 26:101-9. Noyel G, Angulo J, Jeulin D (2007b). On distances, paths and connections for hyperspectral image segmentation. In: Banon G et al. Proc 8th Int Symp Math Morpho 1:399-410. Noyel G (2008). Filtrage, reduction de dimension, classification et segmentation morphologique hyperspectrale. PhD Thesis. Mines ParisTech, France. Sethian JA (1996). A marching level set method for monotonically advancing fronts. Proc Natl Acad Sci USA 93:1591-5. Sethian JA (1999). Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision and Materials Science. Cambridge: Cambridge University Press. Soille P (1991). Spatial distributions from contour lines: an efficient methodology based on distance transformations. J Vis Commun Image Rep 2(2):138-50. Soille P (1992) Morphologie Mathematique: du relief a la dimensionalite - algorithmes et methodes. PhD Thesis. Universite Catholique de Louvain. Soille P (1994). Generalized geodesy via geodesic time. Pattern Recogn Lett 15:1235-40. Soille P (2003). Morphological image analysis, 2nd Ed. Berlin, Heidelberg: Springer-Verlag. Tenenbaum JB (1997). Mapping a manifold of perceptual observations. Adv Neural Inf Process Syst 10:682-8. Tenenbaum JB, de Silva V, Langford JC (2000). A global geometric framework for nonlinear dimensionality reduction. Science 290:2319-23.