Numerical computation of the cut locus via a variational approximation of the distance function

We propose a new method for the numerical computation of the cut locus of a compact submanifold of R 3 without boundary. This method is based on a convex variational problem with conic constraints, with proven convergence. We illustrate the versatility of our approach by the approximation of Voronoi cells on embedded surfaces of R 3 .

using the geodesic equation and deduced an approximation of the cut locus from there. In [14], the authors used the deformable simplicial complexes (DSC) method and finite differences techniques for geodesic computations to compute geodesic circles of increasing radius and their self-intersections, i.e. the cut locus. They applied the method to genus 1 surfaces. These papers contain no proof of convergence of the computed cut locus. Exact geodesic computation on discretized surfaces. This approach was used in [9,13]. In [13], the authors computed the geodesics on a convex triangulated surface. They deduced an approximation of the cut locus of the triangulated surface and filtered it according to the angle formed by the geodesics meeting at a point of the approximated cut locus, to make their approximation stable. They applied the method to ellipsoids. There is no proof of convergence. In [9], the authors computed shortest curves on a graph obtained from a sufficiently dense sample of points of the surface. From there they deduced an approximation of the cut locus and filtered it according to the maximal distance (called spread ) between the geodesics meeting at a point of the approximated cut locus. They proved that the set they compute converges to the cut locus (see [9], Thm. 4.1).
We may also mention [4], where the authors used some more geometric tools to compute (numerically) the cut locus of an ellipsoid or a sphere with some particular metric with singularities. Our method. A natural approach to approximate a cut locus would be to use a fast marching method which provides an efficient way to compute distance functions on a manifold. Unfortunately, classical algorithms do not ensure any convergence result related to the gradient of the approximation. We believe that this absence of estimate makes the approximation of the cut locus by these algorithms difficult to prove. In this article we introduce a new regularized approach designed to fill this gap and to obtain a reliable localization of the cut locus. Given a large constant > 0, let ∈ 1 ( ) be the minimizer of the following variational problem where ∇ denotes the gradient operator on the surface . Intuitively, is a mollification away from of the distance function to the point on . For > 0 to be chosen small, we will use the set as an approximation of Cut ( ). See Figure 1 for an illustration of the sensitivity of and , with respect to parameters and . This is justified by some theoretical results obtained in [12], which will be summarized in Sections 3 and 4. The Sections 2-4 are devoted to explaining how we arrived at such a set , . For now, let us give a bit of intuition about the different terms appearing in , . When perturbing the surface , we expect the same kind of instabilities as the ones observed in [3] in the case of the medial axis. Thus, two kinds of new points may appear in Cut ( ): (1) points where some minimizing geodesics meet with an angle close to zero, (2) points that are near the base point .
Hence, to make Cut ( ) more stable (and so more computable), we need to select points that are not too close to and such that some minimizing geodesics meet with an angle significantly larger than 0. Intuitively, having |∇ ( )| 2 ≤ 1 − 2 for some constant > 0 ensures that we are selecting points where minimizing geodesics meet with an angle significantly larger than 0, and replacing 2 by 2 / 2 ensures that we are selecting points that are away from . Some other definitions of , would have been possible (for instance, the squares are not needed), and this form has been chosen as it corresponds to the -medial axis introduced in [7] (see Sect. 2).
The rest of the paper is organized as follows. In Section 2, we recall the notion of -medial axis that was introduced in [7] and summarize some of its properties. In Section 3, following the strategy of the -medial axis, Figure 1. Approximation of the sets , (green color) for = 10 (first row) and = 50 (second row) for three different values of . The first three columns correspond to values of = 0.6, 0.2 and 0.06 respectively. The last column represents the associated solutions .
we define a " -cut locus" Cut ( ) and show that it can be used as an approximation of the complete cut locus for small enough. In Section 4, we recall the result from [12] which states that the set , defined above is a good approximation of Cut ( ) if is big enough. In Section 5, we discretize problem (1.1) using finite elements, to find a discrete minimizer ,ℎ , where ℎ > 0 is the step of the discretization. From this discrete minimizer ,ℎ , we obtain a function ,ℎ on , and we show that the set , ,ℎ := is a good approximation of , as ℎ → 0. In Section 6, we present the results of some numerical experiments.

The -medial axis
In this section, we recall briefly the notion of -medial axis introduced by Chazal and Lieutier [7]. Given an open subset Ω of R 2 , its medial axis ℳ(Ω) is defined as the set of points of Ω that have at least two closest points on the boundary Ω of Ω: where for any ∈ Ω, Ω ( ) is the distance from to Ω: The medial axis ℳ(Ω) is unstable with respect to small non-smooth perturbations of the boundary of Ω. To deal with this issue, in [7] Chazal and Lieutier defined the so called -medial axis of Ω by setting, for any > 0, It is further proved in Section 3, Theorem 2 of [7] that ℳ (Ω) has the same homotopy type as ℳ(Ω), for small enough. These facts justify that ℳ (Ω) is a good approximation of ℳ(Ω), for small enough. The crucial difference though is that ℳ (Ω) is stable with respect to small variations of the boundary of Ω, whereas ℳ(Ω) is not. We refer the reader to Section 4 of [7] for precise statements and proofs.
To motivate the next section, we will also use an alternative definition of the -medial axis. Given a point ∈ Ω, let Θ( ) be the center of the smallest ball containing all the closest points to on Ω. In Section 2.1 of [7], a vector field ∇ Ω (originally denoted only by ∇) is defined on Ω by: This vector field coincides with the classical gradient of Ω wherever Ω is differentiable, so it can be thought of as a generalized gradient of Ω . Moreover, we have the following relation (see Eq. (1) in Sect. 2.1 of [7]): Therefore, we have the following equivalent definition of the -medial axis:

-Cut locus
We want to define a set similar to the -medial axis in the case of the cut locus Cut ( ). To this end, we need a notion of generalized gradient for the distance function . The notion of generalized gradient we use is presented in Section 1.3 of [16] in the context of Alexandrov spaces and in Section 2.2 of [12]. Here, we introduce the notion omitting the short proofs otherwise needed. First note that, as stated in Proposition 2.7 of [12], the function is locally semiconcave on ∖ { }, which means that for any unit speed geodesic : [0, 1] → ∖ { }, there exists a constant > 0 such that the function ↦ → 2 − ( ( )) is convex on [0, 1]. From there, one gets that for any point ∈ ∖ { } and any direction ∈ , admits a directional derivative where exp denotes the Riemannian exponential map at . Furthermore, for any ∈ ∖ { }, the map ↦ → + ( ) admits a unique maximizer on the closed ball (0, 1) ⊂ . The generalized gradient of is then defined as What is more, we have the following formula See Lemma 3.3 for a geometric interpretation of the generalized gradient.
Analogously to (2.2), for > 0, we define the -cut locus as We have the following proposition from Proposition 2.9 of [12].
Proposition 3.1. The map ↦ → Cut ( ) is nonincreasing, and In addition, the following proposition holds. We recall that Cut ( ) is always connected (see [15] for instance).
If is a real analytic surface, then for > 0 small enough, one of the connected components of Cut ( ) has the same homotopy type as Cut ( ), while the other connected components, if any, are contractible.
These two propositions justify that Cut ( ) is a good approximation of Cut ( ), for > 0 small enough. Before proving Proposition 3.2, we prove the following lemma. Proof. For = 1, 2, let us set = −˙( ( )). Let us denote by exp the Riemannian exponential map at the point . Let 0 ∈ (0, ( )) and = exp ( 0 ). Note that we have / ∈ Cut ( ), so the function is smooth at , and its gradient is − . Given ∈ such that | | = 1, using ( ) = ( ) + ( ), we have Given that the angle between 1 and 2 is , there exists ∈ {1, 2}, such that the angle between and is at most − /2. Thus the last inequality gives + ( ) ≤ cos( /2). This concludes the proof.
Remark 3.4. In the previous lemma, in the case where there are exactly two minimizing geodesics arriving at , one can show that the inequality is an equality.
Using Lemma 3.3, Proposition 3.2 will mainly be a consequence of Proposition 3.4 from [9], which is recalled in Proposition 3.6 below. Following [9], we will use the following terminology. Let be a finite connected graph embedded in . A point of a finite connected graph is called a tree point if is a leaf of or ∖ { } has a connected component whose closure is a tree. Otherwise, is called a cycle point. As shown in the proof of Proposition 3.5 from [9], the closure of the set of cycle points of a finite connected graph is a deformation retract of , hence it is connected. We will also use the following lemma. Proof. Let 1 , . . . , be the connected components of ∖ . As contains all cycle points of , for any 1 ≤ ≤ , is a tree and ∩ is a singleton { }. By contracting all to their roots , we obtain that is a deformation retract of .
Let ∈ be such that there exist two minimizing unit speed geodesics 1 and 2 from to . Following [9], the spread between 1 and 2 is defined as We recall that, as is real analytic, the cut locus Cut ( ) is a finite graph (see [15] in dimension 2 and [5] for the generalization to arbitrary dimensions). In [9], the authors proved the following: If the spread of any two minimizing unit speed geodesics 1 and 2 from to is smaller than the injectivity radius of , then is a tree point of Cut ( ).
Proof of Proposition 3.2. According to Lemma 3.3, given any > 0, if has been taken small enough, then for any point ∈ Cut ( ) ∖ Cut ( ), the angle between any two minimizing unit speed geodesics 1 and 2 from to is smaller than at . As geodesics verify a second order differential equation, if their angle at is small, then their spread is also small. Therefore, applying Proposition 3.6, we deduce that if has been taken small enough, then any point ∈ Cut ( ) ∖ Cut ( ) is a tree point of Cut ( ). Stated otherwise, Cut ( ) contains all cycle points of Cut ( ). Moreover, Cut ( ) is a closed set. Indeed, this is a consequence of the semiconcavity of and the lower semicontinuity of the norm of the gradient of semiconcave functions (see [12], Prop. 7.2). Thus, Cut ( ) contains the closure of the cycle points of Cut ( ), which is connected. In particular, there exists a connected component of Cut ( ) that contains the set of the cycle points of Cut ( ). By Lemma 3.5, is a deformation retract of Cut ( ). This completes the proof.
Therefore, we will use Cut ( ) as an approximation of Cut ( ) for small enough.

Approximation with a variational problem
For > 0, recall that is the minimizer in (1.1). For > 0, let us define the set , by We have the following theorem (see [12], Thms. 1.1 and 1.3): Theorem 4.1. There exists 0 > 0 such that for any > 0 , the function is locally 1,1 on ∖ { }, and = in a neighborhood of . For any > ′ > 0 , Therefore, we can use , as an approximation of Cut ( ). All in all, we will use , as an approximation of Cut ( ).
In the following, we will always assume that we have > 0 .

Finite elements of order on a surface approximation of order
In this section we introduce a discretization framework adapted to the variational problem (1.1), based on finite elements. We follow the notations of [8,11].
Let be a compact oriented smooth two-dimensional surface embedded in R 3 . For ∈ , we denote by ( ) the oriented normal vector field on . Let : R 3 → R be the signed distance function to the surface and = { ∈ R 3 , | ( )| < } the tubular neighborhood of of width > 0. It is well known that if is small enough (for instance 0 < < min =1,2 We consider ℎ a triangular approximation of whose vertices lie on and whose faces are quasi-uniform and shape regular of diameter at most ℎ > 0. Moreover, we will assume that ℎ , the set of triangular faces of ℎ , are contained in some tubular neighborhood such that the map defined by (5.1) is unique. For ≥ 1 and for a triangle ∈ ℎ , we consider the Lagrange basis functions Φ 1 , . . . Φ of degree and define the discrete projection on ℎ by: where 1 , . . . , are the nodal points associated to the basis functions. Now we can define ℎ a polynomial approximation of order of associated to ℎ ℎ = { ( ), ∈ ℎ }. Observe that by definition the image by of the nodal points are both on and on ℎ . Let us now introduce the finite element spaces on ℎ = 1 ℎ and ℎ for ≥ 2. For every integer ≥ 1, let where P is the family of polynomials of degree at most . Analogously, for ≥ 2, let , ℎ = {︀ˆ∈ 0 ( ℎ ),ˆ= ∘ −1 , for some ∈ ℎ }︀ . (5.5) Analogously to (1.1), we will consider the following discrete variational problem: and some fixed nodal point of the mesh ℎ .

Convergence of the lifted minimizers
In order to prove the convergence of our numerical approach, let us first establish that our discrete problem converges in values in the sense of Proposition 5.2. For a function defined on ℎ , we introduce its lifted function defined on , by the relation ( ( )) = ( ). We focus our analysis on the piecewise linear case = = 1. We will use the notation ℎ := 1 ℎ and ℎ := 1 ℎ . For every ℎ > 0, the convex optimization problem (5.6) has a unique solution ,ℎ .
Lemma 5.1. The differential of the projection : → R 3 onto , when restricted to the tangent space of ℎ , is the identity, up to order 1 in ℎ: Proof. The identity estimate on is a direct consequence of [11] equations (4.12), (4.13), (4.10) and (4.11).
In the proof of the next proposition, we will need the following lemma.
We can now establish the convergence of the minimizers: Then, is admissible for problem (1.1), so ( ) ≥ ( ). Moreover, the following algebraic identity holds Therefore, we have 1 2 which proves, together with Proposition 5.2, that In particular, as is compact, the gradient of the function ,ℎ − also goes to 0 in the 1 ( ) norm as ℎ → 0. Recall that the functions (︁ In the proof of Proposition 5.2, we showed that ,ℎ = ,ℎ (1 + (ℎ)) (Eq. (5.14)). Together with (5.18) and (5.19), this concludes the proof. Remark 5.6. The function is generically suspected to be no more than 1,1 regular. This important difficulty makes the extension of the previous convergence analysis to higher order elements not straightforward. However, based on our numerical experiments (see the following paragraph), we believe that a ( , ) approximation improves the order of convergence without being able to provide a complete analysis of this order due to this lack of regularity. Sections 3-5 together justify the claim that the set , ,ℎ is a good approximation of the cut locus of in , if is big enough and and ℎ are small enough. We tried to estimate numerically the order of convergence of our approximation. Unfortunately, the precision of our simulations was not enough to determine if the expected order of convergence 1/4 is optimal or not.

Cut locus approximation
We established the convergence of the minimizers of problem (5.6) when ℎ tends to 0. For a fixed ℎ > 0, this convex discrete problem is of quadratic type with an infinite number of conic pointwise constraints. By the way, it is important to observe that for = = 1, the gradient pointwise bounds for a function of P 1 is equivalent to a single discrete conic constraint on every triangle with respect to the degrees of freedom of P 1 ( ℎ ). In this simplified context, we observed in our experiments that using P 1 elements may lead to approximated cut loci with some tiny artificial connected components (see Fig. 2). Observe that these artifacts do not contradict our convergence estimates in Proposition 5.5. Motivated by this lack of precision, we use in all following illustrations elements of order > 1.
For the general case > 1, the bound constraint on the gradient cannot be easily reduced to a finite set of discrete constraints. In our computations, we approximated the constraint |∇ ℎ | ∞ ( ℎ ) ≤ 1 by forcing the inequality only on a finite number of points of the mesh. In practice, we imposed these constraints on Gauss quadrature points on every triangle of ℎ .
We illustrate in Figures 3-6 the approximation of the cut locus provided by our approach. These computations have been carried out on meshes of approximately 10 5 triangles for = 2 and = 3 using a high precision quadrature formula associated to 17 Gauss points on every element of the mesh. Moreover, for = 3, we imposed the conic gradient constraints on the = 9 Gauss points of every triangle. In order to solve the resulting linear conic constrained quadratic optimization problem, we used the JuMP modeling language and the finite elements library Getfem++ [10,17] combined with the MOSEK optimization solver [2]. For such a precision, the optimization solver identified a solution in less than one hour on a standard computer.
Observe that our approximations of cut loci provide sets with a number of handles equal to twice the genus of the supporting surfaces. This fact is in agreement with Proposition 3.2 since the cut locus has he same homotopy group as the surface (see [6] for instance).

Approximation of the boundary of Voronoi cells
All previous theoretical results still hold if we replace the source point by any compact subset of the surface . For instance, if is replaced by a set of points, the singular set of the distance function can be decomposed as the union of the boundary of Voronoi cells and the cut loci of every point intersected with its Voronoi cell. As a consequence, if the distribution of source points is homogeneous enough, that is every Voronoi cell is small enough, the singular part of the distance function will be exactly equal to the boundary of the Voronoi cells. We illustrate this remark in the following experiments. We used exactly the same framework as in previous sections and just replaced the pointwise condition at with the analogous pointwise Dirichlet conditions at every source  point. Figures 7 and 8 represent the Voronoi diagrams obtained with 10, 30 and 100 points for surfaces of genus 2 and 3. The computational complexity is exactly of the same order as with a single source point.