Medial axis and a moment map

MEDIAL AXIS AND A MOMENT MAP D. SIERSMA AND M. VAN MANEN Abstract. We use Maslov dequantization to compute the medial axis of an embedded hypersurface. 1. Introduction It is well known that the medial axis of a compact embedded submanifold M in Rn can be calculated as an approximation of a Voronoi diagram, see for instance [ABK98]. In this paper we use tropical geometry to give a new interpretation of that fact. Let M be a compact connected smooth manifold of dimension n − 1, embedded by γ : M → Rn . For every x in Rn there exists a distance function gx : M → R defined by gx (s) = kx − γ(s)k2 The medial axis Cut(M ) consists of the closure of the set of those x in Rn for which gx has a non-unique global minimum. Recall that the corner locus Corner(f ) of a convex function f is the set of points where f is not differentiable. The main statement of this paper is: Theorem 1. Let M be a compact connected smooth manifold of dimension n−1, embedded by γ : M → Rn . Write for s ∈ M : 1 def (1) fs (x) = hx, γ(s)i − kγ(s)k2 ; f (x) = sup fs (x) 2 s∈M and rh (x) = logh µZ h fs (x) M ¶ ds . Then: (1) The function f is a convex function, and satisfies the equality: (2) f (x) = lim rh (x) h→∞ (2) The corner locus of f is the medial axis Cut(M ) of M . (3) The derivative of f outside the corner locus of f maps x ∈ Rn \ Cut(M ) to the nearest point on M . Date: July 6, 2009. 2000 Mathematics Subject Classification. Primary 52B55; Secondary 57Q99, 68U05. Key words and phrases. Voronoi diagrams, medial axis. 1 2 D. SIERSMA AND M. VAN MANEN (4) The derivative of the Legendre transform of f maps CH(γ(M )) \ γ(M ) onto the medial axis. (5) The sequence of function rh (x) is monotonely increasing. (6) The gradient of rh (x) is a diffeomorphism from Rn to the interior of the convex hull of γ(M ). We know of no cases where the integrals rh (x) can be evaluated explicitly, but if one uses the integral formula and numerical integration for h sufficiently big one can plot f , as is illustrated in figure 1. Figure 1. On the left hand side the graph of f (x) = sups∈M fs (x) for γ, on the right hand side the graph of g(x) = inf s∈M 21 kx − γ(s)k2 . In the first section we recall the definition of Voronoi diagrams and relate this to affine linear functions. In the second section we proof the main theorem. 2. Voronoi diagrams and piecewise affine convex functions 2.1. The Voronoi Diagram. Take a point set {P1 , · · · , PN } ⊂ Rn . Throughout this article we assume that dim(CH({P1 , · · · , PN })) = n. We want to divide Rn according to the nearest distance principle. Define the functions: 1 (3) gi : Rn → R gi (x) = kx − Pi k2 g(x) = min gi (x) 1≤i≤N 2 Here we have used the following notation: n X 2 kxk = x2i i=1 Definition 1. For a non-void subset α ⊂ {P1 , · · · , PN } the set Vor(α) is the closure of {x ∈ Rn | gi (x) = g(x) Pi ∈ α and gj (x) > g(x) Pj 6∈ α} We put Vor(∅) = Rn . MEDIAL AXIS AND A MOMENT MAP 3 For a subset α ⊂ {P1 , · · · , PN } we have Vor(α) ⊂ ∩Pi ∈α Vor(Pi ). If the intersection of Vor(α) and Vor(α′ ) is not empty then there is a β ⊂ {P1 , · · · , PN } such that Vor(α) ∩ Vor(α′ ) = Vor(β). The union over non-void α of all Vor(α) is Rn . In this situation we can also consider the affine functions and their maximum 1 (4) fi (x) = hx, Pi i − kPi k2 ; f (x) = max fi . i 2 which are related to the distance function by 1 1 gi (x) = kxk2 − fi (x) ; g(x) = kxk2 − f (x). 2 2 Note that we use here the same notation in the finite case, as we did in the continuous case, as stated in theorem 1, formula (1). Obviously we might have defined the sets Vor(α) using the functions fi , i.e.: For a subset α ⊂ {P1 , · · · , PN } the set Vor(α) is {x ∈ Rn | fi (x) = f (x) Pi ∈ α and fj (x) < f (x) Pj 6∈ α}. The affine definition of Voronoi diagrams is a bit more elegant. The union of all Vor(α), where α contains at least two points is a closed set of codimension 1, and is called the Voronoi Diagram. It is easy to see that the function f is convex and that its corner locus is equal to the Voronoi diagram. We have thus defined the Voronoi diagram as a corner locus of a convex function. Below we will define the medial axis as a corner locus of a convex function. It follows that the non-differentiability locus of the functions g and f are the same and that the Voronoi diagram can thus be defined using only affine functions 2.2. The Legendre transform. Dual to the Voronoi diagram is the Delaunay triangulation. Definition 2. Del(α) = ½ CH(α) ∅ if Vor(α) 6= ∅ otherwise We are also interested in defining the Delaunay triangulation in a very different way. Definition 3. The Legendre transform of a convex function f , with domain D ⊂ Rn is fˆ(ξ) = sup (hξ, xi − f (x)) . x∈D When the supremum does not exists, we put fˆ(ξ) = ∞. The domain Dom(fˆ) of fˆ are those ξ for which fˆ(ξ) < ∞. Let us proceed to give a few examples. The Legendre transform of h(x) = 12 kxk2 is the function h itself. The Legendre transform of a linear function fi = hx, Pi i − 12 kPi k2 is < ∞ only when ξ = Pi . 4 D. SIERSMA AND M. VAN MANEN Let us now calculate the Legendre transform of f defined by (1). The domain where fˆ(ξ) < ∞ is the convex hull of the points {P1 , · · · , PN }. In fact, we have ! à X λi 1 f̂i (Pi ) = kPi k2 ⇒ fˆ(ξ) = inf kPi k2 2 2 i where the infimum is taken over all λi such that N X X λi Pi = ξ and ∀i : λi ≥ 0 and λi = 1. 1≤i≤N i=1 We can be more precise: for each ξ in the convex hull CH({P1 , · · · , PN }) there exist a subset α = α(ξ), of {P1 , · · · , PN } such that f (x) = fi (x) for exactly those i ∈ α. This requires x ∈ Vor(α), which should be non-void in this case. This means ξ ∈ Del(α). It follows à ! X λi fˆ(ξ) = kPi k2 2 i∈α and fˆ is an affine function on each of the cells in the Delaunay triangulation. Note that the domain of definition of maxi∈α fi is exactly Del(α). One can combine f and fˆ in the following function (ξ, x) 7→ F (x, ξ) = f (x) + fˆ(ξ) − hx, ξi. This function hides the Gateau differential. Namely, let h : Rn → R be a function. If the limit h(x + tξ) − h(x) h′ (x; ξ) = lim t↓0 t exists, it is called the Gateau differential. See [Hör94], theorem 2.1.22. For a convex set K ⊂ Rn the supporting function of K is the function ξ 7→ suphx, ξi x∈K Theorem 2.2.11 in [Hör94] says that the Gateau differential ξ → f ′ (x; ξ) is the supporting function of {µ | F (x, µ) = 0} . The Gateau differential is called Clarke’s generalized derivative in [APS97]. The set of which it is the support function is called ∂f (x). It is stated in that article that ∂f (x) is the convex hull of the gradients of the functions fi for which fi (x) = f (x). That last statement and the theorem 2.1.22 are all equivalent to what is neatly formulated in proposition 1 of [PR04]: Theorem 2. We have Del(α) = {ξ | F (x, ξ) = 0 ∀x ∈ Vor(α)}. Reversely Vor(α) = {x | F (x, ξ) = 0 ∀ξ ∈ Del(α)}. MEDIAL AXIS AND A MOMENT MAP 5 In other words, if x is a point in the relative interior of Vor(α) then ∂f (x) = Del(α), and reversely, if ξ is a point in the relative interior of Del(α) then ∂ fˆ(ξ) = Vor(α). We see that there are two ways of looking at the Delaunay triangulation. One simply through the duality, and the other through the generalized derivative of f . 2.3. Maslov dequantization and the moment map. In this section we approximate the function f with a parameterized smooth family of convex functions. The function f = max fi is convex, it can be approximated by convex functions x → f (h, x) = P fi (x) ), whose properties are similar to those of f . logh ( N i=1 h Note that f (x) = lim f (h, x) h→∞ due to Maslov dequantization, which is a fancy name for the identity: lim logh (ha + hb ) = max(a, b) (5) h→∞ that holds for any two real numbers a and b. Proposition 1. Let 0 < h′ < h, then • the functions f (x) and f (h, x) are convex, • ∀x ∈ Rn f (h, x) > f (h′ , x) > f (x), • the gradient of f (h, x) is a diffeomorphism from Rn to CH({P1 , · · · , PN }). Proof. We need to show that 1 1 1 f (h, x + y) ≤ (f (h, x) + f (h, y)), 2 2 2 or (6) N X fi ( 21 x+ 12 y) h )≤ i=1 i=1 1 à N X 1 hfi (x) ! 12 à N X hfi (x) i=1 1 ! 21 . 1 Put ai = h 2 fi (x) > 0 and bi = h 2 fi (y) > 0 then ai bi = hfi ( 2 x+ 2 y) , a2i = hfi (x) and b2i = hfi (y) . Hölder’s inequality says that ! 21 à N ! 12 à N N X X X b2i ai bi ≤ a2i i=1 i=1 i=1 which is equivalent to (6). The second part is proven using Jensen’s inequality, which asserts that for a convex function q: N N 1 X 1 X q( a1 ) ≤ q(ai ). N i=1 N i=1 We need to show that 0 < h′ < h ⇒ f (h′ , x) < f (h, x) 6 D. SIERSMA AND M. VAN MANEN We put q(y) = hlogh′ (y) = y logh′ (h) and ai = (h′ )fi (x) The function q(y) = y a , for some a > 1, thus it is convex. Use Jensen’s inequality and take logh at both sides to get the desired result. For the third statement consider the moment map from toric geometry: PN Pi efi (x) n n He : R → R He (x) = Pi=1 . N fi (x) i=1 e It is proved in [Ful93] that He is a real analytic diffeomorphism from Rn to the interior of CH({P1 , · · · , PN }), when that interior is not empty. One recognizes He as a gradient: !! à à N X = ∇x (f (e, x)) H(x) = ∇x log efi (x) i=1 Now consider the gradient of f (h, x): à N à !! P N X Pi hfi (x) Hh (x) = ∇x logh = Pi=1 hfi (x) . N fi (x) h i=1 i=1 The fi are of the form hx, Pi i − 21 kPi k2 . Put f˜i (u) = hu, Pi i − know that PN f˜i (u) i=1 Pi e u −→ H̃(u) = PN ˜ fi (u) i=1 e log(h) kPi k2 . 2 We already is a real analytic diffeomorphism from Rn to the interior of CH({P1 , · · · , PN }). Because we have Hh (x) = H̃(x log(h)) the proof is complete. ¤ When x does not lie on the Voronoi diagram there is an i, with 1 ≤ i ≤ N , fi (x) = f (x) and fi (x) > fj (x) for all j 6= i. It follows that limh→∞ Hh (x) = Pi . Note that limh→∞ Hh is no longer a diffeomorphism ; full Voronoi cells collapse to a point! We continue to describe the inverse of Hh . If t is smooth convex function that the Legendre transform t̂ of t is ¡ ¢ t̂(ξ) = h(∇x t)−1 (ξ), ξi − t (∇x t)−1 (ξ) As a consequence ∇t̂(ξ) = (∇x t)−1 (ξ) + (7) = ∂(∇x t)−1 ξj j=1 ∂ξj Pn − ¡ ¢ ∂(∇x t)−1 −1 t) (∇ t) (ξ) (∇ = x x j j=1 ∂ξ Pn (∇t)−1 (ξ). Thus, the inverse of Hh is the gradient of Legendre transform fˆ(h, x) of f (h, x). MEDIAL AXIS AND A MOMENT MAP 7 3. The medial axis and Maslov dequantization We have seen in the above that the Voronoi diagram can also be calculated using the affine functions fi , instead of the distance function. We have also seen that the corner locus of the function f defined in (1) is the Voronoi diagram. We will now follow the same procedure with the embedding γ : M → Rn . We will show next that the medial axis is a corner locus of a convex function. Let us prove theorem 1. Proof. We switch back from the function f as defined in (4) (the finite case) , to the function f defined in equation (2) (the continuous case). The function f (x) is convex, because it is the supremum of a number of convex functions. Let us prove part 2: the corner locus of f is the medial axis Cut(M ) of M . Because M is compact the supremum is attained for one or more points s = s(x) in M . We see that x ∈ Rn and s ∈ M , for which f (x) = fs (x), are related by hx − γ(s), γ̇(s)i = 0. (8) Here γ̇(s) is the tangent vector to γ and the inner product is taken as vectors with source in γ(s). So x lies on the normal of one or more s ∈ M . Of those that satisfy (8) we choose one that has the highest value of fsi (x), and thus the lowest value of 1 1 kxk2 − fsi (x) = kx − γ(si )k2 2 2 The s we get from sups∈M fs corresponds to the point closest to x on the submanifold M , because it is the s we get from the minimal distance 1 x → inf kx − γ(s)k2 . s∈M 2 If x does not lie on the medial axis there is a unique closest point s on M and it follows that f is smooth in some neighborhood of x. If x lies on the medial axis there are either two or more points that realize the supremum f (x) or x lies on the caustic. In both cases f is not smooth. Hence the points where f (x) is not differentiable form exactly the medial axis of M . Next we prove part 1: the integral formula for f in (2). The core of the argument is: approximate the integral on the right hand side by a point set on M. In the Riemann sum only the highest value of fs (x) counts as h → ∞. More precisely, let ι : Σ → M be a triangulation of M whose top-dimensional simplices are Γ (assume all of the same volume). Compute the lower and upper Riemann sums ! à ! à ¶ µZ X X sup vol(Γ)hfs (x) . hfs (x) ds ≤ logh inf vol(Γ)hfs (x) ≤ logh logh Γ∈Σ s∈ι(Γ) M Γ∈Σ s∈ι(Γ) Evaluate the lower Riemann sum ! à X ¢ ¡P fs (x) = inf vol(Γ)hfs (x) = logh ((#Γ) vol(Γ)) + logh logh Γ∈Σ inf s∈ι(Γ) h Γ∈Σ (9) s∈ι(Γ) = logh (vol(M )) + logh ¡P ¢ inf s∈ι(Γ) fs (x) h . Γ∈Σ 8 D. SIERSMA AND M. VAN MANEN From the Maslov dequantization identity (5) we know that the last expression tends to as h → ∞ (10) max( inf fs (x)). Γ∈Σ s∈ι(Γ) The upper Riemann sum tends to (11) max( sup fs (x)). Γ∈Σ s∈ι(Γ) Sub-dividing the simplices Γ further both (10) and (11) tend to sups∈M fs (x), which is what we needed to prove. We will now prove part 5, namely that the functions rh (x) are convex and monotonic in h: rh′ (x) ≤ rh (x) if h′ < h. The proof is completely analogous to the proof of proposition 1. To prove the convexity of rh (x) put 1 1 p(s) = h 2 fs (x) and q(s) = h 2 fs (y) and apply the continuous version of the Hölder inequality. The monotonicity of the sequence rh (x) follows from the continuous version of Jensen’s inequality, just as the monotonicity of the sequence f (h, x) in the proof of proposition 1 followed from the discrete version of Jensen’s inequality. Next we show part 6 that the differential of rh is diffeomorphism from Rn to the interior of CH(γ(M )). The gradient of rh is R γ(s)hfs (x) ds . (12) ∇rh (x) = MR hfs (x) ds M Consider again a triangulation ι : Σ → M of M , where each top-dimensional simplex Γ has the same volume. In each simplex Γ take a point PΓ . Consider the quotient P hx,PΓ i− 12 kPΓ k2 Γ∈Σ vol(Γ)PΓ h (13) . P hx,PΓ i− 21 kPΓ k2 Γ∈Σ vol(Γ)h From the remarks on Riemann sums above it follows that taking barycentric subdivisions the quotient (13) tends to (12). From proposition 1 we know that the image of (13) is the convex hull CH({PΓ | Γ ∈ Σ}). Hence the image of ∇rh (x) is CH(γ(M )). Part 3. We show that the derivative of f outside the corner locus of f maps x ∈ n R \ Cut(M ) to the nearest point on M . If x lies outside the medial axis then there is a point s0 on x such f (x) = fs0 (x). The γ(s0 ) is the unique closest to x on M . For all Rn ∋ y 6= x we have f (y) ≥ fs0 (x), so that certainly γ(s0 ) ∈ ∂f (x). We are done if we show that f is smooth at x. Again, because x does not lie on the medial axis there is a neighborhood U of x in Rn \ Cut(M ) there is a smooth map h : U → M that assigns to y the closest point h(y) on M . Hence, on U the f can be written f (y) = fh(y) (y), which is a smooth function. Before we prove part 4 we discuss the situation at the medial axis. At the medial axis the gradient of f is not defined. The generalized derivative can be calculated for points on MEDIAL AXIS AND A MOMENT MAP 9 the medial axis. For each point x on the medial axis there is a subset A(x) ⊂ M defined by A(x) = {s ∈ M | f (x) = fs (x)} p The point x is the center of a sphere of radius kxk2 − 2f (x) that is tangent to γ(M ) at all points of A(x). The generalized derivative ∂f (x) is CH(γ(A(x))). One can think of it as the Delaunay cell dual to a point on the medial axis. The union of all the A(x) for all x on medial axis Cut(M ) is M itself. The sets ∂f (x) fill CH(γ(M )). For each ξ ∈ CH(γ(M )) \ γ(M ) there is a unique x ∈ Cut(M ). Part 4: The derivative of the Legendre transform of f maps CH(γ(M )) \ γ(M ) onto the medial axis. The convex functions f and fˆ fulfil the duality relation: F (x, ξ) = f (x) + fˆ(ξ) − hx, ξi = 0. Standard theory of convex functions tells us that this is equivalent with each of the conditions: ∂f (x) = ξ , ∂ fˆ(ξ) = x. As we have seen above ∂f (x) = CH(γ(A(x))). And x is unique for each ξ. So ∂ fˆ(ξ) contains a single point. Now we are ready. ¤ The union of all the A(x) for all x on the interior medial axis IntCut(M ) is M itself. The convex structure of the Legendre transform gives is very intriguing, since it shows the similariry with the Delaunay triangulation. The cells dual to points on the medial axis can be seen if we take a lot of points on M and calculate the Delaunay triangulation of this point set. We get figure 2. This is another way of saying that the dual of the Figure 2. Delaunay triangulation for a lot of points on the curve γ from figure 1. interior medial axis in the sense of the theorem of Passare and Rullgård is the region in Rn bounded by manifold CH(M ) itself. A direct computation of the Legendre transform fˆ of f shows, that its domain is CH(M ) and that it is affine on each of the simplices ∂f (x) = CH(γ(A(x))). For a finite set A(x) P there is the formula fˆ(ξ) = A(x) 12 λi ||y||2 , where λi are coordinates in the hull. 10 D. SIERSMA AND M. VAN MANEN References [ABK98] [APS97] [Ful93] [Hör94] [Man39] [PR04] Nina Amenta, Marshall Bern, and Manolis Kamvysselis, A new voronoi-based surface reconstruction algorithm, SIGGRAPH ’98: Proceedings of the 25th annual conference on Computer graphics and interactive techniques (New York, NY, USA), ACM Press, 1998, pp. 415–421. A. A. Agrachev, D. Pallaschke, and S. Scholtes, On Morse theory for piecewise smooth functions, J. Dynam. Control Systems 3 (1997), no. 4, 449–469. MR MR1481622 (98m:49034) William Fulton, Introduction to toric varieties, Annals of Mathematics Studies, vol. 131, Princeton University Press, Princeton, NJ, 1993, The William H. Roever Lectures in Geometry. MR MR1234037 (94g:14028) Lars Hörmander, Notions of convexity, Progress in Mathematics, vol. 127, Birkhäuser Boston Inc., Boston, MA, 1994. MR MR1301332 (95k:00002) Szolem Mandelbrojt, Sur les fonctions convexes, C. R. Acad. Sci. Paris 209 (1939), 977–978. MR MR0000846 (1,137b) Mikael Passare and Hans Rullgård, Amoebas, Monge-Ampère measures, and triangulations of the Newton polytope, Duke Math. J. 121 (2004), no. 3, 481–507. MR MR2040284 (2005a:32005) Mathematisch Instituut, Universiteit Utrecht, PO Box 80010, 3508 TA Utrecht The Netherlands. E-mail address: [email protected] Department of Mathematics, Hokkaido University Kita 10, Nishi 8, Kita-Ku, Sapporo, Hokkaido, 060-0810, Japan E-mail address: [email protected]