Fast computation of all pairs of geodesic distances
SSubmitted to Image Anal Stereol, 10 pages
Original Research Paper
FAST COMPUTATION OF ALL PAIRS OF GEODESIC DISTANCES G UILLAUME N OYEL , J ESÚS A NGULO AND D OMINIQUE J EULIN
MINES ParisTech, CMM - Centre de Morphologie Mathématique, Mathématiques et Systèmes, 35 rue SaintHonoré - 77305 Fontainebleau cedex, Francee-mail: {guillaume.noyel, jesus.angulo, dominique.jeulin}@ensmp.fr (Submitted)
ABSTRACT
Computing an array of all pairs of geodesic distances between the pixels of an image is time consuming. Inthe sequel, we introduce new methods exploiting the redundancy of geodesic propagations and compare themto an existing one. We show that our method in which the source point of geodesic propagations is chosenaccording 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, geodesic propagation, fast marching.
INTRODUCTION
An array of all pairs of geodesic distances, betweennodes of a graph, is very useful for several applicationssuch as clustering by kernel methods or graph-basedsegmentation or data analysis. However, computingthis array of distances is time consuming. That is thereason why two new methods, that fill in a fast waythe distances array, are presented and compared in thispaper.Our methods are available on general graphs. Inthis paper we present them on images of which thepixels are considered as the nodes of the graph and thelinks between the neighbors corresponds to the edgesof the graph.In image processing, it is interesting to use thegeodesic distance between pixels in place of anotherdistance. The array of all pairs of geodesics distancesallows to segment an image in geodesically connectedregions (i.e., geodesic balls) (Noyel et al. , 2007b). Allpairs of geodesic distances are also useful on definingadaptive neighborhoods of filters used for edge-preserving smoothing (Lerallut et al. , 2007; Bertelliand Manjunath, 2007; Grazzini and Soille, 2009).Nonlinear dimensionality reduction techniques aremostly based on multidimensional scaling on a Grammatrix of distances between the pairs of variables.Particularly interesting for estimating the intrinsicgeometry of a data manifold is the Isometric featuremapping
Isomap (Tenenbaum, 1998; Tenenbaum et al. , 2000). After defining a neighborhood graphof variables, Isomap calculates the shortest pathbetween every pairs of vertices, which is then low-dimensional embedded via multidimensional scaling.The application of Isomap to hyperspectral imageanalysis requires the computation of all pairs of geodesic distances for a graph of all the pixels of animage (Mohan et al. , 2007).However, the computation of all pairs of geodesicdistance in an image is time consuming. In fact,the naive approach consists in repeating N timesthe algorithm to compute the geodesic distance fromeach of the N pixels to all others. Computing thegeodesic distance from one pixel to all others iscalled a geodesic propagation. The pixel at the originof a propagation is called the source point andthe array containing all pairs of geodesic distancesis named the distances array. Therefore, the naiveapproach is of complexity O ( N × M ) , with M the complexity of a geodesic propagation. Severalalgorithms for geodesic propagations are available:the most famous is the Fast Marching Algorithmintroduced by Sethian (1996; 1999) and of complexity O ( N log ( N )) . This algorithm consists in computinggeodesic distances in a continuous domain, using afirst order approximation, to obtain the distance inthe discrete domain. Another algorithm was developedby Soille (1991) for binary images and for greylevel images (Soille, 1992). In order to compareour results, we will use the Soille’s algorithm of“geodesic time function” (Soille, 1994; 2013). Recentimplementations of Soille’s algorithm for binaryimages are in O ( N log ( N )) (Coeurjolly et al. , 2004).Bertelli et al. (2006) have introduced a methodto exploit the redundancy when several geodesicpropagations are computed. In fact, when we perform ageodesic propagation from one pixel, all the geodesicpaths from this pixel are stored in a tree. Using thistree, we know all the geodesic distances between anypairs of points along the geodesic path connecting twopoints. This redundancy is also exploited in earlieralgorithm for computing the propagation function well1 a r X i v : . [ c s . D M ] J u l OYEL G et al. : Fast computation of all pairs of geodesic distances known in mathematical morphology (Lantuejoul andMaisonneuve, 1984).In order to choose the source points of the geodesicpropagations, Bertelli et al. (2006) have proposed toselect them randomly in a spiral like order startingfrom the edges of the image and going to the center.In the sequel, we test several deterministic approachesto select the source points and we show that a methodbased on the filling rate of the distances array canreduce the number of operations up to 32 % comparedto Bertelli et al. (2006) method.After discussing some prerequisites about thedefinition of a graph on an image, the geodesicdistances and the exploitation of redundancy betweengeodesic propagations, we introduce several methodsof computation of all pairs of distances and wecompare them.
PREREQUISITES
An image f is a discrete function defined on thefinite domain E ⊂ N , with N the set of positiveintegers. The values of a gray level image belongs to T ⊂ R . For a color image (i.e. with 3 channels) thevalues are in T = T × T × T and for a multivariateimage of L channels the values are in T L . In whatfollows, we consider T ⊂ R + . The whole resultspresented in the current paper are directly extendableto color or multivariate images (Noyel et. al. 2007a,2007b), (Noyel, 2008).An image is represented on a grid on which theneighborhood relations can be defined. Therefore, animage is seen as a non oriented graph G = { V G , E G } in which the vertices V G correspond to the coordinatesof pixels, V G ∈ Z × Z , and the edges, E G ∈ Z , givethe neighborhood relations between the pixels. Then,the notion of neighborhood of a pixel p in the gridis introduced as the set of pixels which are directlyconnected to it: ∀ p , q ∈ V G , p and q are neighbors ⇔ ( p , q ) ∈ E G , (1)where the ordered pair ( p , q ) is the edge which joinsthe points p and q . We assume that a pixel is notits own neighbor and the neighboring relations aresymmetrical. The neighborhood of pixel p , N G ( p ) ,defined a subset of V G of any size, such as: ∀ p , q ∈ V G , N G ( p ) = { q ∈ V G , ( p , q ) ∈ E G } . (2) Usually in image processing, the followingneighborhoods are defined: 4-neighborhood, 8-neighborhood or 6-neighborhood. In the sequel, weuse the 8-neighborhood. For our study, the choiceof the neighborhood has no influence since wecompare several methods using for each one the sameneighborhood.A path between two points x and y is a chain ofpoints ( x , x , . . . , x i , . . . , x l ) ∈ E such as x = x and x l = y , and for all i , ( x i , x i + ) are neighbours. Therefore, apath ( x , x , . . . , x l ) can be seen as a subgraph in whichthe nodes corresponds to the points and the edges arethe connections between neighbouring points.The geodesic distance d geo ( x , x l ) , or geodesictime, between two points of a graph, x and x l , isdefined as the minimum distance, or time, betweenthese two points, inf P { t P ( x , x l ) } . The geodesic path P geo is one of the paths linking these two points withthe minimum distance, P geo ( x , x l ) = ( x , . . . , x l ) : P geo ( x , x l ) = ( x , . . . , x l ) (3)such as d geo ( x , x l ) = inf P { t P ( x , x l ) } If the edges of the graph are weighted by thedistance between the nodes, t ( x i , x i (cid:48) ) , the geodesic pathis one of the sequences of nodes with the minimumweight.To generate a geodesic distance, several measuresof dissimilarities can be considered between twoneighbors pixels of position x i and x i (cid:48) and of positivegrey values f ( x i ) and f ( x i (cid:48) ) : – Pseudo-metric L d L ( x i , x i (cid:48) ) = | f ( x i ) − f ( x i (cid:48) ) | (4) – Pseudo-metric sum of grey levels: d + ( x i , x i (cid:48) ) = f ( x i ) + f ( x i (cid:48) ) (5) – Pseudo-metric mean of grey levels (similar to theprevious pseudo-metric): d + ( x i , x i (cid:48) ) = f ( x i ) + f ( x i (cid:48) ) P =( x , . . . , x l ) are the sum of the pseudo-metric along thispath P : t P ( x , x l ) = l ∑ i = d ( x i − , x i ) (7)2 mage Anal Stereol ?? (Please use \ volume ):1-10 The geodesic distance is defined as the distance alongthe geodesic path. It is also the minimum of distancesover all paths connecting two points x and x l : d geo ( x , x l ) = l ∑ i = { d ( x i − , x i ) | x i − , x i ∈ P geo } = inf P { t P ( x , x l ) } (8)In this paper, we only use the pseudo-metricsum of grey levels d + which presents the advantage(compared to d L ) not to be null when two pixels havethe same strictly positive value (if f ( x i ) = f ( x i (cid:48) ) > d + ( x i , x i (cid:48) ) = f ( x i ) > d + geo ( x , x l ) = l ∑ i = d + ( x i − , x i )= l ∑ i = f ( x i − ) + f ( x i ) (9) = f ( x ) + f ( x l ) + l − ∑ i = f ( x i ) In order to reduce the number of geodesicpropagations, Bertelli et al. (2006) have used thefollowing properties.
Property 1
Given a geodesic path ( p , p , . . . , p n ) , thegeodesic distance, d geo ( p i , p j ) , between two points p i et p j along the path, i < j, is equal to the differenced geo ( p , p j ) − d geo ( p , p i ) . In figure 1, if C and D belongs to a geodesic pathconnecting A and B the geodesic distance between C and D is equal to: d geo ( C , D ) = d geo ( A , D ) − d geo ( A , C ) = −
12 (10)with d geo ( C , D ) = d geo ( A , D ) =
20 and d geo ( A , C ) =
12. Fig. 1.
Discrete geodesic path on a 8-neighborhoodgraph. The points C and D belong to a geodesicpath between A et B. Therefore the geodesic distancebetween C and D is also known.
Using the property 1, when the geodesic distancebetween two points of the image is computed, thedistances between all pairs of points along theassociated geodesic path are known. Consequently, thedistances array is filled faster using this redundancy.In order to compute the distances between thepoints along the geodesic path, Bertelli et al. (2006)proposed to build a geodesic tree which has three kindsof nodes:1. the root, which is the source point for a geodesicpropagation. Its distance is null and it has noparents ;2. the nodes, which are points having both parentsand children ;3. the leaves, which are points without children.Starting from the leaves to the nodes (or theopposite), the distances between points belonging tothe same geodesic path are easily computed.The following additional property is very useful tocompute all pairs of geodesic distances.
Property 2
The longer the geodesic paths are, thehigher numbers of pairs of geodesic distances arecomputed.
In order to have the benefit of the property 2,Bertelli et al. (2006) have chosen as sources, of thegeodesic propagations, random points in a spiral-likeorder: starting from the edges of the image and goingto the center. Indeed, the points on the edges ofthe image tends to have longer geodesic paths thanpoints located at the center. Consequently, we want totest their remark by comparing their method to someothers.In order to make this comparison, we measure thefilling rates of the distances array D . For an imagecontaining N pixels, the distances array is a square3 OYEL G et al. : Fast computation of all pairs of geodesic distances matrix of size N × N = N elements. By symmetry, thenumber of geodesic paths to compute is equal to: A = N − N a is determined bycounting the unfilled elements of the distances array D .In practice, it is useful to use a boolean matrix D mrk , ofsize N × N , with elements equal to 1 if the distancebetween the pixel located by the line number andthe pixel located by the column number is computed,and 0 otherwise. By convention in the algorithm, weimpose to each element of the diagonal of the booleanmatrix to be equal to one, ∀ i D mrk ( i , i ) =
1, because thedistance from one pixel to itself is equal to zero. Dueto the symmetric properties of array D , the number ofcomputed paths is equal to: a = (cid:34)(cid:32) N ∑ k = N ∑ l = D mrk ( k , l ) (cid:33) − trace ( D mrk ) (cid:35) = (cid:34)(cid:32) N ∑ k = N ∑ l = D mrk ( k , l ) (cid:33) − N (cid:35) (12)Consequently, the filling rate of the distances arrayis defined by: τ = aA (13)The proportion of paths to compute for a given pixel x i ,is named the filling rate of the point x i , and is definedby: τ ( x i ) = ∑ Nk = D mrk ( k , i ) − D mrk ( i , i ) N − = ∑ Nk = D mrk ( k , i ) − N − τ ( x i ) = INTRODUCTION OF NEW METHODSFOR FAST COMPUTATION OF ALLPAIRS OF GEODESIC DISTANCES
In the current section, we initially present Bertelli et al. (2006) method, named the “spiral method”, andthen we introduce two new methods before makingcomparisons : 1) a geodesic extrema method and 2)a method based on the filling rate of all distancespairs array. The empirical comparisons are made onthree different images of size 25 ×
25 pixels: “bumps”,“hairpin bend”, “random” (fig. 2).
Image “bumps” Image “hairpin bend” Image “random”
Fig. 2.
Images of size × pixels “bumps”,“hairpin bend” and “random” whose grey levels arebetween 1 and 255. In order to use homogeneous measures for allmethods, the Soille’s algorithm (Soille, 2013), called“geodesic time function” is used, with a discreteneighborhood of size 3 × et al. ,2004). SPIRAL METHOD
Bertelli et al. (2006) affirms that the source pointsof the geodesic propagations, with the longest paths,are in general situated on the borders of the image.As these points are useful to reduce the number ofoperations, the source points are chosen in a randomway on concentric spiral turns of image pixels. Aconcentric spiral turn is a frame, of one pixel width,in which the top left corner is at position (1,1) or (2,2)or (3,3) or etc. The figure 3 gives an example. Whilenot all the pixels of the spiral turns have been selected,we draw one pixel, among them, in a uniform randomway ; otherwise we switch to the next spiral turn.Fig. 3.
The spiral turns of an image × pixels. mage Anal Stereol ?? (Please use \ volume ):1-10 Table 1.
Algorithm: Spiral method Given D the distances array of size N × N while D is not filled do Select the most exterior spiral turn not yetfilled Determine S the list of pixels of the spiral turnnot yet filled while S is not empty do Select a source point s randomly in S Compute the geodesic tree from s Fill the distances array D Remove the points of S which are filled end while end while For each image, the filling rate is plotted versusthe number of propagations (fig. 4). The number ofpropagations which are necessary to fill the distancesarray by the spiral method and the naive method arealso given in this figure. The relative difference of thenumber of propagations of the spiral method comparedto the number of propagations of the naive method iswritten ∆ r ( naive ) . We notice that the spiral methodreduces the number of propagations by a factor rangingbetween 13.4 % and 25.6 %, as compared to the naiveone. Consequently, it is very useful to exploit theredundancy in the propagations by building a geodesictree. (a) “bumps” (b) “hairpin bend” ∆ r ( naive ) = . ∆ r ( naive ) = . ∆ r ( naive ) = . Fig. 4.
Comparison of the filling rates τ of thedistances array, between the naive method (in green) and the spiral method (in blue) for the images“bumps”, “hairpin bend” and “random”. The relativedifferences ∆ r ( naive ) of the number of propagationsof the spiral method compared to the number ofpropagations of the naive method are given on thebottom line. SPIRAL METHOD WITH REPULSION
Two neighbours points have a high probabilityto have similar geodesic trees. Consequently, a firstimprovement of the spiral method is to introduced arepulsion distance between the points drawn randomlyin a spiral like-order (algorithm of table 2). Severaltests have shown us that a repulsion distance of threepixels on both sides of a source point gives the bestfilling rates.These tests are empirical tests. In fact, severalrepulsion distances were tried and it has been noticedthat the value of three pixels gives the best results inorder to fill the array of distances. This value of threepixels is certainly related to the image size, because,generally, the farther the source points of the geodesicpropagations are the faster the array of all pairs ofgeodesic distances is filled.Table 2.
Algorithm: Spiral method with repulsion Given D the distances array of size N × N Given h the repulsion distance: h ← while D is not filled do Select the most exterior spiral turn not yetfilled Determine S the list of pixels of the spiral turnnot yet filled while S is not empty do Select a source point s randomly in S Remove in the list S the two left points of s and the two right points of s if they are still in S Compute the geodesic tree from s Fill the array D Remove the points of S which are filled end while end while The figure 5, shows that the relative differencesin the number of propagations necessary to fill thedistances array are larger than 15 % for images“bumps” and “hairpin bend” and than 5.3 % for theimage “random”. Therefore, the spiral method with therepulsion distance is faster than the spiral method to fillthe distances array.5
OYEL G et al. : Fast computation of all pairs of geodesic distances (a) “bumps” (b) “hairpin bend” ∆ r ( spiral ) = . ∆ r ( spiral ) = . ∆ r ( naive ) = . ∆ r ( naive ) = . ∆ r ( spiral ) = . ∆ r ( naive ) = . Fig. 5.
Comparison of the filling rates τ of thedistances array, between the spiral method withrepulsion (in red) and the spiral method (in blue) forthe images “bumps”, “hairpin bend” and “random”.The relative differences ∆ r ( spiral ) (resp. ∆ r ( naive ) ) ofthe number of propagations of the spiral method withrepulsion compared to the number of propagations ofthe spiral (resp. naive) method are given on the bottomlines. GEODESIC EXTREMA METHOD
As the longest geodesic paths are those whichfill the most the distance array, we look for thegeodesic extrema of the image. To compute thegeodesic extrema of the image, we use two geodesicpropagations: – a first propagation starts from the edges of theimage. Then we select one of the points c withthe longest distance from the edges. This point iscalled the geodesic centroid ; – a second propagation from the geodesic centroidgives the farthest points from the centroid, i.e. thegeodesic extrema of the image.On the figure 6, which shows the results from thesetwo propagations, we notice that the geodesic extremaare mainly on the edges of the image, which ascertainsthe motivation to use a spiral method.Then the pixels are sorted by geodesic distancefrom 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 samegeodesic distance from the centroid, then the less filledis selected (algorithm of table 3). “bumps” “hairpin bend” “random”Geodesic distance from the edges of the imageGeodesic 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 D the distances array of size N × N Compute the geodesic centroid c Compute the decreasing list of geodesic extrema
Ext whose first element
Ext [ ] is the greatestgeodesic extrema not yet filled Initialise to zero, the list, of size N , of the fillingrate of points τ . while D is not filled do B ← { x ∈ Ext | d geo ( x , c ) = d geo ( Ext [ ] , c ) } s ← x Compute the geodesic tree from s Fill the distances array D Remove the points of the list
Ext which arefilled
Update the list of filling rates τ end while mage Anal Stereol ?? (Please use \ volume ):1-10 (a) “bumps” (b) “hairpin bend” ∆ r ( spiral ) = . ∆ r ( spiral ) = . ∆ r ( naive ) = . ∆ r ( naive ) = . ∆ r ( spiral ) = . ∆ r ( naive ) = . Fig. 7.
Comparison of the filling rates τ of thedistances array, between the geodesic extrema method(in red) and the spiral method (in blue) for the images“bumps”, “hairpin bend” and “random”. The relativedifferences ∆ r ( spiral ) (resp. ∆ r ( naive ) ) of the numberof propagations of the the geodesic extrema methodcompared to the number of propagations of the spiral(resp. naive) method are on the bottom lines. The filling rates of the spiral method and thegeodesic extrema method versus the number ofpropagations are plotted for each image (fig. 7). Wenotice that the filling rates of the geodesic extremamethod are at the beginning inferior or similar to theseof the spiral method. However, at the end the fillingrates of the geodesic extrema method are better thanthese of the spiral method. In fact, we are lookingfor an approach filling totally the distances array inthe fastest way. As the number of propagations of thegeodesic extrema method necessary to fill the distancesarray are lower than in the spiral method, the geodesicextrema method fills faster the distances array than thespiral method.In order to get an exact comparison it is necessaryto generate two propagations, corresponding to thedetermination of the geodesic extrema, to the numberof propagations necessary to fill the distances array.Even, with this modification, the extrema method fillsthe distances array faster than the spiral method (therelative differences ∆ r ( spiral ) are between 3.4 % and6 %). METHOD BASED ON THE FILLINGRATE OF THE DISTANCES ARRAY
In place of selecting the source points from theirdistance from the geodesic centroid propagation, weselect first the less filled points. To this aim, after eachgeodesic propagation the filling rate is computed foreach point. Then the less filled point is selected as asource of the propagation. If several points are amongthe less filled, then the greatest geodesic extrema ischosen among these points (algorithm of table 4).Table 4.
Algorithm: Method based on the filling rate ofthe distances array Given D the distances array of size N × N Compute the decreasing list of geodesic extrema
Ext Initialise to zero, the list, of size N , of the fillingrate of points τ . while D is not filled do B ← { argmin x ∈ E τ [ x ] } if Card { B } > then s ← argmin x ∈ B Ext [ x ] else s ← argmin x ∈ E τ [ x ] end if Compute the geodesic tree from s Fill the distances array D Remove the points of the list
Ext which arefilled
Update the list of filling rates τ end while As for the geodesic extrema method, we comparethe filling rates of the method based on the filling rateof the distances array and the spiral method versusthe number of propagations (fig. 8). We notice thatthe method based on the filling rate of D reducesof 32.3 % the number of propagations of the spiralmethod on the image “bumps” and 18.7 % on theimage “hairpin bend”. Consequently this method fillsfaster the distances array than the spiral method. Evenfor the image “random” the method based on the fillingrate of D still improves the spiral method of 2.7 %.However it is not very common to compute all pairsof geodesic distances on a strong unstructured imagesuch as the “random” one.By comparison to the naive approach, the methodbased on the filling rate of D reduces the number ofpropagations by a 49.6 % rate (resp. 29.6 %) on theimage “bumps” (resp. “hairpin bend”).7 OYEL G et al. : Fast computation of all pairs of geodesic distances
As for the previous method, in order to get an exactcomparison, it is necessary to add two propagations,corresponding to the determination of the geodesicextrema, to the number of propagations necessary tofill the distances array. Even, with this modification,the number of propagations necessary to fill thedistances array is still lower for the method based onthe filling rate of D . (a) “bumps” (b) “hairpin bend” ∆ r ( spiral ) = . ∆ r ( spiral ) = . ∆ r ( naive ) = . ∆ r ( naive ) = . ∆ r ( spiral ) = . ∆ r ( naive ) = . Fig. 8.
Comparison of the filling rates τ of thedistances array, between the method based onthe filling rate of the distances array (in red)and the spiral method (in blue) for the images“bumps”, “hairpin bend” and “random”. The relativedifferences ∆ r ( spiral ) (resp. ∆ r ( naive ) ) of the numberof propagations of the method based on the filling ratecompared to the number of propagations of the spiral(resp. naive) method are given on the bottom lines. DISCUSSION
After having presented and tested several methodsto fill the distances array, we have compared themfor the three test images “bumps”, “hairpin bend”and “random” on the figure 9 and in the table 5. Forthe images “bumps” and “hairpin bend”, the methodbased on the filling rate of the distances array isfaster than the other algorithms. For the “randomimage” (an extreme case presenting no texture) themethods introduced here give similar results, sincethe relative difference between the maximum numberof propagations is less than 2.7 %. Therefore, we conclude that the method based on the filling rate of thedistances array is the best one to calculate the array ofall pairs of distances. According to their performances,the others are ranked in the following order: 1) thespiral method with repulsion distance, 2) the geodesicextrema method and 3) the spiral method of Bertelli et al. (2006).Moreover, we have shown that the method basedon the filling rate of D reduces the number ofoperations between 19 % and 32 %, as compared tothe spiral approach and between 30 % and 50 %, ascompared 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 severalpoints have the minimum filling rate, seems to be thebest criterion to fill efficiently the distances array. (a) “bumps” (b) “hairpin bend”(c) “random” Fig. 9.
Comparison of the filling rates, of the array ofdistances, for the methods: spiral, geodesic extrema,the method based on the filling and the spiral methodwith a repulsion distance of 3 pixels for the images“bumps”, “hairpin bend” and “random”. mage Anal Stereol ?? (Please use \ volume ):1-10 Spiral spiral geodesic filling ∆ r ( naive ) method method with extrema raterepulsion method method“Bumps” 25.6 % 38.6 % 30.1 % “Hairpin bend” 13.4 % 26.4 % 17.9 % “Random” 16.2 % (a) spiral geodesic filling ∆ r ( spiral ) method with extrema raterepulsion method method“Bumps” 17.4 % 6.0 % “Hairpin bend” 15.0 % 5.2 % “Random” (b)Table 5. Relative difference values of filling rates (a) between the different methods and the naive approach,or (b) between the different methods and the spiralapproach. For each image is given in bold the bestrelative difference value.
CONCLUSION
From a comparison between different approaches,it turns out that a method based on the optimizationof the filling rate of the distances array is the mostefficient to compute the geodesic distances between allpairs of pixels in an image.Besides, in the current paper, we have shown ourresults on grey images. They can be directly extendedto hyperspectral images using appropriate pseudo-metrics (Noyel et al. , 2007a). This can be a useful stepfor a subsequent clustering by kernel methods or datareduction approaches on multivariate images.The main motivation for our developments is thecomputation of all pairs of geodesic distances forthe pixels of an image, which is usually a graph ofthousands of vertices arranged spatially. Nevertheless,our approach is valid on more general graphs, thanthose associated to bitmap images, after determiningthe “boundary vertices” of the graph, since then,the computation of the geodesic centre (and thegeodesic extremities) of a graph can be obtained bya first propagation from the “boundary vertices”. The“boundary vertices” can be defined, for instance, asthe vertices having less neighbouring vertices than theaverage number of connectivity in the graph.
REFERENCES
Bertelli L, Manjunath BS (2007). Edge preservingfilters using geodesic distances on weightedorthogonal domains. In: 2007 IEEE InternationalConference on Image Processing, vol. 1.Bertelli L, Sumengen B, Manjunath BS (2006).Redundancy in all pairs fast marching method.In: 2006 IEEE International Conference on ImageProcessing.Coeurjolly D, Miguet S, Tougne L (2004). 2d and3d visibility in discrete geometry: an applicationto discrete geodesic paths. Pattern RecognitionLetters 25:561 – 570. Discrete Geometry forComputer Imagery (DGCI’2002).Grazzini J, Soille P (2009). Edge-preservingsmoothing using a similarity measure in adaptivegeodesic neighbourhoods. Pattern Recognition42:2306 – 2316. Selected papers from the14th IAPR International Conference on DiscreteGeometry for Computer Imagery 2008.Lantuejoul C, Maisonneuve F (1984). Geodesicmethods in quantitative image analysis. PatternRecognition 17:177 – 187.Lerallut R, Decencière E, Meyer F (2007). Imagefiltering using morphological amoebas. Image andVision Computing 25:395 – 404. InternationalSymposium on Mathematical Morphology 2005.Mohan A, Sapiro G, Bosch E (2007). Spatiallycoherent nonlinear dimensionality reduction andsegmentation of hyperspectral images. IEEEGeoscience and Remote Sensing Letters 4:206–10.Noyel G (2008). Filtering, dimensionality reduction,classification and morphological segmentation ofhyperspectral images. Theses, École NationaleSupérieure des Mines de Paris.Noyel G, Angulo J, Jeulin D (2007a). Morphologicalsegmentation of hyperspectral images. ImageAnalysis Stereology 26:101–9.Noyel G, Angulo J, Jeulin D (2007b). Ondistances, paths and connections for hyperspectralimage segmentation. In: Banon G, et al. ,eds., International Symposium on MathematicalMorphology, vol. 1. Rio de Janeiro, Brazil:Instituto Nacional de Pesquisas Espaciais (INPE).ISBN 978-85-17-00035-5.Sethian J (1999). Level Set Methods and FastMarching Methods: Evolving Interfaces inComputational Geometry, Fluid Mechanics,Computer Vision, and Materials Science.Cambridge Monographs on Applied andComputational Mathematics. CambridgeUniversity Press.9