Compensated Convexity, Multiscale Medial Axis Maps and Sharp Regularity of the Squared Distance Function

We introduce a new stable mathematical model for locating and measuring the medial axis of geometric objects, called the quadratic multiscale medial axis map of scale $\lambda$, and prove a sharp regularity result for the squared-distance function to any closed non-empty subset $K$ of $\mathbb{R}^n$. Our results exploit properties of the function $C^l_{\lambda}(dist^2(\cdot;\, K))$ obtained by applying the quadratic lower compensated convex transform of parameter $\lambda$ to $dist^2(\cdot;\, K)$, the Euclidean squared-distance function to $K$. Using an estimate for the tight approximation of $dist^2(\cdot;\, K)$ by $C^l_{\lambda}(dist^2(\cdot;\, K))$, we prove $C^{1,1}$-regularity of $dist^2(\cdot;\, K)$ outside a neighbourhood of the closure of the medial axis $M_K$ of $K$, and give an asymptotic formula for $C^l_{\lambda}(dist^2(\cdot;\, K))$ in terms of the scaled squared distance to $K$ and to the convex hull of the set of points that realize the minimum distance to $K$. The multiscale medial axis map, $M_{\lambda}(\cdot;\, K)$, is a family of non-negative functions whose limit as $\lambda \to \infty$ exists and is called the multiscale medial axis landscape map, $M_{\infty}(\cdot;\, K)$. We show $M_{\infty}(\cdot;\, K)$ is strictly positive on the medial axis $M_K$ and zero elsewhere. We give conditions to ensure $M_{\lambda}(\cdot;\, K)$ keeps a constant height along parts of $M_K$ generated by two-point subsets with the height dependent on the distance between the generating points, so giving a hierarchy between different parts of $M_K$ that enables subsets of $M_K$ to be selected by thresholding. Given a compact subset $K$ of $\mathbb{R}^n$, while it is well known that $M_K$ is not Hausdorff stable, we prove $M_{\lambda}(\cdot;\, K)$ is stable under Hausdorff distance, and deduce implications for localization of the stable parts of $M_K$. Examples are included.


Introduction
The medial axis of an object is a geometric structure that was introduced by Blum [14] as a means of providing a compact representation of a shape. Initially defined as the set of the shock points of a grass fire lit on the boundary and allowed to propagate uniformly inside the object, closely related definitions of skeleton [17] and cut-locus [53] have since been proposed, and have served for the study of its topological properties [3,22,41,44,51], its stability [23,21] and for the development of fast and efficient algorithms for its computation [1,12,11,39,46]. Applications of the medial axis are ample in scope and nature, ranging from computer vision to image analysis, from mesh generation to computer aided design. We refer to [50] and the references therein for applications and accounts of some recent theoretical developments.
An inherent drawback of the medial axis is, however, its sensitivity to boundary details, in the sense that small perturbations of the object (with respect to the Hausdorff distance) can produce huge variations of the corresponding medial axis. This observation has prompted a large body of research that has roughly followed two lines, both aimed at the definition of some stable modification of the medial axis: one consists of reducing the complexity of the medial axis by pruning the less important parts of the domain [49], the other considers the definition of filter conditions that identify subsets of the medial axis which are stable to perturbations of the sets and retain some of its topological properties, for instance, homotopy equivalence with the object. Within this second line of research, we mention, among others, the λ−medial axis introduced in [20], the θ−homotopy preserving medial axis introduced initially in [31] and subsequently modified in [52] to ensure the homotopy equivalence, and the power crust method [8]. The λ-medial axis and the θ-homotopy preserving medial axis are explicitly defined as subsets of the medial axis, being collections of those points of the medial axis that meet some geometrical criteria. Such criteria are expressed in terms of a bound either on the distance to the boundary of the object or on the separation angle θ (see Definition 3.16 below), respectively. The power crust, in contrast, and also the algorithm discussed in [25,26], provide a continuous approximation of the medial axis constructed using a subset of the vertices, called poles, of the Voronoi diagram of a finite point sample of the object boundary. In all such works, the stable modifications are sought by identifying directly points of the medial axis or of an approximation of it. The excellent survey paper [10] contains a thorough discussion of such approaches and of the related stability issues.
We adopt in this paper a fundamentally distinct strategy which, if compared with the works mentioned above, represents an indirect approach relying on the use of the compensated convex transforms [59]. The theory of compensated convex transforms has been introduced and applied in the calculus of variations for finding the quasiconvex envelope of a function [55,56,57,58] and for finding tight smooth approximations of the maximum function and the squared-distance function [60]. Compensated convex transforms, however, also provide a natural and stable global method to extract geometric singularities, such as ridges, valleys and edges, from a given function by manipulating its 'landscape' [62,63], and it is in this way that the transforms, in particular the lower compensated convex transform (hereafter, called also the lower transform), will be used in this paper. Whether one applies the lower compensated convex transform or the upper compensated convex transform depends on the type of geometric singularities to be extracted. The works [62,63] present a systematic study on the use of these basic transforms to extract singularities from the graph of functions in general, or from the characteristic functions of compact sets, whereas the patent application [61] contains various applications including our method for extracting the multiscale medial axis map. The key properties that are exploited to highlight and/or to design a specific singularity are: the tight approximation of the compensated transforms, their regularity and the manner in which they respond to the type of curvature. More specifically, [62] focuses on the basic use of these transforms to detect ridges, valleys and saddle points of graph of functions, whereas [63] presents the design of a transform which is capable of filtering out the 'regular points' and the 'regular directions' on manifolds.
The application of the lower transform to study the medial axis of a set is motivated funda-mentally by the identification of the medial axis with the singularity set of the distance function [37,Lemma 8.5.12] and by the geometric structure of this set [3,18,43]. On the other hand, the distance function, its regularity and its geometric structure, are well-studied both in geometric measure theory [30] and in the theory of partial differential equations [19,29,34,37]. If the set K is a smooth compact submanifold of R n , there are many local regularity results of the distance function near K [40,32,27,28], whereas, for a general bounded open set Ω, some results by Albano [2] imply that the distance function dist(·; Ω c ) is locally C 1,1 in Ω \ M Ω c in the sense that if x ∈ Ω \ M Ω c , there is a δ > 0 such that dist(·; Ω c ) ∈ C 1,1 (B(x; δ)).
In the following, however, it is more convenient to refer to the squared-distance function and use the identification of the singular set of the distance function with the set of points where the squared distance function fails to be locally C 1,1 . Here, we just note that the advantage of referring to the squared distance function rather than to the distance function has also been realized in other contexts, such as, in the study of the motion of surfaces by its mean curvature represented by manifolds with codimension greater than one [24,5]. We refer to [4] for a detailed study on the properties of the squared distance function and on its applications in the geometric evolution problems.
Using properties of the lower transform, we apply the lower compensated convex transform to the Euclidean squared-distance function which gives a smooth (C 1,1 ) tight approximation outside a neighbourhood of the closure of the medial axis (see Theorem 3.3), and define our multiscale medial axis map as a scaled difference between the squared-distance function and its lower transform. From the property of the tight approximation of the lower transform of the squared-distance function, we also deduce a sharp C 1,1 -regularity result (see Corollary 3.8 and Example 3.10) of the squareddistance function outside a neighbourhood of the closure of the medial axis of K, which can be viewed as a weak Lusin type theorem for the squared-distance function and extend regularity results of the squared-distance function to any closed non-empty subset of R n . This result also offers an instance of application of the compensated convex transform to obtain a fine result of geometric measure theory and is, somehow, related to the behaviour of semiconcave functions (see [19] and Remark 2.7(c) below). We observe that, in general, the regularity of dist 2 (·; K) cannot be better than C 1,1 even for a compact convex set K where M K = ∅. A simple example is given by the square K = [0, 1] × [0, 1] ⊂ R 2 . In this case, it can be easily verified that dist 2 (x; K) is globally C 1,1 but not C 2 .
The application of the lower compensated convex transform of scale λ to the squared-distance function produces a continuous function in R n that remains strictly positive on the medial axis and tends to zero outside of it as a positive parameter λ becomes very large (see Proposition 3.20). We will, in fact, characterize the limit of the multiscale medial axis map of scale λ as λ approaches to infinity (see Theorem 3.23) and refer to this geometric structure as the quadratic multiscale medial axis landscape map of K. The values of this map are well separated, in the sense that they are zero outside the medial axis and remain strictly positive on it. Furthermore, we will give conditions (see Proposition 3.18 and Section 5) that ensure that the multiscale medial axis map of scale λ actually keeps a constant height along the parts of the medial axis generated by two-point subsets, with the value of the height dependent on the distance between the two generating points. Such values can, therefore, be used to define a hierarchy between different parts of the medial axis and we can thus select the relevant parts through simple thresholding, that is, by taking suplevel sets of the multiscale medial axis map. To reflect this property, we use the word "multiscale". For each branch of the medial axis, the multiscale medial axis map automatically defines a scale associated with it. In other words, a given branch has a strength which depends on some geometric features of the part of the set that generates that branch.
Given a closed non-empty subset K of R n , we will also prove that, despite the medial axis of K not being Hausdorff stable, the quadratic multiscale medial axis map, is indeed Hausdorff stable (see Theorem 4.3). It follows that the graph of the medial axis map carries more information than the medial axis itself, which allows the definition of a hierarchy between the parts of the medial axis and the selection of the relevant ones through simple thresholding, that is, by taking suplevel sets of the medial axis map. In this manner, it is possible to choose the main parts that reflect genuine geometric features of the object and remove minor ones generated by noise.
In conclusion, we observe that while our method seems to share similarities with those based on the extraction of ridges of the distance transform [9,16,39,45,54], (that require, however, an a-priori definition of ridge, based usually on an approximation of the derivative of the distance transform), the method we propose is, in fact, substantially different from such approaches, given that we obtain a neighbourhood of the singularities as the difference between the squared-distance transform and its smooth tight approximation. In this manner, as mentioned above, we provide an indirect definition of the singularity, which does not require any derivative approximation or any differentiability assumption.
After this brief introduction, the next section will introduce the relevant notation and recall basic results in convex analysis and lower compensated convex transforms. Section 3 contains the definition of the multiscale medial axis map, and some of its principal properties, such as the tight approximation of the lower compensated transform to the squared-distance transform (see Theorem 3.3) and as an application, we deduce a sharp regularity result of the squared-distance function to any non-empty closed subset of R n (see Corollary 3.8). Section 4 presents the Hausdorff stability of the multiscale medial axis map whereas Section 5 discusses some mathematical prototype models of explicitly calculated medial axis maps for a simple four point set to some more complicated three dimensional objects. Finally, Section 6 concludes the paper with the proofs of the main results.

Notation, Basic Definitions and Preliminary Results
Throughout the paper R n denotes the n-dimensional Euclidean space, and |x| and x·y the standard Euclidean norm and inner product, respectively, for x, y ∈ R n . In some cases, we will also make use of the notation (x, y) to denote the point of R n given by xe 1 + y 2 e 2 + . . . + y n e n , where {e 1 , . . . , e n } is an orthornormal basis of R n , (x, y 2 , . . . , y n ) ∈ R n and y = y 2 e 2 + . . . + y n e n . Given a non-empty subset K of R n , K c denotes the complement of K in R n , i.e. K c = R n \ K, K its closure and co[K] the convex hull of K, that is, the smallest (with respect to inclusion) convex set that contains the set K. For x ∈ R n and r > 0, B(x; r) indicates the open ball with center x and radius r whereas S(x; r) denotes the sphere with center x and radius r and is the boundary of B(x; r). The distance transform of a non-empty set K ⊂ R n is the function that, at any point x ∈ R n , associates the distance of x to K, which is defined as inf{|x − y|, y ∈ K} and is denoted as dist(x; K). We use the notation Df to denote the derivative of f . Across the current literature, there is no uniform definition of the medial axis, with its meaning changing from one author to another. What the medial axis is for one, becomes the skeleton for another, and in some cases subtle differences are present, especially in the continuum case, where the closure of such sets is considered. In this paper we adopt the definition given by Lieutier in [41], but it is reformulated here to include a non-empty closed set K ⊂ R n with K = R n as well as a non-empty bounded open set Ω.
Definition 2.1. For a given non-empty closed set K ⊂ R n , with K = R n , we define the medial axis M K of K as the set of points x ∈ R n \ K such that x ∈ M K if and only if there are at least two different points y 1 , y 2 ∈ K, satisfying dist(x; K) = |x − y 1 | = |x − y 2 |. For a non-empty bounded open set Ω ⊂ R n , the medial axis of Ω is defined by M Ω := Ω ∩ M ∂Ω . (b) Other frequently used notions are those of the skeleton of K, denoted as skl(K), and the cut locus of a manifold, denoted as cl(K), which applies to the more general case of Riemannian geometry. Here, for a non-empty closed subset K of the Euclidean space R n , we define the skeleton of K to be the set of the centers of the maximal (with respect to inclusion) open balls contained in K c , whereas the cut locus of K is taken to be the set of the cut locus of the points of the boundary of K in R n , where the cut locus of p in ∂K is the set of points in the manifold where the geodesics starting at p stop being minimizing. It can then be shown that cl(K) = M K . As a result, the notions of medial axis, skeleton and cut locus are related but are not the same; for example, in general, [44,25] Our definition of the medial axis M K is, in particular, consistent with our main convergence result (Theorem 3.23) which recovers M K as the set on which the limit of the medial axis map M λ (·; K) (Definition 3.1) as λ tends to infinity is strictly positive.
Next we collect definitions and results from convex analysis for functions f taking finite values, i.e. for f : R n → R, which will be used in this paper, and refer to [36,47] for details and proofs.
Given a function f : R n → R bounded below, the convex envelope co[f ] is the largest convex function not greater than f . We will often make use of the following characterization. Proposition 2.3. Let f : R n → R be coercive in the sense that f (x)/|x| → ∞ as |x| → ∞, and x 0 ∈ R n . Then (i) The value co [f ] (x 0 ) of the convex envelope of f at x 0 ∈ R n is given by If, in addition, f is lower semicontinuous, the infimum is reached by some (λ * i , x * i ) for i = 1, 2, . . . , n + 1 with (x * i , f (x * i ))'s lying in the intersection of a supporting plane of the epigraph of f , epi(f ), and epi(f ).
(ii) The value co [f ] (x 0 ), for f taking only finite values, can also be obtained as follows: with the sup attained by an affine function * ∈ Aff(R n ).
We then introduce the following definition [13], which is needed in Lemma 2.5.
Definition 2.4. Assume x 0 ∈ R n . We say that a function f : R n → R is upper-semidifferentiable at x 0 if there exists a ∈ R n such that The following lemma, concerning the existence and properties of an optimal affine function, will be needed for the proofs of Proposition 2.14 and Theorem 3.3.
Definition 2.6. Let f : R n → R be a lower semicontinuous function [47,36] satisfying for some constants A 1 , A 2 ≥ 0. The (quadratic) lower compensated convex transform (lower transform for short) for f with scale λ > A 1 is defined for x ∈ R n by If f is bounded below, we may set λ ≥ 0.
Remark 2.7. (a) The requirement of the lower semicontinuity of f is to guarantee that C l λ (f )(x) → f (x) as λ → ∞ for all x ∈ R n , since otherwise, the lower transform will converge to the lower semicontinuous envelope of f . (c) Recalling from [19] that a function f : R n → R is called c−semiconvex if, for some constant c > 0, the function f (x) + c/2|x| 2 is convex, we observe that the lower compensated convex transform for f with scale λ, C l λ (f ) is a 2λ-semiconvex function. In fact, C l λ (f ) represents the 2λ-semiconvex envelope of f . We sometimes use such a property to extend some properties of semiconvex functions to C l λ (f ).
(d) To gain further geometric insight into the lower compensated convex transform defined by (2.6), in Figure 1 we display the steps of the construction of C l λ (f ) for f (x) = dist 2 (x; K) with K = {−1, 1} and λ = 2. The graph of the augmented function f + λ|x| 2 is displayed in Figure 1(b) along with f , whereas Figure 1(c) shows the convex envelope of the augmented function. Figure 1(d) displays finally the graph of C l λ (f ) which is compared with that of f . Note that the convex envelope of the augmented function is different from f + λ|x| 2 only in a neighbourhood of the singular point 0 of f , so that, when then we subtract the weight, the final effect is a smoothing of f only in such neighbourhood. This simple example, along with the ones discussed in Section 5, enables one also to understand the role of the parameter λ and our meaning of scale. The parameter λ acts as a scale parameter in the sense that it controls the curvature of the lower compensated convex transform in the neighbourhood of the singularity of the function and allows the extraction of the singularity with a value which gives somehow a measure of its strength. Also one may observe the so-called 'tightness' of the lower compensated convex transform approximation of the original function from below (see Proposition 2.10), which agrees with the original function except in the neighbourhood near the singular point. (c) Graph of the convex envelope of f + λ| · | 2 compared to that of f + λ| · | 2 ; (d) Graph of C l λ (f ) compared to that of f .
The following properties of C l λ (f ) will also be used. Proposition 2.8. Given f : R n → R that satisfies (2.5), then for all A 1 < λ < τ < ∞, we have

8)
and, for λ > A 1 Proposition 2.9. If f ≤ g in R n and satisfy (2.5), then The transform C l λ (f ) realizes a 'tight' approximation of the function f , in the following sense (see [59,Theorem 2.3(iv)]).
Such a property motivates the definition of the multiscale ridge transform which was introduced in [62] to extract ridges of general functions and shown to be invariant with respect to translation. This multiscale ridge transform will be used in Section 3 to define the multiscale medial axis map (see Definition 3.1).
Definition 2.11. Given λ > 0, the ridge transform of scale λ, for a given function f : R n → R satisfying (2.5), is defined as: We now present some regularity properties of C l λ (dist 2 (·; K)), which will be exploited to analyze the behaviour of the multiscale medial axis map. We recall first the following result given in [59,Lemma 4.3].
The next property is a useful inequality for the derivative of the lower transform, DC l λ (dist 2 (·; K)).
We next introduce the sets K(x) and K 2,λ (x), which will be used to gain insight into the geometric structure of C l λ (dist 2 (·; K)).
Definition 2.16. Let K ⊂ R n be a non-empty closed set. For any x ∈ R n , let r(x) = dist(x; K). We then define the following sets: and for λ > 0, The following result, obtained in the proof of [59, Theorem 3.7], gives an explicit expression of the lower transform of dist 2 (·; K(x)), the squared distance to the set K(x), and will be used to produce a bound on the multiscale medial axis map (see Theorem 3.15(i)).
Proposition 2.18. Let K ⊂ R n be a non-empty closed set and M K the medial axis of K. Assume x ∈ M K and denote by K(x) and K 2,λ (x) the sets defined by (2.15) and by (2.16), respectively. Then, for all y ∈ R n , We will also need, for the proof of Theorem 3.23, the following explicitly calculated formula of the lower transform for compact sets contained in a sphere S(0; r) = {x ∈ R n , |x| = r} centred at 0 ∈ R n with radius r > 0. The formula is easy to derive following similar calculations to those in the proof of [60, Theorem 1], or of [59, Theorem 5.1].
Proposition 2.19. Let K ⊂ S(0; r) be a non-empty compact set. Then for every x ∈ R n , We will invoke the following technical lemma several times (see Lemma 3.2 in [59]).
In the next lemma, which generalizes slightly [59, Lemma 3.3], we give the expression of the lower transform of the squared distance to a set of two points. The two points, without loss of generality, are assumed to lie along a basis vector of R n , specifically, along the basis vector e 1 ∈ R n . This lemma will be used extensively when we investigate the behaviour of the multiscale medial axis map with respect to perturbations of the boundary of K.
Lemma 2.21. Assume n ≥ 2 and let {e 1 , . . . , e n } be an orthonormal basis of the Euclidean space R n . Let K = {−αe 1 , αe 1 }, where α > 0. We write y = e 2 y 2 + · · · + e n y n ∈ R n−1 and represent the point xe 1 + y ∈ R n as the pair (x, y), which therefore denotes the point (x, y 2 , . . . , y n ) ∈ R n . Then for every λ > 0, we have (2.20) In particular, We conclude this section with the definition of δ−neighbourhood of a set, of Hausdorff distance between two sets [6], and of −sample of a set [10].
Definition 2.22. Given a non-empty subset E of R n and δ > 0, we define the δ-neighbourhood E δ of E by Note that E δ is an open subset of R n .
Definition 2.23. Let E, F be non-empty subsets of R n . The Hausdorff distance between E and F is defined in [6] by This definition is also equivalent to saying that Definition 2.24. Let K be a compact subset of R n . A sample S of the boundary of K is a finite set of points of the boundary of K, i.e. S ⊂ ∂K and #(S) ∈ N where #(S) denotes the cardinality of S. An −sample of ∂K is a sample whose Hausdorff distance to ∂K is less than , that is, where the diameter of K, diam(K), is defined as

The Multiscale Medial Axis Map
In this section, we define the quadratic multiscale medial axis map M λ (·; K), characterize some of its properties, and establish its relation to the medial axis M K . As a by-product, we also infer sharp regularity results for the squared distance function dist 2 (·; K), which are of independent interest.
Definition 3.1. Let K ⊂ R n be a non-empty closed set. The quadratic multiscale medial axis map of K (medial axis map for short) with scale λ > 0 is defined for x ∈ R n by For a bounded open set Ω ⊂ R n with boundary ∂Ω, we define the quadratic multiscale medial axis map of Ω with scale λ > 0 as The convergence of the lower transform to the original function as λ → ∞ yields that lim λ→∞ R λ (dist 2 (·; K))(x) = 0, implying that the values of the ridge transform can be very small when λ > 0 is large. To make the height of our medial axis map on the medial axis bounded away from zero, we thus need to scale the ridge transform. The factor (1 + λ) turns out to be the "right" scaling factor, as will be justified in Theorem 3.15 below, where it will be shown that on the medial axis M K , the medial axis map M λ (x; K) is bounded both above and below by quantities independent of λ.
(b) The quadratic multiscale medial axis map can also be seen as a morphological operator [48], equal to the scaled top-hat transform of the squared distance transform with quadratic structuring function. Letting f (x) = dist 2 (x; K) and b λ (x) = −λ|x| 2 , it can be shown that the lower transform corresponds to the grayscale opening operator with quadratic structuring function [62]; i.e., Notwithstanding such an interpretation, it is convenient to view Definition 3.1 in terms of the lower compensated convex transform. The exploitation of properties of such transforms permits a relatively easy evaluation of the geometrical properties of M λ (·; K) and also permits an easy numerical realization of M λ (·; K). This relies on the availability of numerical schemes for computing the lower transform of a given function, which entails the availability of schemes to compute the convex envelope of a function. We refer to [64] for the algorithmic and implementation details of the schemes for realizing the lower transform of a function.
We begin with a key quantitative estimate of the tight approximation of the squared distance function dist 2 (·; K) by its lower transform C l λ (dist 2 (·; K)). This result not only underpins our study of the rôle of M λ in characterizing the medial axis M K , but also yields improved locality and regularity properties of C l λ (dist 2 (·; K)) and dist 2 (·; K) respectively (see Corollaries 3.6, 3.8 and 3.13) which are of interest in their own right. Theorem 3.3. Let K ⊂ R n be a non-empty closed set and denote by M K the medial axis of K.
Now assume λ > 0 and introduce the set Clearly M K ⊂ V λ, K , so this defines a "neighbourhood" of the medial axis M K of K (note that it is possible that M K ∩ K = ∅, so V λ, K is not necessarily a neighbourhood in the strict sense), and V λ, K is a closed set. Moreover, as λ > 0 increases, V λ,K describes a family of shrinking sets such that and if we take the support of the multiscale medial axis map, Theorem 3.3 yields that With the help of (3.8), we can show the following result that characterizes M K in terms of spprt(M λ (·; K)).
An important consequence of Theorem 3.3 is the following locality property of the lower transform of the squared distance function. This result is also of independent interest, in particular because it quantifies the size of neighbourhood needed to evaluate C l λ (dist 2 (·; K)), and also because it will be exploited in the proofs of Theorems 3.15 and 3.15 to establish results characterizing the properties of M λ (·; K). Corollary 3.6. (Locality Property) Suppose K ⊂ R n is a non-empty closed set. Then for every Remark 3.7. (a) Corollary 3.6 improves the result in [38], where the radius of the ball for the locality property is r( (b) In [38], it was also established that x 0 ∈ R n is a stationary point of C l λ (dist 2 (·; K)) if and only if x 0 ∈ co[K(x 0 )]. We will see that the 'only if ' part of this result is a consequence of arguments from the proof of Theorem 3.3: see Remark 6.1. Theorem 3.3 can also be combined with Proposition 2.13 to yield a regularity property of the distance transform, which can be viewed as a weak version of the Lusin theorem for the squareddistance function.
Corollary 3.8. Assume λ > 0. Let K ⊂ R n be a non-empty closed set and V λ,K the neighbourhood of M K defined by (3.6). Then Remark 3.9. It follows from Corollary 3.8 and Stepanov's Theorem [35,42] that the Hessian of dist 2 (·; K) exists almost everywhere in R n \ M K .
Both estimate (3.3) in Theorem 3.3 and estimate (3.12) for the Lipschitz constant in Corollary 3.8 (when λ ≥ 1) are, in fact, sharp, as the following example shows.
We cover next the case of a bounded open set Ω ⊂ R n , giving a precise statement about equality of the medial axis maps M λ (x; Ω c ) and M λ (x; ∂Ω), followed by a modification of Corollary 3.8.
Proposition 3.11. Suppose Ω ⊂ R n is a non-empty bounded open set and let λ > 0. Then and consequently, (3.14) Remark 3.12. (a) Property (3.14) of the medial axis map is important in many practical situations. For example, in image processing, the objects Ω of which we wish to find the medial axis might be defined by taking a threshold from a greyscale image, that is, as a suplevel set of the image function. The object is then represented by a binary image rather than by its boundary. Therefore, in this case, it might be more convenient for us to compute numerically the medial axis map M λ (x; Ω c ) rather than M λ (x; ∂Ω).
(b) It is worth noting the different qualitative behaviour of the convex envelope and the compensated convex transform that appears in (3.13). For λ > 0, the left hand side of (3.13) is always positive in R n \ Ω whereas the right hand side equals zero in Ω c . By setting λ = 0, the left hand side of (3.13) reduces to the the convex envelope of dist 2 (x; ∂Ω), which vanishes in the convex hull co[Ω] of the closure of Ω, whereas the right hand side of (3.13) gives the convex envelope of dist 2 (x; Ω c ), which is identically zero in R n .
For a bounded open set Ω, Corollary 3.8 modifies as follows.
Corollary 3.13. Let Ω ⊂ R n be a bounded non-empty open set. Then

16)
and diam(Ω) is the diameter of Ω. Furthermore, for all x, y ∈ Ω \ W λ,Ω , A consequence of Corollary 3.13 is that outside any neighbourhood of W λ,Ω , dist 2 (·; Ω c ) is a C 1,1 function. However, we also notice that M K can have positive n-dimensional Lebesgue measure. Therefore the measure of W λ,Ω might not be small even when λ > 0 is large. Corollary 3.8 and Corollary 3.13 also demonstrate that the lower transform can be viewed as a C 1,1 extension of the squared-distance function from the set V c λ,K , on which dist 2 (·; K) = C l λ (dist 2 (·; K)), to R n and from Ω \ W λ,Ω to Ω, respectively. Theorem 3.3 showed that if x ∈ M K , then M λ (x; K) = 0 when λ is sufficiently large. We now further explore the relationship between the medial axis map M λ (·; K) and the medial axis M K , both establishing λ-independent positive upper and lower bounds on M λ (x; K) whenever x ∈ M K , and fully characterizing the limit of M λ (·; K) as λ → ∞. The following geometric structure will play a key rôle in both Theorem 3.15 and Theorem 3.23.
Definition 3.14. Let K ⊂ R n be a non-empty closed set and for x ∈ R n , denote by K(x) the set defined by (2.15), that is, K(x) = B(x; r(x)) ∩ K, and denote by co[K(x)] the convex hull of K(x). The quadratic multiscale medial axis landscape map of K is defined for x ∈ R n by Theorem 3.15. Let K ⊂ R n be a non-empty closed set, and denote by M K the medial axis of K and by M ∞ (x; K) the quadratic multiscale medial axis landscape map defined by (3.18).
(i) For every λ > 0 and every (ii) For every λ > 0 and for every x ∈ R n , The lower bound in (3.19) can be expressed in terms of the separation angle which has been used, for instance, in [52], for a local geometrical characterization of the medial axis.
Definition 3.16. Let K ⊂ R n be a non-empty closed set and denote by M K the medial axis of K. For x ∈ M K , let y 1 , y 2 ∈ K(x) and denote by ∠[y 1 − x, y 2 − x] the angle between the two non-zero vectors y 1 − x and y 2 − x, taken between 0 and π, i.e. ∠[y 1 − x, y 2 − x] = cos −1 (y 1 −x)·(y 2 −x) |y 1 −x||y 2 −x| . We then define the separation angle θ x for x ∈ M K as follows: Since for x ∈ M K , there are at least two different points y 1 , y 2 ∈ ∂B(x; r(x)) ∩ K, it follows that θ x > 0, and hence, comparing with This implies that the limit of M λ (x; K) can extract exactly the medial axis of K. If we apply a slightly weaker scaling to the ridge transform, say (1 + x) α for 0 < α < 1, and define M α λ (x; that is, M α λ (x; K) approaches the indicator function of M K [36] as λ becomes large.
We can now characterize the limit of M λ (x; K) as λ → +∞ for all x ∈ R n .
Theorem 3.23. Suppose K ⊂ R n be a non-empty closed set. Then for every x ∈ R n , From Theorem 3.3, Corollary 3.8 and Theorem 3.23, it follows that when λ > 0 is increasing, the support of M λ (·; K) is contained in a shrinking neighbourhood of M K and approaches the multiscale medial axis landscape map M ∞ (x; K). The numerical advantage of studying M λ (·; K) as an approximation of the multiscale medial axis landscape map M ∞ (·; K) is that it relies only on the computation of the lower compensated convex transform of the squared distance transform, whose construction is local by virtue of Corollary 3.6, whereas the computation of M ∞ (·; K) is difficult because we need to evaluate the convex hull co[K(x)]).
Remark 3.24. A further consequence of Theorem 3.23 is that for every fixed x ∈ R n and for every non-empty closed set K ⊂ R n , the family of lower transforms λ → C l λ (dist 2 (·; K))(x) is 'differentiable' at infinity. If we let = 1/λ, g( ; x) = C l 1/ (dist 2 (·; K))(x), and g(0; x) = lim →0 g( ; x) = dist 2 (x; K), then so we have the asymptotic expansion is not continuous in R n as x approaches the medial axis M K . But from Theorem 3.23, we can show that for every x ∈ R n , using the equality For large λ > 0, (3.27) can be viewed as an approximation of dist 2 (x; co[K(x)]) by continuous functions. The function dist(x; co[K(x)]) has been used, for instance, for surface reconstruction when K ⊂ R 3 is finite [25]. While it is difficult in general to calculate dist 2 (x; co[K(x)]) directly, we see from (3.27) that the numerical computation of C l λ (dist 2 (·; K))(x), whose evaluation involves only local convex envelope calculations because of Corollary 3.6, offers an easy approximation of dist 2 (x; co[K(x)]).
We conclude this section by observing briefly that, based on the estimates of Theorem 3.15 and Proposition 3.18, it is reasonable to define an alternative medial axis map by taking the square root of M λ (x; K). Definition 3.26. We define the multiscale medial axis map of linear growth (linear medial axis map for short) by From Proposition 3.18, we obtain that the height of this linear medial axis map is 'proportional' to the distance function itself; that is, for λ > 0, we have Note that the linear medial axis map M 1 λ (x; K) is different from a definition based on the lower compensated convex transform for the distance function dist(x; K) itself, i.e. based on R λ (dist(·; K)). Of course, we can define such maps using the p-distance function dist p (x; K) for any 1 ≤ p < ∞. But in this paper, we focus mainly on the medial axis map M λ (·; K) defined using the squared-distance function, i.e., for p = 2, in which case the geometry of M λ (x; K) is easy to control. For instance, as we will see in the next section, M λ (x; K) has the same height along the parts of the medial axis generated by two points. This is a key property when one looks for approximate medial axes by applying the Voronoi diagram method of finite -samples.

Hausdorff Stability
Quantifying the instability of the medial axis is of fundamental importance for both theory and computation. This aspect becomes more and more relevant in practice nowadays, given that point clouds are increasingly being used for geometric modeling over a wide range of applications. Moreover, there are computational approaches, such as the Voronoi diagram method, which search for a continuous approximation of the medial axis of a shape starting from subsets of the Voronoi diagram of a sample of the shape boundary. The presence of noise on the boundary, and/or the discrete character of samples of the boundary shape thus call for methods that permit the control of the parts of the medial axis which are not stable. In this section, we will discuss how this aspect is tackled by the multiscale medial axis map. In the first part of the section, we examine the values of M λ (·; K) when the distance of the point to the boundary of the set is achieved by two points, whereas in the second part we discuss the Hausdorff stability of M λ (·; K). Proposition 4.1. Assume n ≥ 2 and let {e 1 , . . . , e n } be an orthonormal basis of the Euclidean space R n . Let K = {−αe 1 , αe 1 }, where α > 0. We write y = e 2 y 2 + · · · + e n y n ∈ R n−1 and represent the point xe 1 + y ∈ R n as the pair (x, y) ∈ R × R n−1 , which therefore denotes the point (x, y 2 , . . . , y n ) ∈ R n . Then for every λ > 0, we have (4.1) Remark 4.2. The medial axis map M λ ((x, y); K) reaches its maximum on the medial axis of K, at the point (0, y) = y 2 e 2 +· · ·+y n e n ∈ R n−1 , attaining the value M λ ((0, y); K) = α 2 . Note that α > 0 is half the distance between the two points −αe 1 and αe 1 of K. Another important observation is that M λ ((x, y); K) is a function of the x-variable only, and does not change its height along the y direction. Therefore, on branches of the medial axis generated by two points, the height remains the same. If α is small (equivalently, the two points in K are close to each other), the values of M λ ((x, y); K) will be uniformly small.
We next give the Hausdorff stability property of the multiscale medial axis map, followed by some comments on implications of this property for the localization of the medial axis of a domain. Theorem 4.3. Assume λ > 0. Let K, L ⊂ R n be non-empty compact sets. Then as L → K under the Hausdorff distance, M λ (·; L) → M λ (·; K) uniformly in every fixed bounded set in R n . More precisely, if we let µ := dist H (K, L) be the Hausdorff distance between K and L, then for x ∈ R n , C l λ (dist 2 (·; K))(x) − C l λ (dist 2 (·; L))(x) ≤ µ (dist(x; K) + µ) 2 + 1 + µ ,
which shows that the medial axis map is uniformly continuous on compact sets with respect to the Hausdorff metric.
(b) While the medial axis of K is not a stable structure with respect to the Hausdorff distance, its medial axis map M λ (x; K) is by contrast a stable structure. This result complies with (4.3) which shows that as λ becomes large, the bound in (4.3) becomes large.
As an immediate consequence of Theorem 4.3 we have the following result, which relates the medial axis map of the boundary of a domain Ω with that of its -samples K .
Remark 4.6. (a) If we consider an -sample K of ∂Ω, that is, a discrete set of points such that dist H (∂Ω, K ) ≤ , Corollary 4.5 yields a simple criteria that permits the suppression of those parts of the Voronoi diagram of K that are not related in the limit, as → 0, to the stable parts of the medial axis of Ω.
(b) Since the medial axis of K is the Voronoi diagram of K , if V denotes the set of all the vertices of the Voronoi diagram Vor(K ) of K , and P is the subset of V formed by the poles of Vor(K ) introduced in [7], (i.e. those vertices of Vor(K ) that converge to the medial axis of Ω as the sample density approaches infinity), then as a result of Proposition 3.20, for λ > 0, we conclude that Since as → 0+, K → ∂Ω, and knowing that P → M Ω [8,15], then on the vertices of Vor(K ) that do not tend to M Ω , M λ (x ; K ) must approach zero in the limit because of Proposition 3.20. As a result, in the context of the methods of approximating the medial axis starting from the Voronoi diagram of a sample (such as those described in [8,25,26,50]), the use of the multiscale medial axis map offers an alternative and much easier tool to construct continuous approximations to the medial axis with guaranteed convergence as → 0+.
With the aim of giving insight into the implications of the Hausdorff stability of M λ (x; ∂Ω) and Corollary 4.5, we display in Figure 2 the graph of the multiscale medial axis map of a nonconvex domain Ω and of an -sample K of its boundary. Inspection of the graph of M λ (x; ∂Ω) and M λ (x; K ), displayed in Figure 2(a) and Figure 2(b), reveals that both functions take comparable values along the main branches of M Ω . Also, M λ (x; K ) takes small values along the secondary branches, generated by the sampling of the boundary of Ω. These values can therefore be filtered out by simple thresholding so that a stable approximation of the medial axis of Ω can be computed. This can be appreciated by looking at Figure 2(d), which displays a suplevel set of M λ (x; K ) that appears to be a reasonable approximation of the support of M λ (x; ∂Ω) shown in Figure 2

Examples of Exact Medial Axis Maps and Their Supports
In this section we illustrate the behaviours of our multiscale medial axis map for some 2d geometric objects K, for which it is possible to obtain an explicit analytical expression for M λ (x; K). Thanks to the translation and the partial rotation invariance property of the convex envelope [62, Proposition 2.3, 2.10], it is then possible to derive an explicit analytical expression for M λ (x; K) in the case that K is a 3d-solid obtained by, for instance, rotations or translations of the models considered in this section. For the sake of conciseness, we leave the derivations to interested readers.
Though the derivation here is limited only to 2d geometric models, these models retain, nevertheless, their basic geometric features, because they are able to show that M λ (x; K) can, in fact, provide an accurate and stable way to find M K , the medial axis of K, and represents likewise an effective tool to analyze the geometry and structure of M K . We will also see how it is possible to select either the main stable parts of M K or to locate its fine parts by using suplevel sets of M λ (x; K). and, after some lengthy calculations based on the construction of affine functions, we can show that the lower transform can be expressed as follows where the auxiliary function g = g(x, y) is a continuous piecewise quadratic function defined as follows The multiscale medial axis map M λ (x; K) is then computed using the definition (3.1). In particular, for this example, after some algebraic rearrangements, it is possible to show that since all four points in K lie on a circle centerd at the origin, the medial axis map of K can be expressed as M λ ((x, y); K) = (1 + λ) 2 dist 2 ((x, y); K/(1 + λ)) − dist 2 ((x, y); co(K)/(1 + λ)) (x, y) ∈ R 2 , (5.4) where for Z ⊂ R n and α ∈ R \ {0}, we use the notation Z/α to denote the set {w ∈ R n : w = z/α for z ∈ Z}. By a closer inspection of (5.4), we can make then the following observations: The 'thickness' of the support for the main branch y-axis is, therefore, 2b/(1 + λ) while that for the minor branch x-axis is 2 b/(1 + λ).   with r > 0, whose results will be used to construct the multiscale medial axis map of a rectangular domain. By inspection, we can easily infer that whereas the lower transform, obtained after lengthy calculations based on the construction of affine functions, is given, for (x, y) ∈ Ω s by where the auxiliary function g = g(x, y) is a continuous piecewise quadratic function defined as follows The multiscale medial axis map of Ω s , M λ ((x, y); Ω s ), is obtained by applying definition (3.1) and by taking into account (5.5) and (5.6). By exploiting properties of the lower transform with respect to symmetry and translation of axis, we can then easily obtain the analytical expression for the multiscale medial axis map of a rectangular domain. If, for instance, we consider the open bounded set Ω = − r + r 2 , r + r 2 × (−r, r), then it is not difficult to show that M λ ((x, y); (5.8) Figure 4(a) displays the support of M λ ((x, y); Ω) which is a neighbourhood of the medial axis, whereas Figure 4(b) depicts the graph of M λ ((x, y); Ω). For the points (x, y) ∈ M Ω with θ x = π, it follows that M λ ((x, y); Ω) = dist 2 ((x, y); ∂Ω), so implying that in this sense, the upper bound in  , y); Ω). Example 5.3. We consider now the oval shaped domain Ω ⊂ R 2 made by the union of two semicircles with center at the points (−r/2, 0) and (r/2, 0) and radius r/2, respectively, and the rectangle (−r/2, r/2) × (−r/2, r/2). For this domain, it is not difficult to verify that whereas the lower transform, obtained after some lengthy calculations based on the construction of affine functions, is given, for (x, y) ∈ Ω, by (5.10)  , y); Ω), for λ = 10 and of the squared distance function, dist 2 ((x, y); ∂Ω) whose zero level set gives the boundary of the domain. Only the restriction to y ≥ 0 is displayed, given that M λ ((x, y); Ω) ≤ dist 2 ((x, y); ∂Ω); (b) The oval shaped domain Ω displayed with the support of its medial axis map M λ ((x, y); Ω).
In the following example we describe the behaviour of M λ (x; K) for the case of a discrete set K, sampled from a connected set, and evaluate the structure of M λ (x; K) as the sample density approaches infinity.
Example 5.4. We consider the geometric model of the uniform sampling of two parallel lines at distance b to each other. The points are taken equally spaced over each line at distance c = b with ∈ (0, 1) measuring the sampling density, in the sense that as → 0, the sampling density on the two lines tends to infinity. The sampling of the two parallel lines is defined so that the discrete points are aligned along the x− and y−axis as displayed in Figure 6(a). For such a sample K, we can then use the results obtained in Example 5.1 for the four-point set which we refer to as K 4p . It is not difficult to show that, for (x, y) ∈ R 2 , and i ∈ Z such that |y + i b| ≤ b, M λ ((x, y); K) = M λ ((x, y + i b); K 4p ) , (5.11) where M λ ((x, y + i b); K 4p ) is the multiscale medial axis map for the four-point set discussed in the Example 5.1.  The comparison of the graphs of M λ ((x, y); K) for the two different sample densities of the same object, shows how the value of M λ ((x, y); K) along the minor branches of the medial axis of the discrete set K attenuates as → 0. Consistently with the finding obtained for the four-point set, we have that the value of M λ ((x, y); K) along the minor branches is proportional to 2 . It then follows that by setting a threshold not lower than such a value, we can single out the stable part of M K , which provides an approximation for (and in this case, is in fact coincident with) the medial axis of the two parallel lines.
In this last example, we analyze a model of perturbations of the boundary domain represented by staircase-like piecewise affine curves. This effect is very common, for instance, in digital images and is the source of unrealistic medial axis branches, which are usually not desirable. Common practice in this case is to perform a boundary smoothing prior to any image processing operation. We will show that this is not needed with the multiscale medial axis map. The fine structure of the medial axis corresponding to the irregularities of the boundary is indeed captured by the multiscale medial axis map and can be filtered out. We will verify this statement on a prototype model of this boundary domain perturbation, by showing that the height of the medial axis map on such branches can be very small if the stair like effect is small.
Example 5.5. Assume c > 0 and let us consider the set K = {(x, y) ∈ R 2 : x ≥ 0, y ≥ 0, x + y ≤ c} ∩ {(x, y) ∈ R 2 , x + y ≥ c} displayed in Figure 7, which is used as a prototype of one single step perturbation of a boundary domain. The squared distance function dist 2 ((x, y); K c )) to the complement of K and the lower transform C l λ (dist 2 (·; K c ))(x, y) are then given, respectively, by The medial axis map M λ ((x, y); K c ) is then obtained from (3.1) using (5.12) and (5.13). The graph of M λ ((x, y); K c ) for λ = 9 and step size c = 1 is shown in Figure 8(a), whereas Figure 8(b) displays its support. By inspecting the graph of M λ ((x, y); K c ), we observe that after an initial increase near the corner tip, M λ ((x, y); K c ) keeps a constant value along M K , with this value proportional to the square of the step size. It is not difficult to verify the following with uniform convergence if y = x and x ≥ 0.
Despite its simplicity, this basic model elucidates the behaviour of M λ ((x, y); K c s ) for a set K s with a stair like boundary profile as, for instance, the one displayed in Figure 9(a).
For such a set K s , it is not difficult to verify that for (x, y) with M λ ((x, y); K c ) corresponding to the one step boundary domain perturbation discussed at the beginning of this example. Figure 9(b) contains the graph of M λ ((x, y); K s ), whereas Figure 9(c) shows its support, displayed together with the set K s . The height of the ridges along M Ks depends only on the gap size c, in particular, it is proportional to c 2 . It follows, therefore, that by setting the threshold larger than c 2 , the corresponding suplevel set of M λ ((x, y); K s ) will filter out all minor branches of M Ks , generated by the step-stair like boundary.   Multiscale medial axis map as obtained by the numerical implementation of M λ ((x, y); K): (a) Support of the multiscale medial axis map M λ ((x, y); K) for the digital image of a maple leaf, for λ = 10. All the fine branches generated by the steps on the boundary are displayed.; (b) Suplevel set of M λ ((x, y); K) corresponding to a threshold equal to one displaying only stable parts of the medial axis.
As an application of these concepts, we show in Figure 10 the results of the numerical realization of M λ ((x, y); K) for the digital image of a maple leaf, where we can note the effects just discussed. In particular, Figure 10(a) depicts the support of M λ ((x, y); K) with the display of all fine branches created by the step-like irregularities of the boundary domain, whereas Figure 10(b) shows the suplevel set of M λ ((x, y); K) corresponding to a threshold equal to one which singles out only the neighbourhood of stable parts of M K .
Proof of Corollary 3.8: Since M K ⊂ V λ,K , it follows from Theorem 3.3 that dist 2 (x; K) = C l λ (dist 2 (·; K))(x) for all x in the set R n \ V λ,K , which is an open set because V λ,K is closed. Thus the result is immediate from Proposition 2.13.
Proof of Proposition 3.18: Let x ∈ M K , r(x) = dist(x; K) > 0, and denote by x 1 , x 2 ∈ B(x; r(x)) ∩ K the points of K(x) that realize the separation angle θ x at the point x.
and hence as required.
The aim of the following technical construction is to show that for x in a small neighbourhood of 0, dist 2 (x; K 0,2 ) is a lower bound for dist 2 (x; K). For δ > 0, define the closed neighbourhood K δ 0, = {(r 0 + t)y/|y|, y ∈ K 0, , 0 ≤ t ≤ δ} and note that K δ 0, is clearly a compact set. Then it can easily be proved, using a contradiction argument, that for every 0 < < 1, there exists 0 < δ ≤ 2 such that K ∩ B(0; r 0 + δ) ⊂ K δ 0, . Define also another compact set by where K 0,2 will be used to 'shadow' K δ 0, , and the unbounded closed set W ,δ = V ,δ ∪ B c (0; r 0 + δ).
Proof of Proposition 4.1: This follows from the definition of M λ ((x, y); K) and Lemma 2.21.