John Chaussard
email: j.chaussard@esiee.fr

Michel Couprie
email: m.couprie@esiee.fr

Hugues Talbot
email: h.talbot@esiee.fr

Robust skeletonization using the discrete lambda-medial axis

Medial axes and skeletons are notoriously sensitive to contour irregularities. This lack of stability is a serious problem for applications in e.g. shape analysis and recognition. In 2005, Chazal and Lieutier introduced the λ-medial axis as a new concept for computing the medial axis of a shape subject to single parameter filtering. The λ-medial axis is stable under small shape perturbations, as proved by these authors. In this article, a discrete λ-medial axis (DLMA) is introduced and compared with the recently introduced integer medial axis (GIMA). We show that DLMA provides measurably better results than GIMA, with regard to stability and sensibility to rotations. We give efficient algorithms to compute the DLMA, and we also introduce a variant of the DLMA which may be computed in linear-time.

Introduction

The notion of medial axis has been introduced by Blum in the 60s [START_REF] Blum | An associative machine for dealing with the visual field and some of its biological implications[END_REF][START_REF] Blum | A transformation for extracting new descriptors of shape[END_REF]. It has proved its usefulness in many practical applications, and numerous works were devoted to its properties and implementations. The original definition of the medial axis by Blum was based on a fire propagation analogy. However, its simplest definitions only require elementary geometry. In the continuous Euclidean space, the two following definitions can be used to formalise this notion: let X be a bounded subset of R n ; -Interpretation (a) of the medial axis of X consists of the centers of the ndimensional balls that are included in X but that are not included in any other n-dimensional ball included in X.

-Interpretation (b) of the medial axis of X consists of the points x ∈ X that have more than one nearest points on the boundary of X.

These two definitions differ only by a negligible set of points (see [START_REF] Matheron | Examples of topological properties of skeletons[END_REF]), in general interpretation (a) of the medial axis is a strict subset of interpretation (b). Notice that in some works, the term "skeleton" is used to refer to both interpretations, especially in the continuous framework. In this paper, we shall restrict the use of this term to the cases where the skeleton is topologically equivalent to the original shape.

To compute the medial axis approximately or exactly, different methods have been proposed, relying on different frameworks: discrete geometry [START_REF] Borgefors | The euclidean distance transform: finding the local maxima and reconstructing the shape[END_REF][START_REF] Ge | On the generation of skeletons from discrete euclidean distance maps[END_REF][START_REF] Malandain | Euclidean skeletons[END_REF][START_REF] Rémy | Exact medial axis with euclidean distance[END_REF][START_REF] Hesselink | Euclidean skeletons of digital image and volume data in linear time by the integer medial axis transform[END_REF], digital topology [START_REF] Davies | Thinning algorithms: a critique and a new methodology[END_REF][START_REF] Vincent | Efficient computation of various types of skeletons[END_REF][START_REF] Talbot | Euclidean skeletons and conditional bisectors[END_REF][START_REF] Pudney | Distance-ordered homotopic thinning: a skeletonization algorithm for 3D digital images[END_REF], mathematical morphology [START_REF] Serra | Image analysis and mathematical morphology[END_REF][START_REF] Soille | Morphological image analysis[END_REF], computational geometry [START_REF] Attali | Delaunay conforming iso-surface, skeleton extraction and noise removal[END_REF][START_REF] Ogniewicz | Hierarchic voronoi skeletons[END_REF][START_REF] Attali | Modelling noise for a better simplification of skeletons[END_REF], partial differential equations [START_REF] Siddiqi | The hamiltonjacobi skeleton[END_REF], and levelsets [START_REF] Kimmel | Skeletonization via distance maps and level sets[END_REF]. In this paper, we focus on medial axes in the discrete grid Z 2 or Z 3 , which are centered in the shape with respect to the Euclidean distance.

A major difficulty when using the medial axis in applications (e.g., shape recognition), is its sensitivity to small contour perturbations, in other words, its lack of stability. A recent survey [START_REF] Attali | Stability and computation of the medial axis -a state-of-the-art report[END_REF] summarises selected relevant studies dealing with this topic. This difficulty can be expressed mathematically: the transformation which associates a shape to its medial axis is only semi-continuous. This fact, among others, explains why it is usually necessary to add a filtering step (or pruning step) to any method that aims at computing the medial axis.

Hence, there is a rich literature devoted to medial axis pruning, in which different criteria were proposed in order to discard "spurious" medial axis points or branches: see Attali et al. [START_REF] Attali | Pruning discrete and semicontinuous skeletons[END_REF], Ogniewicz and Kübler [START_REF] Ogniewicz | Hierarchic voronoi skeletons[END_REF], Attali and Montanvert [START_REF] Attali | Modelling noise for a better simplification of skeletons[END_REF], Malandain and Fernández-Vidal [START_REF] Malandain | Euclidean skeletons[END_REF], Attali and Lachaud [START_REF] Attali | Delaunay conforming iso-surface, skeleton extraction and noise removal[END_REF], Svensson and Sanniti di Baja [START_REF] Svensson | Simplifying curve skeletons in volume images[END_REF], Bai et al. [START_REF] Bai | Skeleton pruning by contour partitioning with discrete curve evolution[END_REF], Couprie et al. [START_REF] Couprie | Discrete bisector function and Euclidean skeleton in 2D and 3D[END_REF], Hesselink and Roerdink [START_REF] Hesselink | Euclidean skeletons of digital image and volume data in linear time by the integer medial axis transform[END_REF], to cite only a few. However, we lack theoretical justification, that is, a formalised argument that would help to understand why a filtering criterion is better than another.

In 2005, Chazal and Lieutier introduced the λ-medial axis and studied its properties, in particular those related to stability [START_REF] Chazal | The λ-medial axis[END_REF]. Consider a bounded subset X of R n , as for example, for n = 2, the region enclosed by the solid curve depicted in Fig. 1 (left). Let x be a point in X, we denote by Π(x) the set of points of the boundary of X that are closest to x. For example in Fig. 1, we have Π(x) = {a, b}, Π(x ′ ) = {a ′ , b ′ } and Π(x ′′ ) = {a ′′ }. Let λ be a non-negative real number, the λ-medial axis of X is the set of points x of X such that the smallest n-dimensional ball1 including Π(x) has a radius greater than or equal to λ. Notice that the 0-medial axis of X is equal to X, and that any λ-medial axis with λ > 0 is included in the medial axis according to definition b). We show in Fig. 1 (right) two λ-medial axes with different values of λ.

A major outcome of Chazal and Lieutier [START_REF] Chazal | The λ-medial axis[END_REF] is the following property: informally, for "regular" values of λ, the λ-medial axis remains stable under perturbations of X that are small with respect to the Hausdorff distance. Typical non-regular values are radii of locally largest maximal n-dimensional balls.

This property is a strong argument in favor of the λ-medial axis, especially in the absence of such result for other proposed criteria.

In the field of computational geometry, the λ-medial axis has been exploited in particular by Samozino et al. [START_REF] Samozino | Reconstruction with voronoi centered radial basis functions[END_REF] to propose a robust method for reconstructing surfaces from point clouds. Also, notions closely related to the λ-medial axis have led Chazal et al. to propose stable approximations of tangent planes and normal cones from noisy samples [START_REF] Chazal | Normal cone approximation and offset shape isotopy[END_REF]. In the discrete grids, namely Z 2 and Z 3 , a similar filtering criterion has been considered in independent works [START_REF] Malandain | Euclidean skeletons[END_REF][START_REF] Hesselink | Euclidean skeletons of digital image and volume data in linear time by the integer medial axis transform[END_REF]. It consists of selecting for each medial axis point two of its closest boundary points, and using the distance between these two points as filtering criterion. The work of Hesselink and Roerdink [START_REF] Hesselink | Euclidean skeletons of digital image and volume data in linear time by the integer medial axis transform[END_REF] provides a linear-time algorithm to compute a filtered medial axis based on this criterion, which exhibits good noise robustness properties in practice.

In this article (which extends Chaussard et al. [START_REF] Chaussard | A discrete lambda-medial axis[END_REF], a preliminary version published in the DGCI conference proceedings), we introduce the definition of a discrete λ-medial axis (DLMA) in Z n . We evaluate experimentally its stability and rotation invariance in 2D and 3D. In this experimental study, we compare it with the previously introduced integer medial axis (GIMA) [START_REF] Hesselink | Euclidean skeletons of 3d data sets in linear time by the integer medial axis transform[END_REF][START_REF] Hesselink | Euclidean skeletons of digital image and volume data in linear time by the integer medial axis transform[END_REF] and show that the DLMA provides measurably better results. Furthermore, we introduce a variant of the DLMA which may be computed in linear time, for which the results are very close to those of the DLMA, and which is only slightly slower than the one proposed by Hesselink and Roerdink [START_REF] Hesselink | Euclidean skeletons of digital image and volume data in linear time by the integer medial axis transform[END_REF].

The λ-medial axis

Let us first recall the original definition of the λ-medial axis given by Chazal and Lieutier.

Let x = (x 1 , . . . , x n ), y = (y 1 , . . . , y n ) ∈ R n , the Euclidean distance between x and y is denoted by d(x, y), in other terms, d(x, y) = (

n k=1 (y k -x k ) 2 ) 1 2 .
Let X be a finite subset of R n or Z n . We set d(y, X) = min x∈X {d(y, x)}. We denote by |X| the number of elements of X, and by X the complement of X.

Let E be either R n or Z n . Let x ∈ E, r ∈ R + , we denote by B r (x) the n-dimensional ball (or ball for simplicity) of radius r centered on x, defined by B r (x) = {y ∈ E | d(x, y) ≤ r}. We also define Br (x) = {y ∈ E | d(x, y) < r}.

Let X ⊆ E. A ball B r (x) ⊆ X is maximal for X if it is not strictly included in any other ball included in X. The Euclidean medial axis of X, denoted by EMA(X), is the set of the centers of all the maximal balls for X.

Let S be a non-empty subset of E, and let x ∈ E. The projection of x on S, denoted by Π S (x), is the set of points y of S which are at minimal distance from x ; more precisely,

Π S (x) = {y ∈ S | ∀z ∈ S, d(y, x) ≤ d(z, x)}. If X is a subset of E, the projection of X on S is defined by Π S (X) = x∈X Π S (x).
Let S ⊂ R n , we denote by R(S) the radius of the smallest ball enclosing S, that is,

R(S) = min{r ∈ R | ∃y ∈ R n , B r (y) ⊇ S}.
The λ-medial axis may now be defined based on these notions.

Definition 1 (Chazal and Lieutier [START_REF] Chazal | The λ-medial axis[END_REF]). Let X be an open bounded subset of R n , and let λ ∈ R + . The λ-medial axis of X is the set of points x in X such that R(Π X (x)) ≥ λ.

A discrete λ-medial axis

Transposing the original definition of the λ-medial axis directly to Z n would lead to an unsatisfactory result. For instance, consider a horizontal ribbon in Z 2 with constant, even width and infinite length. Clearly, the projection of any point of this set on its complementary set is reduced to a singleton. This is why, if we keep the same definition, any λ-medial axis of this object with λ > 0 would be empty.

In order to avoid such unwanted behaviour, we replace the projection by the so-called extended projection [START_REF] Couprie | Discrete bisector function and Euclidean skeleton in 2D and 3D[END_REF]. The extended projection was originally introduced in order to propose a discrete definition of the bisector function, another indicator used to filter skeletons.

For each point x ∈ Z n , we define the direct neighborhood of x as N (x) = {y ∈ Z n | d(x, y) ≤ 1}. The direct neighborhood comprises 2n + 1 points.

The set Π X (X) is composed of the points of X that are in the direct neighborhood of a point of X, that is, Π X (X) = {y ∈ X | N (y) ∩ X = ∅}. We call Π X (X) the (external) boundary of X.

Let X ⊆ Z n , and let x ∈ X. The extended projection of x on X, denoted by Π e X (x), is the union of the sets Π X (y), for all y in N (x) such that

d(y, X) ≤ d(x, X). ( condition 1)
Using the extended projection allows us to cope with the problem pointed out previously. Considering the example of an horizontal ribbon in Z 2 with constant even width and infinite length, points close to the "middle" of the ribbon will each have an extended projection consisting of a pair of points. Informally, the extended projection permits to "see" both sides of the object from a point x close to the "middle" of the object, whereas only one side is taken into account when considering the projection.

As illustrated in Fig. 2, condition 1 avoids producing multiple medial axis points when only one is sufficient, in other words, it yields a thinner axis according to definition 2.

Fig. 2 helps us to illustrate this definition and the following one. Consider the object X (white pixels) depicted in Fig. 2a, and the two pixels x, y in X. Fig. 2b shows the squared Euclidean distance map of X, that is, the value in each pixel is the square Euclidean distance between this pixel and the nearest background pixel. We can see that d(x, X) = 4 and d(y, X) = 5. The projection of x (resp. y) on X is {b} (resp. {d}). We have Π e X (x) = {a, b, c}, and Π e X (y) = {d, b, c}. We are now ready to introduce the definition of our discrete λ-medial axis.

Definition 2. Let X be a finite subset of Z n , and let λ ∈ R + . We define the function F X which associates, to each point x of X, the value

F X (x) = R(Π e X (x)). The discrete λ-medial axis (or DLMA) of X is the set of points x in X such that F X (x) ≥ λ.
Back to our illustration in Fig. 2, the function F X is displayed in Fig. 2c and the 4-medial axis of X (threshold of F X for λ = 4) is displayed in Fig. 2d. Notice that the function F X can be computed once and stored in an n-dimensional array (a grayscale image), and that any DLMA of X is obtained by thresholding the values of this array at a particular value λ.

In Fig. 2e, we show how the function F X would look like if we omit condition 1 from the definition of the extended projection. In this case, the extended projection of x would be {a, b, c, d} instead of {a, b, c}, this explains the higher value of F X (x). Fig. 2f shows how the 4-medial axis of X would be.

In Fig. 3(b,c), we show two examples of DLMAs of a bigger shape. More examples are given in Fig. 4 and Fig. 5. Notice that DLMA in general does not exhibit the same topological characteristics as the original shape.

The integer medial axis

The integer medial axis was introduced by [START_REF] Hesselink | Euclidean skeletons of digital image and volume data in linear time by the integer medial axis transform[END_REF]. Conceptually, this notion is close to the DLMA, and also defined in the framework of discrete grid Z n .

Let X ⊂ Z n , and x ∈ X. We define Π ′ X (x) as the element of Π X (x) that is the smallest with respect to the lexicographic ordering of its coordinates.

Definition 3 (Hesselink and Roerdink [START_REF] Hesselink | Euclidean skeletons of digital image and volume data in linear time by the integer medial axis transform[END_REF]). Let X be a finite subset of Z n , and let γ ∈ R + . The γ-integer medial axis (or GIMA) of X is the set of points x in X such that at least one y

∈ N (x) verifies : i) d(Π ′ X (x), Π ′ X (y)) > γ, and ii) d(m -Π ′ X (y)) ≤ d(m -Π ′ X (x))
, where m = x+y 2 . The integer medial axis of X is the γ-integer medial axis of X for γ = 1.

(a)

x y b a d c 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 1 4 2 1 1 1 0 0 0 0 0 1 4 5 4 4 2 1 1 1 0 0 1 2 4 4 5 5 4 4 2 1 0 0 1 1 1 2 4 4 5 4 1 0 0 0 0 0 1 1 1 2 4 1 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 Notice that the GIMA is not invariant under coordinate permutation.

There are indeed some links between the GIMA and the DLMA. In 2D, in the case of a point x that is, conceptually, a "regular medial axis point" (neither a branch extremity nor a branch junction), condition i) is similar to the condition R(Π X (x)) ≥ λ in Def. 1.

Condition ii) plays a role analogous to condition d(y, X) ≤ d(x, X) in the definition of the extended projection, that is to get a thinner axis.

However, these definitions differ, leading to sensible differences in the results of these two transformations (see Fig. 5). For example, in definition 3, Π ′ X (x) is the only element of the projection Π X (x) that is taken into account into the computation of the integer medial axis of X, while all the elements of Π e X (x) are taken into account into the computation of the lambda medial axis. Moreover, in definition 3, the result depends on the longest distance between two elements of the (reduced) projection, while in definition 2, the result depends on the size of the smallest ball including the whole projection.

Generally speaking, one can say that the lambda medial axis "takes into account" more elements than the integer medial axis. We analyse quantitatively the consequences of these differences in Sec. 7 and Sec. 8. These differences lead to different axes, as illustrated on Fig. 2.

Topology preservation issues

The illustration on Fig. 5 (right) is sufficient to demonstrate that a DLMA or a GIMA of a given shape X may have a homotopy type different from the one of X.

Nearly identical methods for filtering a medial axis while preserving homotopy were proposed by several authors [START_REF] Davies | Thinning algorithms: a critique and a new methodology[END_REF][START_REF] Vincent | Efficient computation of various types of skeletons[END_REF][START_REF] Talbot | Euclidean skeletons and conditional bisectors[END_REF][START_REF] Pudney | Distance-ordered homotopic thinning: a skeletonization algorithm for 3D digital images[END_REF]. They consist of iteratively removing simple points [START_REF] Kong | Digital topology: introduction and survey[END_REF][START_REF] Bertrand | Simple points, topological numbers and geodesic neighborhoods in cubic grids[END_REF][START_REF] Couprie | New characterizations of simple points in 2D, 3D and 4D discrete spaces[END_REF] from the shape X, that do not belong to the set M of its medial axis points (called constraint set). A priority function is needed in order to specify which points must be considered at each step of the thinning. In the general case, the choice of this priority function is not obvious [START_REF] Talbot | Euclidean skeletons and conditional bisectors[END_REF][START_REF] Couprie | Discrete bisector function and Euclidean skeleton in 2D and 3D[END_REF].

An example is shown on Fig. 6, where the filtered DLMA of an object is used as a constraint set during skeletonization. Choosing the exact Euclidean distance map as a priority function for removing simple points from the object may lead to geometric distortions [START_REF] Talbot | Euclidean skeletons and conditional bisectors[END_REF]. In some cases, "extra branches" may even appear (see Fig. 6b). Choosing the map F X as priority function yields more satisfying results (see Fig. 6c), as it guides the thinning process towards elements that belong to the different nested discrete λ-medial axes.

Implementation issues

In this section, we discuss the algorithms that should be used in order to efficiently implement DLMA, and we indicate their complexity. We also introduce a variant of the DLMA that can be computed in linear time.

The set of all projections Π X (x), for all x ∈ X, can be computed in optimal time and space (that is, in O(n) where n = x∈X |Π X (x)|), thanks to an algorithm introduced by Coeurjolly in [START_REF] Couprie | Discrete bisector function and Euclidean skeleton in 2D and 3D[END_REF]. After this, to compute the extended projection of a particular point x, we just have to perform the union of the projections of all points y in N (x) that satisfy condition 1.

In order to avoid computing the map Π X , which requires a data structure like an array of lists to store all the sets Π X (x), we propose a variant of the DLMA that consists of replacing Π X by Π ′ X in the definition of the extended projection. We refer to this variant as the discrete linear λ-medial axis (DLLMA): it is illustrated on Fig. 2g andh. Notice that the DLLMA, as the GIMA, is not invariant under a permutation of the coordinates (on the contrary of the DLMA). One advantage of this variant is that the map Π ′ X can be computed in linear time and stored in an array of integers [START_REF] Coeurjolly | Algorithmique et géométrie discrète pour la caractérisation des courbes et des surfaces[END_REF][START_REF] Maurer | A linear time algorithm for computing exact euclidean distance transforms of binary images in arbitrary dimensions[END_REF][START_REF] Hesselink | Euclidean skeletons of digital image and volume data in linear time by the integer medial axis transform[END_REF]. We included DLLMA in our experimental comparative evaluation (see Sec. 7 and Sec. 8). The smallest d-dimensional ball enclosing a given set of n points in R d can be computed in expected O(n) time if d is considered as a constant, thanks to an algorithm proposed by Welzl [START_REF] Welzl | Smallest enclosing disks (balls and ellipsoids)[END_REF]. This algorithm is simple and easy to implement.

For each object point x, the computational cost is in O(n) where n = |Π e X (x)| (see Def. 2). In the variant DLLMA, we have n ≤ 5, since only one projection point is considered for x and each of its four direct neighbors, at worse. Hence, computing the radius of the smallest enclosing d-dimensional ball can be performed in constant time.

Overall, the DLLMA procedure runs in linear time with regard to the number of object points.

Using local characterisations, in 2D and 3D, checking whether a point is simple or not can be done in constant time. The homotopic thinning described in section 4 may be implemented to run in O(n log n) complexity, where n is the number of object points, using e.g. a balanced binary tree to store and retrieve the candidate points and their priority.

Methodology for the comparison of medial axes

In order to compare the robustness of DLMA, DLLMA and GIMA, we measure the dissimilarity between the medial axis of a shape X and the one of a shape X ′ derived from X by a small perturbation. In Sec. 7 this perturbation is obtained through a random process altering the contour of X, and in Sec. 8 it is obtained via a discrete rotation. We also include in this comparison the filtered Euclidean Medial Axis (EMA) defined as the set of the centers of maximal balls that have a radius greater than a parameter ρ.

We first introduce some notions that will help us to ensure a fair comparison between methods, for which the filtering parameter does not have the same meaning.

Let X be a finite subset of Z n . Let Y be a subset of X, we set REDT X (Y ) = y∈Y Bd(y,X) (y). The transformation REDT X is sometimes called reverse Euclidean distance transform [START_REF] Coeurjolly | d-dimensional reverse Euclidean distance transformation and Euclidean medial axis extraction in optimal time[END_REF][START_REF] Coeurjolly | Optimal separable algorithms to compute the reverse euclidean distance transformation and discrete medial axis in arbitrary dimension[END_REF].

It is well known that any object can be fully reconstructed from its medial axis, more precisely, we have X = REDT X (M ) whenever M is the (exact and non-filtered) medial axis of X. This property holds if M is the set of the centers of all maximal balls of X (EMA of X with ρ = 0), the DLMA of X with λ = 1, or the integer medial axis of X with γ = 1. However, it is no longer true if we consider filtered medial axes, e.g. EMA, DLMA or GIMA with arbitrary ρ, λ or γ.

Then, it is interesting to measure how much information about the original object is lost when we raise the filtering parameter. We set

L X (λ) = |X \ REDT X (LM X (λ))| |X|
where LM X (λ) is the DLMA of X. In other words, L X (λ) is the area of the difference between X and the set reconstructed from its DLMA, divided by the area of X. We call L X (λ) the (normalised) residuals of the DLMA of X filtered at value λ. The same way, we define L ′ X (λ ′ ) (resp. I X (γ), M X (ρ)) as the (normalised) residuals of the DLLMA (resp. GIMA, EMA) of X filtered at value λ ′ (resp. γ, ρ).

Fig. 8 shows the evolution of the residuals for the EMA, GIMA, DLMA and DLLMA for a given shape. Since differences are not negligible, to ensure a fair evaluation we compare the results of methods for approximately equal values of their residuals, rather than for equal values of their parameters. For comparing shapes or medial axes, we use the Hausdorff distance (see below), and also a dissimilarity measure proposed by Dubuisson and Jain [START_REF] Dubuisson | A modified hausdorff distance for object matching[END_REF]. The drawback of Hausdorff distance for measuring shape dissimilarity is its extreme sensibility to outliers, the latter measure performing better in this respect.

Let X, Y be two subsets of R n . We set

H(X|Y ) = max

Stability

In this section, we compare the stability of DLMA, DLLMA, GIMA and EMA with respect to contour noise.

We based our 2D experiments on Kimia's database of 216 shapes [START_REF] Sharvit | Symmetry-based indexing of image databases[END_REF]. For the 3D case, we used 20 three-dimensional objects gathered on the Internet or that we created ourselves. Fig. 9 shows a sample of these shapes . Medial axes are notoriously sensitive to border noise. Since the λ-medial axis is supposed to cope reasonably well with shape deformation, it is useful to test how it fares in practice.

To introduce noise to the boundary of an object we propose deforming it using a process derived from the Eden's (accretion) process [START_REF] Eden | A two-dimensional growth process[END_REF]. The Eden's process is an iterative random cellular automaton which, in its simplest form, attributes an equal probability to all the outer border points to be set to 1 at each step. That is, at each step, a neighbour of the object is chosen randomly and added to the shape. In spite of its simplicity, the Eden's process exhibits good asymptotic isotropy [START_REF] Lee | A stochastic tesselation of digital space[END_REF]. In our process, we required accretion steps to deal only with simple points. This way at each step, the object's homotopy type remains unchanged.

We denote by E(X, n) the result of applying n steps of this process to the shape X. In this experiment, we compare the (filtered) medial axis of an original shape with the one of a deformed shape. 2: Average dissimilarity between M (X) and M (E(X, n)). Noise level (parameter n) is expressed as a percentage of object area, results shown are averages of all tested residuals. Lowest values are highlighted.

In Table 1 (resp. Table 2) we give the average Hausdorff distance (resp. dissimilarity) between M (X) and M (E(X, n)) on the 216 shapes of Kimia's database and the 20 shapes of our 3d database for different definitions of M (filtered EMA, GIMA, DLMA or DLLMA of X). Filtering parameters ρ, γ, λ and λ ′ were chosen in order to give a residual varying from 5% to 30% in 2d, and a residual varying from 3% to 10% in 3d. Results are given for values of n equal to 3%, 6% and 9% of object's surface/volume (values shown are averages of results obtained for all tested residuals), and the last column indicates all the results averaged together.

Note that the tested residual values are not the same in the 2d and 3d case. Indeed, in 3d, choosing a residual higher than 10% gives a medial axis which contains very little information about the original object (surfacic parts are mostly deleted). We constrained our study to medial axes holding enough information on the original objects, and ignored residual superior to 10% in the 3d case.

On Fig. 11, we show the evolution of the Hausdorff distance and the dissimilarity between M (X) and M (E(X, n)) against the quantity of noise added to the shape, for DLMA, DLLMA, GIMA and EMA algorithms. The results shown are average of the results obtained on the 216 shapes of Kimia's 2D image database. Fig. 12 shows the same experimental results on the 20 shapes of our 3d database.

As we expected, DLMA fares better on the average than other axes, with regard to Hausdorff distance criterion. Its faster variant DLLMA ranks second. According to the dissimilarity criterion, DLMA and DLLMA yield the best results in 3d, and in 2d, DLMA performs nearly as well as EMA.

Sensitivity to rotations

In the continuous framework, it is well known that the medial axis is a rotation-invariant shape descriptor. More precisely, if we denote by R θ the rotation of angle θ about the origin, and by M the medial axis transform, the rotation invariance property states that M (R θ (X)) = R θ (M (X)), whatever the shape X and the angle θ.

In a discrete framework, this property holds only in particular cases (e.g., when θ is a multiple of 90 degrees). Nevertheless, we can experimentally measure the dissimilarity between M (R θ (X)) and R θ (M (X)) for different instances, and different definitions of the medial axis. The lower this dissimilarity, the more stable under rotation the method is.

For each shape of Kimia's database and of our database, we computed the Hausdorff distance and the dissimilarity between M (R θ (X)) and R θ (M (X)) for values of the parameters yielding 5%, 10%, 15%, 20%, 25% and 30% residuals, and for rotation angles varying from 0 to 90 degrees by 3 degrees steps. In 3d, rotations were successively performed around the Y and Z axis. Moreover, in 3d, only parameters yielding 3%, 5% and 10% residuals were tested: higher residual values yielded skeletons that were not keeping enough information from the original object (see discussion in Sec. 7).

Fig. 13 shows detailed results of such experiment in 2d, while Fig. 14 and Fig. 15 show results in 3d. For both cases, an inverse truncated real rotation algorithm was used to perform images rotation. The results are summarised in Tables 3 and4. In Table 3 (resp. Table 4) we give the average Hausdorff distance (resp. dissimilarity) between M (R θ (X)) and R θ (M (X)) on the 216 shapes of Kimia's database and on the 20 3d shapes of our database, for residual values 5%, 10%, 15%, 20%, 25% and 30% (for 3d shapes, only 3%, 5% and 10% were tested) and for θ varying from 0 to 89 degrees (in 3d shapes, rotation was performed independently around Y and Z axes). We show, in the first three columns, the average value obtained for all tested rotation angles for a given residual, while the fourth column indicates the average value obtained for all residuals and all rotation angles.

In 3d, DLMA and DLLMA yield, in average, the best results for both measures. In 2d, they are better than other axes with regards to Hausdorff distance, and are only outperformed by EMA with regards to dissimilarity.

Computational cost

In addition to complexity analysis reported in Sec. 5, we made some experiments to compare actual computation costs of the presented methods.

In Fig. 16, we show the results of computing time measurements that we performed on an Intel Core 2 Duo processor at 1.83 GHz. Computing times for the GIMA are the lowest, but DLLMA is only slightly slower (and also linear).

Conclusion

We introduced the discrete λ-medial axis (DLMA), and its variant DLLMA, a transposition in discrete grids of the continuous λ-medial axis recently proposed and known for its stability property. The experimental study reported in this paper shows that both DLMA and DLLMA are robust w.r.t. boundary noise and rotations in practice. For our experimental study of rotational invariance and border noise robustness, both in 2D and 3D, we introduced an original methodology that ensures a fair comparison between different methods, under the mere assumption that their result decreases in size under the control of a single parameter. This methodology could be of interest for comparing other medial axis filtering approaches.

2D

3D Residual 10% 20% 30% all 3% 5% 10% all EMA 1,21 1,34 3,88 1,71 0,72 0,72 0,73 0,72 GIMA 1,50 2,41 6,09 2,27 1,40 1,43 1,52 1,45 DLMA 1,46 1,89 3,72 1,94 0,66 0,66 0,76 0,69 DLLMA 1,37 1,78 3,60 1,85 0,67 0,68 0,78 0,71

Figure 1 :

 1 Figure 1: Illustration of the λ-medial axis. Left: Points x, x ′ and x ′′ and their respective closest boundary points. Top right: λ-medial axis with λ = ǫ, a very small positive real number. Bottom right: λ-medial axis with λ = d(a ′ , b ′ )/2 + ǫ.

  0 0.0 0.0 0.0 0.0 0.0 0.0 2.2 4.0 2.8 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 (d) (e) 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.8 4.1 2.2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.6 4.5 4.1 1.0 3.2 2.2 0.0 0.0 0.0 0.0 0.0 0.0 3.2 5.2 5.4 5.1 5.3 2.2 1.0 3.2 2.2 0.0 0.0 0.0 2.2 4.2 5.3 5.1 5.2 5.2 5.1 5.3 4.2 2.2 0.0 0.0 0.0 2.2 3.2 1.0 2.2 5.3 5.1 5.4 5.2 3.2 0.0 0.0 0.0 0.0 0.0 0.0 2.2 3.2 1.0 4.1 4.5 3.6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.2 4.1 2.8 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 (f) (g) 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.8 2.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.2 4.5 1.4 1.0 2.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 2.0 5.1 5.4 5.1 3.2 1.4 1.0 2.0 1.0 0.0 0.0 0.0 1.0 1.4 3.2 5.1 5.2 5.2 5.1 3.2 1.4 2.2 0.0 0.0 0.0 2.2 3.2 1.0 0.0 3.2 5.1 5.4 3.2 3.2 0.0 0.0 0.0 0.0 0.0 0.0 2.2 3.2 1.0 0.0 4.5 2.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.2 4.0 2.8 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 (h) (i) 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.2 1.0 1.4 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 3.6 0.0 1.0 1.0 1.4 0.0 0.0 0.0 0.0 0.0 0.0 1.0 4.5 5.1 5.0 2.2 0.0 1.0 1.0 1.4 0.0 0.0 0.0 1.4 0.0 2.2 1.0 5.0 5.1 5.0 2.2 0.0 2.2 0.0 0.0 0.0 2.2 1.0 1.0 0.0 2.2 1.0 4.5 4.5 1.0 0.0 0.0 0.0 0.0 0.0 0.0 2.2 1.0 1.0 0.0 4.1 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.2 0.0 2.2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 (j)

Figure 2 :

 2 Figure 2: (a): An object X (white pixels). We have Π e X (x) = {a, b, c}, and Π e X (y) = {d, b, c}. (b): The squared Euclidean distance map of X. (c): The function F X (truncated to the first decimal). (d): The λ-medial axis of X, for λ = 4 (black disks). (e,f): Same as (c,d) assuming that condition 1 is omitted from the definition of the extended projection. (g): The map, for each x ∈ X, of the minimal value λ such that x ∈ λ ′ -medial axis of X. (h): The λ ′ -medial axis of X, for λ ′ = 4 (black disks). (i,j): Same as (g,h) for the integer medial axis.

Figure 3 :

 3 Figure 3: (a): The function F X superimposed to the shape X. Darkest colors represent highest values of F X (x). Any DLMA of X is a threshold of this function at a particular value λ. (b): Discrete 5-medial axis. (c): discrete 10-medial axis of X.

Figure 4 :

 4 Figure 4: Top-right: The function F X superimposed to the top-left shape X. Red color represents highest values of F X (x), and blue color represents lowest ones. Bottom left: discrete 11-medial axis. Bottom right: discrete 15-medial axis of X.

Figure 5 :

 5 Figure 5: Results of DLMA (first row) and GIMA (second row) for parameter values yielding similar reconstruction ratios (see Sec. 6).

Figure 6 :

 6 Figure 6: (a): The original object X (in white). (b): Result of the skeletonization using the distance map as a priority function, and a filtered DLMA as constraint set. (c): Result of the skeletonization using F X as a priority function, and a filtered DLMA as constraint set (same set than in b).

Figure 7 :

 7 Figure 7: Results of homotopic thinning of X constrained by the DLMA for parameter values 5, 10, 20 and 30, and with F X as priority function.

Figure 8 :

 8 Figure8: Residuals L X (λ) (2 variants) and I X (γ), for the set X depicted in Fig.3. Horizontal axis: the value of the parameter (λ or γ). Vertical axis: the value of the residual. Notice that curves corresponding to DLMA and DLLMA are superimposed.

  y)}}, and d H (X, Y ) = max{H(X|Y ), H(Y |X)} is the Hausdorff distance between X and Y . We set D(X|Y ) = 1 |X| x∈X min y∈Y {d(x, y)}, and d D (X, Y ) = max{D(X|Y ), D(Y |X)} is the Dubuisson and Jain's dissimilarity measure between X and Y (called dissimilarity in the sequel for the sake of brevity).

Figure 9 :

 9 Figure 9: A sample of the 216 shapes of Kimia's database, and of the 3d shapes from our database.

Figure 10 :

 10 Figure 10: Illustration of the border deformation process, after 0 (original image), 50, 100, 200 and 400 steps. The function F X is superimposed to each shape X.

Figure 11 :

 11 Figure 11: Hausdorff distance (on the left) and dissimilarity measure (on the right) between M (X) and M (E(X, n)) on the 216 shapes of Kimia's 2D image database for residual values 15% and 30%. Noise level (parameter n) is expressed as a percentage of object area.

Figure 12 :

 12 ABBAAADCBCDF9DDCBA

Figure 13 :

 13 Figure 13: Hausdorff distance (on the left) and dissimilarity measure (on the right) between M (R θ (X)) and R θ (M (X)) on the 216 shapes of Kimia's 2D image database for residual values 15% and 25%. Noise level (parameter n) is expressed as a percentage of object area.

4 Figure 14 :

 414 Figure 14: Hausdorff distance (on the left) and dissimilarity measure (on the right) between M (R θ (X)) and R θ (M (X)) on the 20 shapes of our 3D image database for residual values 3%, 5% and 10% (the horizontal axis is proportional to the parameter n).

Figure 15 :

 15 Figure 15: Hausdorff distance (on the left) and dissimilarity measure (on the right) between M (R θ (X)) and R θ (M (X)) on the 20 shapes of our 3D image database for residual values 3%, 5% and 10% (the horizontal axis is proportional to the parameter n).

  C#$DBC"A#D#%AF

Figure 16 :

 16 Figure 16: Computing times (in seconds, vert. axis) versus image sizes (horiz. axis).

Table 3 :

 3 Average Hausdorff distance between M (R θ (X)) and R θ (M (X)) (in 3d, rotations shown were performed around the Y and around the Z axis). Results shown are averages of all tested angles. Lowest values are highlighted.

Table 1 :

 1 .97 8.78 9.88 7.58 7.07 7.42 7.72 6.92 GIMA 6.30 8.05 8.97 6.77 7.75 8.46 8.77 7.56 DLMA 5.61 7.57 8.32 6.23 5.17 6.01 6.78 5.36 DLLMA 5.89 8.05 9.00 6.66 4.92 5.97 7.72 5.61 Average Hausdorff distance between M (X) and M (E(X, n)). Noise level (parameter n) is expressed as a percentage of object area, results shown are averages of all tested residuals. Lowest values are highlighted.

			2D				3D
	Noise	3%	6%	9%	all	3%	6%	9%	all
	EMA 62D				3D
	Noise	3%	6%	9%	all	3%	6%	9%	all
	EMA 1.01 1.42 1.75 1.21 0.39 0.52 0.60 0.44
	GIMA 1.32 1.83 2.20 1.53 0.80 0.99 1.11 0.85
	DLMA 1.02 1.51 1.75 1.23 0.28 0.37 0.43 0.31
	DLLMA 1.12 1.66 1.94 1.36 0.28 0.36 0.44 0.31

Table

  ,23 8,13 10,71 8.55 4,86 5,71 6,89 5.82 GIMA 6,94 9,00 13,17 8.65 5,24 5,48 6,64 5.79 DLMA 7,60 9,07 10,97 8.28 3,86 4,65 6,59 5.03 DLLMA 7,41 8,67 10,62 8.02 4,14 4,80 7,26 5.40

	2D				3D
	Residual 10% 20%	30%	all	3%	5% 10%	all
	EMA 8

Table 4 :

 4 Average dissimilarity between M (R θ (X)) and R θ (M (X)) (in 3d, rotations shown were performed around the Y and around the Z axis). Results shown are averages of all tested angles. Lowest values are highlighted.

The center of this ball is also the projection of x onto the convex hull of Π(x) (see[START_REF] Chazal | Normal cone approximation and offset shape isotopy[END_REF]).

(Hugues Talbot) This work has been partially supported by the "ANR BLAN07-2 184378 MicroFiss" project.