Voronoi Cells in Metric Algebraic Geometry of Plane Curves
Madeline Brandt Madeleine Weinstein


 
Voronoi cells of varieties encode many features of their metric geometry. 
We prove that each Voronoi or Delaunay cell of a plane curve appears as the limit of a sequence of cells obtained from point samples of the curve.
We use this result to study metric features of plane curves, including the medial axis, curvature, evolute, bottlenecks, and reach. In each case, we provide algebraic equations defining the object and, where possible, give formulas for the degrees of these algebraic varieties. We show how to identify the desired metric feature from Voronoi or Delaunay cells, and therefore how to approximate it by a finite point sample from the variety.



Designer fermion models in functionalized graphene bilayers
    J. Fernández-Rossier [On leave from Departamento de Física Aplicada, Universidad de Alicante,  Spain] 
    December 30, 2023
==========================================================================================================












§ INTRODUCTION





Metric algebraic geometry addresses questions about real algebraic varieties involving distances. For example, given a point x on a real algebraic plane curve X ⊂ℝ^2, we may ask for the locus of points which are closer to x than to any other point of X. This is called the Voronoi cell of X at x <cit.>.
The boundary of a Voronoi cell consists of points which have more than one nearest point to X. So we may ask, given a point in ℝ^2, how close must it be to X in order to have a unique nearest point on X? This quantity is called the reach, and was first defined in <cit.>. 


We use Voronoi cells and their duals, Delaunay cells (see Definition <ref>), to study metric features of plane curves. The following theorem makes precise the idea behind Figures <ref> and <ref>. 


Let X be a compact algebraic curve in ℝ^2 and {A_ϵ}_ϵ↘ 0 be a sequence of finite subsets of X containing all singular points of X such that every point of X is within distance ϵ of some point in A_ϵ.

    
  * Every Voronoi cell is the Wijsman limit (see Definition <ref>) of a sequence of Voronoi cells of {A_ϵ}_ϵ↘ 0.
    
  * If X is not tangent to any circle in four or more points, then every maximal Delaunay cell is the Hausdorff limit (see Definition <ref>) of a sequence of Delaunay cells of {A_ϵ}_ϵ↘ 0.









Voronoi diagrams of finite point sets are widely studied and have seen applications across science and technology, most notably in the natural sciences, health, engineering, informatics, civics and city planning. For example in Victoria, a state in Australia, students are typically assigned to the school to which they live closest. Thus, the catchment zones for schools are given by a Voronoi diagram <cit.>. 
Metric features of varieties, such as the medial axis and curvature of a point, can be detected from the Voronoi cells of points sampled densely from a variety.  Computational geometers frequently use Voronoi diagrams to approximate these features and reconstruct varieties <cit.>. 



The reach of an algebraic variety is an invariant that is important in applications of algebraic topology to data science. For example, the reach determines the number of sample points needed for the technique of persistent homology to accurately determine the homology of a variety <cit.>. For an algebraic geometric perspective on the reach, see <cit.>. The medial axis of a variety is the locus of points which have more than one nearest point on X. This gives the following definition of the reach. 




The reach τ(X) of an algebraic variety X ⊂ℝ^n is the infimum of the set of distances from any point on X to a point on the medial axis of X.


The paper <cit.> describes how the reach is the minimum of two quantities. 
We have 

    τ(X) = min{ q , ρ/2},

where q is the minimum radius of curvature (Definition <ref>) of points in X and ρ is the narrowest bottleneck distance (Definition <ref>). An example is depicted in Figure <ref>.




The paper is organized as follows. We begin with a systematic treatment in Section <ref> of convergence of Voronoi cells of increasingly dense point samples of a variety. Here, we also introduce Delaunay cells, which are dual to Voronoi cells but are compact, and thus able to exhibit Hausdorff convergence. This gives the proof of Theorem <ref>, split among Theorems <ref> and <ref>, as well as Proposition <ref>, which treats the singular case separately. Theorem <ref> is robust because it is not affected by the distribution of the point sample.
Theorem <ref> provides the theoretical foundations for estimating metric features of a variety from a point sample. We do this for the medial axis (Section <ref>), curvature and evolute (Section <ref>), bottlenecks (Section <ref>) and reach (Section <ref>). For each of these metric features, we first give defining equations and where possible a formula for the degree. We then turn our attention to detecting information about a real plane curve X from its Voronoi cells. For each metric feature, we state a theoretical result about how to detect the feature from the Voronoi cells of X or a subset of X. Corollaries to Theorem <ref> provide convergence results for these features. The overall aim is to give a path to compute the metric features of an algebraic plane curve X from Voronoi cells of dense point samples of X. We use the butterfly curve

    b(x,y) = x^4 - x^2 y^2 + y^4 - 4 x^2 - 2 y^2 - x - 4 y + 1

in our examples. 
 In computational geometry and data science, these problems are often considered when there is noise in the sample. In this paper we assume that our samples lie precisely on the curve X. 










§ VORONOI AND DELAUNAY CELLS OF VARIETIES AND THEIR LIMITS


Let X ⊂ℝ^n be a nonempty real algebraic variety, that is, the zero locus of a set of polynomial equations with real coefficients. We note that this definition allows a variety to be reducible. We call X a curve if it is of dimension 1 and say X is smooth if the real locus is smooth.  Let d(x,y) denote the Euclidean distance between two points x,y ∈ℝ^n.

The Voronoi cell of x ∈ X is

    Vor_X(x) = {y ∈ℝ^n | d(y,x) ≤ d(y,x')  for all x' ∈ X}.


An example of a Voronoi cell is given in Figure <ref>. 
This is a convex semialgebraic set whose dimension is equal to codim(X) so long as x is a smooth point of X. It is contained in the normal space to X at x:

    N_X(x) = {u ∈ℝ^n  | u-x is perpendicular to the tangent space of X at x}.

The topological boundary of the Voronoi cell Vor_X(x) consists of the points in ℝ^n that have two or more closest points in X, one of which is x. The collection of boundaries of Voronoi cells is described as follows.


The medial axis M(X) of an algebraic variety X ⊂ℝ^n is the collection of points in ℝ^n that have two or more closest points in X. An example of the medial axis is given in Figure <ref>. 

 
Let B(p,r) denote the open disc with center p ∈ℝ^n and radius r > 0. We say this disc is inscribed with respect to X if X ∩ B(p,r) = ∅ and we say it is maximally inscribed if no disc containing B(p,r) shares this property.
Each inscribed disc gives a Delaunay cell, defined as follows.


Given an inscribed disc B of an algebraic variety X ⊂ℝ^n, 
the Delaunay cell Del_X(B) 
is
conv(B∩ X). 
An example of a Delaunay cell and the corresponding maximally inscribed disc is given in Figure <ref>. 





For plane curves, the collection of centers of all inscribed spheres which give maximal Delaunay cells (Delaunay cells which are not contained in any other Delaunay cell) is the Euclidean closure of the medial axis. Points of an algebraic plane curve X which are themselves maximal Delaunay cells are points of X with locally maximal curvature. In this case, the maximally inscribed circle is an osculating circle, see Definition <ref>.






We now describe two convex sets whose face structures encode the Delaunay and Voronoi cells of X.
We embed ℝ^n in ℝ^n+1 by adding a coordinate. We usually imagine that this last coordinate points vertically upwards.
So, we say that x ∈ℝ^n+1 is below y ∈ℝ^n+1 if x_n+1≤ y_n+1 and all other coordinates are the same. 
 Let 

    U = { x ∈ℝ^n+1 | x_n+1 = x_1^2 + ⋯ + x_n^2 }

be the standard paraboloid in ℝ^n+1. 
If p ∈ℝ^n, then let p_U = (p,||p||^2) denote its lift to U. 


Given a convex set C ⊂ℝ^n+1, a convex subset F ⊂ C is called a face of C if for every x ∈ F and every y,z ∈ C such that x ∈ conv(y,z), we have that y,z ∈ F. We say that a face F is exposed if there exists an exposing hyperplane H such that C is contained in one closed half space of the hyperplane and such that F = C ∩ H. We call an exposed face F a  lower exposed face of C if there is an exposing hyperplane lying below C.



The Delaunay lift of an algebraic variety X ⊂ℝ^n is the convex set 

    P^*_X = conv(x_u | x ∈ X) + {(0,…,0,λ) : λ∈ℝ_≥ 0}⊂ℝ^n+1,

where we recall that x_u = (x , || x || ^2) and use + to denote the Minkowski sum.
 The Delaunay lift of the butterfly curve is shown in Figure <ref>.





We now study how the lower exposed faces of the Delaunay lift P^*_X project to conv(X), and give the Delaunay cells of X.




Let X ⊂ℝ^n be an algebraic variety. Let π: ℝ^n+1→ℝ^n be the projection onto the first n coordinates. A subset F ⊂ P^*_X is a lower exposed face if and only if π(F) is a Delaunay cell of X. Furthermore, if H_F is the hyperplane exposing F, then π(U∩ H_F) is an inscribed sphere of X and π(F) = Del_X(π(U∩ H_F)).



The map from ℝ^n →ℝ^n+1 defined by x ↦ x_U= (x,||x||^2) lifts every sphere in ℝ^n to the intersection of a hyperplane H with U <cit.>. Moreover, the projection of the intersection of any hyperplane H (satisfying π(H) = ℝ^n) with 
U gives a sphere in ℝ^n <cit.>.

Given a Delaunay cell Del_X(B) for some inscribed sphere B, we have that P^*_X lies above the corresponding hyperplane H. This is because any points below H would project to points in X lying inside of B, contradicting the condition that X ∩ B = ∅ for an inscribed disc B. So, H is the exposing hyperplane of the face (Del_X(B))_U.

Suppose F ⊂ P^*_X is a lower exposed face with exposing hyperplane H_F. The interior of the sphere π(U∩ H_F) contains no points of X, because if it did contain a point x, then x_U would lie in the lower half-space of H_F, which does not intersect P^*_X. Then π(U∩ H_F) is the inscribed disc corresponding to a Delaunay cell. 

Since π(U∩ H_F) is a sphere, we have Del_X(π(U∩ H_F)) = conv(π(U∩ H_F)∩ X). Let X_U denote the lift of X to ℝ^n+1. Then π(U∩ H_F)∩ X = π(U∩ H_F∩ X_U) = π(H_F ∩ X_U), and so Del_X(π(U∩ H_F)) = conv(π(H_F ∩ X_U)) = π(conv(H_F ∩ X_U)) = π(F).





We may define a convex set whose faces project down to the Voronoi cells as follows. For any point x ∈ X, let T(x) denote the hyperplane in ℝ^n+1 through x_U =  (x,||x||^2) tangent to the paraboloid U. Let T(x)^+ be the closed half-space consisting of all points in ℝ^n+1 lying above the hyperplane T(x).


The Voronoi lift of an algebraic variety X ⊂ℝ^n is the convex set
P_X = ∩_x ∈ X T(x)^+.
The Voronoi lift of the butterfly curve is shown in Figure <ref>. 





The lower exposed faces of the Voronoi lift P_X project to Voronoi cells of X, as we now show.



Let X ⊂ℝ^n be an algebraic variety. Let π: ℝ^n+1→ℝ^n be the projection onto the first n coordinates. A subset F of the Voronoi lift P_X is an exposed face of P_X if and only if π(F) is a Voronoi cell of X. Furthermore, if H_F is the hyperplane exposing F and U∩ H_F ≠∅,  then U∩ H_F is a point and 
π(F) = Vor_X(π(U ∩ H_F)).



For some point x ∈ X, consider P_X ∩ T(x). Let p ∈ℝ^n. There exists p' ∈ T(x) with π(p_U) = π(p'). The distance from  p_U to the point p' is the square of the distance d_ℝ^n(π(p),π(x)) <cit.>.
Therefore, P_X ∩ T(x) consists of those points p' for which the distance d_ℝ^n(p,x) is minimal over all x ∈ X. In other words, π(P_X ∩ T(x)) = Vor_X(x). 

Suppose F ⊂ P_X is an exposed face with exposing hyperplane H_F such that U∩ H_F ≠∅. Let p ∈ U∩ H_F. Since U ⊂ P_X we have that p ∈ P_X. Then, p ∈ F = H_F ∩ P_X. This implies H_F is the tangent hyperplane to U at the point p, so in particular, p=U ∩ H_F. Since p is on the boundary of P_x, we have π(p)∈ X and T(π(p)) = H_F. We have π(F)= π(P_X ∩ T(π(p))) =  Vor_X(π(p)) = Vor_X(π(U ∩ H_F)), where in the second equality we use the result in the preceding paragraph. 



There is a sense in which the Voronoi lift P_X and the Delaunay lift P^*_X are dual. We now describe this relationship. Suppose that X  is not contained in any proper linear subspace of ℝ^n. This implies that P_X is pointed, meaning it does not contain a line. Therefore, it is projectively equivalent to a compact set <cit.>. Embed ℝ^n+1 into ℙ^n+1 by the map 

    ι (x_1, …, x_n+1) = (1:x_1 : ⋯ : x_n+1).

Let l be the transformation of ℙ^n+1 defined by the following (n+2) × (n+2) matrix

    [  1  0  ⋯  0  1;  0  2  0  ⋯  0;  ⋮  0  ⋱  0  ⋮;  0  ⋯  0  2  0; -1  0  ⋯  0  1 ].

Then by <cit.> the projective transformation l maps U to the sphere S ⊂ℝ^n+1. The tangential hyperplane at the north pole (1: 0 : ⋯ : 0 : 1) is the image of the hyperplane at infinity. Moreover, the topological closure of l(P_X) is a compact convex body so long as the origin is in the interior of P^*_X. In this case, we call the convex body l(P_X) the Voronoi body. The Voronoi body is full dimensional and contains the origin in its interior. Its polar dual

    l(P_X)^∘ := { y ∈ℝ^n+1 : ∑_i=1^n x_i y_i ≤ 1 for all  x ∈π(P_X)  }
 
is also full dimensional and has the origin in its interior. If we apply l^-1 to l(P_X)^∘ we obtain an unbounded polyhedron, which is exactly the Delaunay lift P^*_X of X. For more details, see <cit.>. 


We now study convergence of Voronoi and Delaunay cells. More precisely, given a real algebraic curve X and a sequence of samplings A_N ⊂ X with |A_N| = N, we show that Voronoi (or Delaunay) cells from the Voronoi (or Delaunay) cells of the A_N limit to Voronoi (or Delaunay) cells of X. We begin by introducing two notions of convergence which describe the limits.

The Hausdorff distance of two compact sets B_1 and  B_2 in ℝ^n is defined as

    d_h(B_1,B_2) := sup{sup_x ∈ B_1inf_y ∈ B_2 d(x,y), sup_y ∈ B_2inf_x ∈ B_1 d(x,y) }.

More intuitively, we can define this distance as follows. If an adversary gets to put your ice cream on either set B_1 or B_2 with the goal of making you go as far as possible, and you get to pick your starting place in the opposite set, then d_h(B_1,B_2) is the farthest the adversary could make you walk in order for you to reach your ice cream.



A sequence {B_ν}_ν∈ℕ of compact sets is Hausdorff convergent to B if d_h(B,B_ν) → 0 as ν→∞.




Given a point x ∈ℝ^n and a closed set B ⊂ℝ^n, define

    d_w(x,B) = inf_b ∈ B d(x,b).

A sequence {B_ν}_ν∈ℕ of compact sets is  Wijsman convergent to B if for every x ∈ℝ^n, we have that 

    d_w(x,B_ν) → d_w(x,B).




An ϵ-approximation of a real algebraic variety X is a discrete subset A_ϵ⊂ X such that for all y∈ X there exists an x ∈ A_ϵ so that d(y,x) ≤ϵ.
By definition, when X is compact a sequence of ϵ-approximations as ϵ approaches 0 is Hausdorff convergent to X. For all X, a sequence of ϵ-approximations as ϵ approaches 0 is Wijsman convergent to X.
We use Wijsman convergence as a variation of Hausdorff convergence which is well suited for unbounded sets. Delaunay cells are always compact, while Voronoi cells may be unbounded. 


We now study convergence of Delaunay cells of X, and introduce a condition on real algebraic varieties which ensures that the Delaunay cells are simplices.



We say that an algebraic variety X⊂ℝ^n is Delaunay-generic if X does not meet the closure of any maximally inscribed disc at greater than n+1 points.



 The standard paraboloid U in any dimension n+2 is not Delaunay-generic because it contains n-spheres. 


Although the focus of this paper is on algebraic curves in ℝ^2, we state the following theorem for curves in ℝ^n because the proof holds at this level of generality.



Let X ⊂ℝ^n be a Delaunay-generic compact algebraic curve, and let {A_ϵ}_ϵ↘ 0 be a sequence of ϵ-approximations of X. Every maximal Delaunay cell is the Hausdorff limit of a sequence of Delaunay cells of A_ϵ. 


Consider a sequence {A_ϵ}_ϵ↘ 0 of ϵ-approximations of X, where ϵ↘ 0 indicates a decreasing sequence of positive real numbers ϵ_ν for ν∈ℕ. We will study the convex sets
P^*_A_ϵ = conv(a_U | a ∈ A_ϵ),
where a ↦ a_U = (a,||a||^2) lifts a to the paraboloid U. 
The lower faces of P^*_A_ϵ project to Delaunay cells of A_ϵ <cit.>.

We now apply <cit.> to our situation.
This result says the following.
Let C be a curve and B_ϵ be a sequence of ϵ-approximations of C.
Suppose every point on C which is contained in the boundary of conv(C) is an extremal point of conv(C), meaning it is not contained in the open line segment joining any two points of conv(C). Let F be a simplicial face of conv(C) which is an exposed face of conv(C) with a unique exposing hyperplane.
Then F is the Hausdorff limit of a sequence of facets of conv(B_ϵ). We apply this result in the case when C= X_U = {x_U ∈ℝ^n+1 | x ∈ X} and B_ϵ = (A_ϵ)_U = {a_U ∈ℝ^n+1 | a ∈ A_ϵ}. 

Since every point on U is extremal in conv(U) and conv(X_U) ⊂ conv(U), every point on X_U which is contained in the boundary of conv(X_U) is also extremal in conv(X_U). 
A maximal Delaunay cell of X is a simplex because X is Delaunay-generic.
Consider a maximal Delaunay cell of X which is not a vertex. It has a unique description as Del_X(B) for a disc B.
Proposition <ref> establishes a one-to-one correspondence between such Delaunay cells and lower exposed faces of P^*_X, which are uniquely exposed by the hyperplane containing (∂B)_U.

 In this case, <cit.> holds, so the result is proved.


If a maximal Delaunay cell is a vertex, then it is a point x ∈ X. It is then also an extremal point of conv(X_U). Since conv(B_ϵ) is sequence of compact convex sets converging in the Hausdorff sense to P_X^*, by <cit.> there exists a sequence of points of B_ϵ which are extremal points of conv(B_ϵ) converging to x_U. So, their projections are Delaunay cells of A_ϵ converging to x, since every point in a finite point set is a Delaunay cell of that point set. 




We will now study limits of Voronoi cells, using results from <cit.>, which studies convergence of Voronoi cells of r-nice sets (for a definition, see <cit.>). In the plane, these are open sets whose boundary satisfies some properties. In particular, open sets whose boundaries are an algebraic curve with positive reach r satisfy the r-nice condition. All closed C^2 submanifolds of ℝ^n have positive reach, and thus in particular, a compact smooth algebraic curve in ℝ^2 has positive reach; we refer the reader to <cit.> for further discussion of which sets have positive reach. 

To study continuity and convergence of closed sets in the plane, Brandt uses the hit-miss topology on the set ℱ of closed subsets of the plane <cit.>. 


In the hit-miss topology, a sequence {F_n} converges to F if and only if

    
  *  for any p ∈ F, there is a sequence p_n ∈ F_n such that p_n → p; and
    
  * 
    if there exists a subsequence p_n_k∈ F_n_k  converging to a point p, then p ∈ F.

Then, to determine if a function with range in ℱ is continuous, we need to examine the above conditions for sequences of sets obtained by applying the function to countable convergent sequences in ℝ^2. If all such sequences satisfy (<ref>) then the function is upper-semicontinuous. If all such sequences satisfy (<ref>) then it is lower-semicontinuous. If a function satisfies both then it is continuous.




Let X ⊂ℝ^2 be a smooth algebraic plane curve. Then the function Vor_X: X →ℱ sending x ↦ Vor_X(x) is continuous in the hit-miss topology.


By <cit.>, the Voronoi function Vor_X:X →ℱ is lower semicontinuous. By <cit.>, if the curve is C^2 and the skeleton (locus of centers of maximally inscribed discs) is closed, then the Voronoi function is continuous. A smooth algebraic curve is C^2. The skeleton is closed because a smooth curve satisfies the r-nice condition, and r-nice curves have closed skeletons <cit.>.


By <cit.>, convergence in the hit-miss topology is equivalent to Wijsman convergence. In what follows, we rephrase the results from <cit.> in the setting of Wijsman convergence of Voronoi cells of plane curves, and extend them to singular curves. 




Let X be a compact smooth algebraic curve in ℝ^2 and {A_ϵ}_ϵ↘ 0 be a sequence of ϵ-approximations of X. Every Voronoi cell is the Wijsman limit of a sequence of Voronoi cells of A_ϵ.



By Lemma <ref>, the function Vor_X: X →ℱ is continuous. Theorem 3.1 from <cit.> states that in this case, if x_ϵ is a sequence such that x_ϵ∈ A_ϵ and x_ϵ→ x, then Vor_A_ϵ(x_ϵ) → Vor_X(x). Such a sequence must exist because for all y ∈ X, there exists a y_ϵ∈ A_ϵ such that d(y,y_ϵ) ≤ϵ.



We now investigate the structure of Voronoi cells of different types of singular points. We show examples in which the Voronoi cell at a singular point is 0-, 1-, and 2- dimensional in Proposition <ref> and Figure <ref>. First, we need a glueing lemma. 



Let C and D be subsets of ℝ^2 containing a point p ∈ C ∩ D. Then 

    Vor_C ∪ D(p) = Vor_C(p) ∩ Vor_D(p).



A point x ∈ Vor_C ∪ D(p) is closer to p than it is to any other point of C or D. On the other hand, a point in Vor_C(p) ∩ Vor_D(p) is closer to p than it is to any other point of C or D.








Let X ⊂ℝ^2 be a real algebraic plane curve and p be a singular point.

    
  * If p is a node, then its Voronoi cell is 0-dimensional and equal to p;
    
  * If p is a tacnode, then its Voronoi cell is 1-dimensional.
    
  * If p is an isolated point, its Voronoi cell is 2-dimensional;







  * If p is a node, then we claim that the only point contained in Vor_X(p) is p. At p, the curve meets in two branches which have distinct tangent directions at p. If we treat this as two separate 1-dimensional subsets X_1 and X_2 and apply Lemma <ref>, we see that Vor_X(p) = Vor_X_1(p) ∩ Vor_X_2(p). But, since p is a smooth point of X_1 and X_2, the Voronoi cells Vor_X_1(p) and Vor_X_2(p) are each contained in their respective normal directions, which are distinct. Therefore, Vor_X_1(p) ∩ Vor_X_2(p) = p.
    

  * If p is a tacnode, two branches of the curve meet at p where they share a tangent direction, and thus a normal line. We can choose ϵ>0 so that we can separate X ∩ B(p, ϵ) into subsets X_1 and X_2 corresponding to the two branches. Both Vor_X_1(p) and Vor_X_2(p) are 1-dimensional subsets of the same normal line. By Lemma <ref>, Vor_X(p) is a 1-dimensional subset of this normal line. 


  * Suppose p is an isolated point. Then there is a ball B(p,r) centered at p such that the ball contains no other points of the curve X. Therefore, the ball B(p,r/2) is entirely contained in Vor_X(p), so it is 2-dimensional.





In this example we illustrate why Theorem <ref> fails when the curve has a singular point. From this example it will be clear that the singular points must be included in the samples A_ϵ, and it turns out that this condition is enough to extend Theorem <ref> to the singular case.

Consider the curve defined by the equation y^2 = x^3. In <cit.> the authors give equations for the Voronoi cell of the cusp at the origin. This region is

    Vor_y^2=x^3((0,0)) = {(x,y) ∈ℝ^2  :  
    27 y^4 + 128 x^3 + 72 x y^2 + 32 x^2 + y^2 + 2 x ≤ 0
    }.

In Figure <ref> we give three  ϵ-approximations of the curve and the corresponding Voronoi decompositions. Let ϵ = 1/n. The points in the ϵ-approximation A_ϵ are given by:

    A_ϵ = { (j/n, ± (j/n)^3/2 )}_j = 1^∞.

As we can see in Figure <ref>, there is no sequence of cells converging to Vor_y^2=x^3((0,0)) because the x-axis, present due to the symmetrical nature of the sample, always divides the Voronoi cell.





We now are able to expand Theorem <ref> to include singular varieties.



Let X ⊂ℝ^2 be a compact algebraic curve and {A_ϵ}_ϵ↘ 0 a sequence of ϵ-approximations with the singular locus Sing(X) ⊂ A_ϵ for all ϵ. Then every Voronoi cell of X is the Wijsman limit of a sequence of Voronoi cells of A_ϵ.



By <cit.>, the Voronoi function is always lower-semicontinuous. So, we must show that condition (<ref>) in Definition <ref> holds. That is, we need that for all p ∈ X, there is a sequence p_ϵ∈ A_ϵ with p_ϵ→ p such that for any x ∈ Vor_X(p) there is an x_ϵ∈ Vor_A_ϵ(p_ϵ) with x_ϵ→ x. We distinguish the cases when p is smooth and singular.

If p is a smooth point on X, and x ∈ Vor_X(p), there exists an ϵ such that x and p are both in the Voronoi cell Vor_A_ϵ(p_ϵ) for some p_ϵ.

Suppose now that p ∈ X is a singular point. We wish to show that there is a sequence of Voronoi cells converging to Vor_X(p), and we take the sequence Vor_A_ϵ(p). To establish convergence, it is now enough to show that for all x ∈ Vor_X(p), there is an x_ϵ∈ Vor_A_ϵ(p) with x_ϵ→ x. Since x ∈ Vor_X(p), we have that x is closer to p than it is to any other point in X. So, in particular, x ∈ Vor_A_ϵ(p). 

Now we have shown that for each p ∈ X, condition (<ref>) in Definition <ref> holds. Therefore, for each p ∈ X, we have sequences of Voronoi cells which are convergent to Vor_p(X) in the hit-miss topology. Since convergence in the hit-miss topology and Wijsman convergence are equivalent, every Voronoi cell of X is the Wijsman limit of a sequence of Voronoi cells of the A_ϵ.



This concludes the proof of Theorem <ref>.

 




§ MEDIAL AXIS
 

Let X = V(F)⊂ℝ^2 be a smooth algebraic plane curve.

We now study the medial axis of X, as defined in Definition <ref>. 

The Zariski closure of the medial axis is an algebraic variety which has the same dimension as the medial axis. We can obtain equations in variables x,y for the ideal I of a variety containing the Zariski closure of the medial axis in the following way.

Let (s,t) and (z,w) be two points on X. Then, s,t,z, and w satisfy the equations 

    F(s,t)=0 and   F(z,w)=0 .
 
If (x,y) is equidistant from (s,t) and (z,w) then

    (x-s)^2+(y-t)^2=(x-z)^2+(y-w)^2 .
 
Furthermore, (x,y) must be a critical point of the distance function from both (s,t) and (z,w). Thus we require that the determinants of the following 2 × 2 augmented Jacobian matrices vanish: 

    [ x-s y-t; F_s F_t ] ,
         [ x-z y-w; F_z F_w ] ,

where F_s,F_t,F_z and F_w denote the partial derivatives of F(s,t) and F(z,w), respectively. 
Let

    I =  ⟨    
    F(s,t), F(z,w), (x-s)^2+(y-t)^2-(x-z)^2-(y-w)^2,
        (x-s)F_t - (y-t)F_s , (x-z)F_w - (y-w)F_z ⟩.


Then,
J =  (  I : (s-z,t-w)^∞ ) ∩ℝ[x,y]
is an ideal whose variety contains the Zariski closure of the medial axis.






We now study the medial axis from the perspective of Voronoi cells.
It has been observed that an approximation of the medial axis arises as a subset of the Voronoi diagram of finitely many points sampled densely from a curve <cit.>. We now discuss theoretical results given in <cit.> about the convergence of medial axes.

Let X be a compact smooth algebraic plane curve, and let A_ϵ be an ϵ-approximation of X. A Voronoi cell Vor_A_ϵ(a_ϵ) for a_ϵ∈ A_ϵ is polyhedral, meaning it is an intersection of half-spaces.



 For sufficiently small ϵ, exactly two edges of Vor_A_ϵ(a_ϵ) will intersect X <cit.>. We call these edges the long edges of the Voronoi cell, and all other edges are called short edges. An example is given in Figure <ref>.




In this case, let Ŝ_ϵ(a_ϵ) denote the union of the short edges and vertices of the Voronoi cell Vor_A_ϵ(a_ϵ). An ϵ-medial axis approximation is the set of all short edges

    Ŝ_ϵ = ⋃_p ∈ A_ϵŜ_ϵ(p).


(<cit.>)
Let X be a compact smooth algebraic plane curve. The medial axis approximations Ŝ_ϵ converge to the Euclidean closure of the medial axis.


The following corollary shows the relationship between the medial axis, parts of Voronoi diagrams, and maximally inscribed circles for ϵ-approximations.


Let {A_ϵ}_ϵ↘ 0 be a sequence of ϵ-approximations of a compact smooth algebraic curve X ∈ℝ^2.


    
  * The collection of vertices of the Voronoi diagrams of the A_ϵ converge to the medial axis. 
    
  * The collection of centers of maximally inscribed discs of the A_ϵ converge to the medial axis.



This is a consequence of Theorem <ref>, Theorem <ref>, and Proposition <ref>.




In Figure <ref> we display the centers of maximally inscribed circles, or equivalently circumcenters of the Delaunay triangles, for an ϵ-approximation of the butterfly curve where 898 points were sampled. In Figure <ref> we show the short edges of Voronoi cells from an ϵ-approximation of the butterfly curve where 101 points were sampled.







The medial axis plays an important role in applications for understanding the connected components and regions of a shape. As such, it is a very well-studied problem in computational geometry to find approximations of the medial axis from point clouds. 
A survey on medial axis computation is given in <cit.>.















§ CURVATURE AND THE EVOLUTE


Curvature of plane curves, osculating circles and evolutes have interested mathematicians since antiquity. 

We refer readers to works of Salmon in the 19th century <cit.> and to modern lectures by Fuchs and Tabachnikov outlining this history <cit.>. 




We now discuss the minimal radius of curvature of a plane curve. This is one of the two quantities which determines the reach, see Equation <ref>. 



Let X⊂ℝ^2 be an algebraic curve and p ∈ X be a smooth point that is not a point of inflection. The osculating circle at p is the circle that passes through p and an additional pair of points infinitesimally close to p. The center of curvature at p is the center of the osculating circle at p. The radius of curvature at p is the distance from p to its center of curvature. 



An alternative definition of center and radius of curvature can be given using envelopes. 



The envelope of a one-parameter family of algebraic plane curves given implicitly by 
F(x,y,t)=0
is a curve that touches every member of the family tangentially. The envelope is the variety defined by the ideal 

    ⟨∂ F/∂ t,   F(x,y,t) ⟩∩ℝ[x,y].
 


The envelope of the family of normal lines parametrized by the points of the curve is called its evolute. Equivalently, the evolute is the locus of the centers of curvature.
A generalization of the evolute to all dimensions is called the ED discriminant, and is studied in <cit.>.
The article shows that for general smooth algebraic plane curves, the degree of the evolute is
3d(d-1) <cit.>. See also <cit.> for a discussion of the degree of the evolute.

We now derive a formula for the square of the radius of curvature of a plane curve at a point. Our derivation follows Salmon <cit.>. 



<cit.>
Let X=V(F(x,y)) ∈ℝ^2 be a smooth curve of degree d. The square of the radius of curvature at a point (x_0,y_0) ∈ V(F) that is not a point of inflection is given by following expression in the partial derivatives of F evaluated at (x_0,y_0):

    R^2=  (F_x^2+F_y^2)^3/(F_xxF_y^2-2F_xyF_xF_y+F_yy F_x^2)^2.




The equation of a normal line to X at a point (x,y) ∈ X  in the variables (α, β) is  

    F_y(α-x)-F_x(β-y)=0.


The total derivative of the equation for the normal line is 

    (F_xy+F_yydy/dx)(α-x)-(F_xx+F_xydy/dx)(β-y)-F_y+F_xdy/dx=0.


The total derivative of F(x,y) is 

    F_x(x,y)+F_ydy/dx = 0.


 The equations (<ref>), (<ref>) are a system of two linear equations in the unknowns {α, β}. We solve this system to obtain expressions for α and β in terms of x, y, and dy/dx. We substitute in for dy/dx the expression given by (<ref>). The center of curvature of X at a point (x,y) ∈ X is given by the coordinates (α, β), which are now expressions in x and y. 
 
 The squared radius of curvature at a point (x,y) is its squared distance to its center of curvature (α, β), so we have 
R^2=(α-x)^2+(β-y)^2.
Substituting in the equations for α and β, we find 
 
    R^2=  (F_x^2+F_y^2)^3/(F_xxF_y^2-2F_xyF_xF_y+F_yy F_x^2)^2.

 We note that the denominator evaluates to zero only at points of inflection.
 




The degree of critical curvature of a smooth algebraic curve X ⊂^2 is the degree of the variety obtained by intersecting the Zariski closure X⊂ℙ_ℂ^2  with the variety of the total derivative of the equation for the squared radius of curvature. If X ⊂ℝ^2 is a smooth, irreducible algebraic curve of degree greater than or equal to 3, then the intersection consists of finitely many points, called the points of critical curvature. Thus the degree of critical curvature of X gives an upper bound for the number of real points of critical curvature of X. 




In the differential geometry literature, points of critical curvature are called vertices.




Let X ⊂ℝ^2 be a smooth, irreducible algebraic curve. Then the degree of critical curvature of X is 6d^2-10d. 


We remark that the critical points of curvature of X give cusps on the evolute <cit.>. 
That is, if a normal line is drawn through a point of critical curvature on a curve, then the normal line will pass through a cusp of the evolute. 
In addition, the evolute of a curve of degree d has d cusps at infinity <cit.>. Thus the evolute of a general plane curve of degree d has 6d^2-10d+d=6d^2-9d cusps <cit.>. In Figure <ref>, we picture the evolute, the butterfly curve, and the pairs of critical curvature points on the butterfly curve with their corresponding cusp on the evolute.





Consider the butterfly curve. Using the above description, we can compute the 56 points of critical curvature using  <cit.>. Twelve of these points are real, and they are plotted in Figure <ref>. The maximal curvature is approximately 9.65. This is achieved at the lower left wing of the butterfly.






We now describe how to recover the curvature at a point from the Voronoi cells of a subset of a curve X. In applications, Voronoi-based methods are used for obtaining estimates of curvature at a point. 
An overview of techniques for estimating curvature of a variety from a point cloud is given in <cit.>. Further, there are also  Delaunay-based methods for estimating curvature of a surface in three dimensions <cit.>. 



Let X be a smooth, irreducible plane curve of degree at least 3 and p∈ X a point that is not a critical point of curvature. 
Let δ be less than the distance to the critical point of curvature nearest to p, and
let B(p, δ) be a ball of radius δ centered at p. Then

    
  * 
    The Voronoi cell Vor_X ∩ B(p, δ)(p) is a ray. The distance from p to the endpoint of this ray is the radius of curvature of X at p.
    
  * 
    Consider a sequence of ϵ-approximations A_ϵ of X ∩ B(p,δ). Let a_ϵ be a point such that p ∈ Vor_A_ϵ(a_ϵ), and let d_ϵ be the minimum distance from a_ϵ to a vertex of Vor_A_ϵ(a_ϵ). Then, the sequence d_ϵ converges to the radius of curvature of X at p.



The Voronoi cell Vor_X ∩ B(p, δ)(p) is a subset of the line normal to X ∩ B(p, δ) at p. This line has an endpoint either at the center of curvature of p or at a point where it intersects the normal space of a distinct point p' in X∩ B(p,δ). The point where the normals at p and p' intersect is contained in the Voronoi cell with respect to X ∩ B(p, δ) of each of them, so in particular X ∩ B(p,δ) has a nonempty medial axis. This medial axis has an endpoint which corresponds to a point of critical curvature <cit.>. This contradicts the constraint on δ. Therefore, the endpoint of the Voronoi cell Vor_X ∩ B(p, δ)(p) is the center of curvature of p. This concludes the proof of (<ref>).

For (<ref>), we know that the sequence Vor_A_ϵ(a_ϵ) is Wijsman convergent to Vor_X ∩ B(p, δ)(p) by Theorem <ref>. Denote by V_ϵ the set of vertices of Vor_A_ϵ(a_ϵ). By Corollary <ref>, we also have that the sets V_ϵ are Wijsman convergent to the endpoint of Vor_X ∩ B(p,δ)(p), which we call p'. By the definition of Wijsman convergence, this means that for any x ∈ℝ^2, d_w(x,V_ϵ) → d_w(x,p'). By the definition of d_w, we have d_w(p,V_ϵ) = d_ϵ and d_w(p,p') is the radius of curvature of p. This concludes the second part of the proof.



The evolute E of a plane curve is the locus of all centers of curvature of the curve. Therefore, to find the evolute using Voronoi cells we may splice the curve into sections and apply Theorem <ref>.
Let X be compact and irreducible of degree greater than or equal to 3. Let C ⊂ X denote the points of locally maximal curvature and E_C denote the centers of curvature corresponding to points in C. Then X \ C consists of finitely many components X \ C= X_1∪⋯∪ X_n. Let τ denote the reach of X, and cover each X_i by balls B_i,j of radius less than τ.
Let E_i,j denote the collection of vertices of Voronoi cells of X_i ∩ B_i,j. Then by Theorem <ref>, E \ E_C= ∪_i,jE_i,j. 
Furthermore, for ϵ-approximations A_ϵ,i,j of X_i ∩ B_i,j, the union over i,j of their Voronoi vertices will converge to E \ E_C by Theorem <ref>. To find the evolute E one need only to add the finite set of points E_C. 








§ BOTTLENECKS







As in the colloquial sense of the word, a bottleneck refers to a narrowing of a variety, or a place where it gets closer to self-intersection. 
Consider a smooth algebraic variety X ⊂^n. We define a⊥ b by ∑_i=1^na_ib_i=0 for
a=(a_1,…,a_n), b= (b_1,…,b_n) ∈^n.
For a point x ∈ X,
let (T_xX)_0 denote the embedded tangent space of X translated to
the origin. Then the Euclidean normal space of X at x is
defined as N_xX={z ∈^n:(z-x) ⊥ (T_xX)_0}.



A bottleneck of a smooth algebraic variety X ⊂^n is a pair of
distinct points (x,y) ∈ X × X such that xy⊆ N_xX ∩ N_yX, where xy is the line spanned by x and y. 

 

We note that bottlenecks are given not only by the narrowest parts of the variety, but also by maximally wide parts of the variety, as our algebraic definition considers all critical points rather than just the minimums. 
The bottlenecks of the butterfly curve are shown in Figure <ref>.



The narrowest bottleneck distance ρ of a variety X ⊂^n is 

    ρ(X)=inf_(x,y)  a bottleneck d(x,y)

where d(x,y) is the Euclidean distance between x and y. 





We will now describe the bottleneck locus in
^2n which consists of the bottlenecks of X <cit.>.  Let (f_1,…,f_k)
⊆[x_1,…,x_n] be the ideal of X. Consider the ring
isomorphism ϕ: [x_1,…,x_n] →[y_1,…,y_n] defined by x_i ↦
y_i and let f'_i=ϕ(f_i). Then f_i and f_i' have gradients
∇ f_i and ∇ f'_i with respect to {x_1,…,x_n} and
{y_1,…,y_n}, respectively. The augmented
  Jacobian J is the following matrix of size (k+1) × n with
entries in R=[x_1,…,x_n, y_1,…, y_n]:

    J = 
    [   y-x; ∇ f_1;     ⋮; ∇ f_k ] ,

where y-x is the row vector (x_1-y_1,…,x_n-y_n). Let N
denote the ideal in R generated by (f_1,…,f_k) and the
(n-dim(X)+1) × (n-dim(X)+1) minors of J. Then the points (x,y) of the
variety defined by N are the points (x,y) ∈ X × X ⊂^2n such
that y ∈ N_xX. In the same way we define a matrix J' and an
ideal N' ⊆ R by replacing f_i with f'_i and ∇ f_i
with ∇ f'_i.

The bottleneck locus B is the variety

    B = V((N + N'): ⟨ x-y⟩ ^∞) ⊂ X × X ⊂^2n.

The saturation removes the diagonal, as (x,y) is not a bottleneck if x=y. 






Next, we give the bottleneck degree, which is 
a measure of the complexity of computing all bottlenecks of an algebraic variety. We refer readers to <cit.> for a discussion of the numerical algebraic geometry of bottlenecks. Under suitable genericity assumptions described in <cit.>, it coincides with twice the number of bottlenecks of the complexification of X. The factor of 2 is attributed to the fact that in the product X × X, the points (x,y) and (y,x) are distinct, though they correspond to the same pair of points in X. 

 

Under certain genericity assumptions, the degree of the bottleneck locus of a smooth algebraic curve X ⊂^2 of degree d is

    d^4-5d^2+4d.
 


A degree formula for the bottleneck locus of varieties of any dimension is provided in <cit.>. The proof applies the double point formula from intersection theory to a map taking the variety to a variety of its normals. 



We now compute the bottlenecks for the quartic butterfly curve b(x,y) = 0. Theorem <ref> predicts that there are 192/2 = 96 bottlenecks. Using the description above and  <cit.>, we obtain the 96 bottleneck pairs. Of these, 22 are real. We show them in Figure <ref>.






We now study bottlenecks from the perspective of Voronoi cells. For a smooth point p in a algebraic curve X ⊂ℝ^2, the Voronoi cell Vor_X(p) is a 1-dimensional subset of the normal line to X at the point p. Therefore, the normal direction can be recovered from the Voronoi cell Vor_X(p).
For sufficiently small ϵ, an ϵ-approximation A_ϵ of X will have Voronoi cells whose long edges approximate the normal direction. More precisely, by Theorem <ref>, if a_ϵ∈ A_ϵ is the point such that p ∈ Vor_A_ϵ(a_ϵ), then the directions of the long edges of Vor_A_ϵ(a_ϵ) converge to the normal direction at p.
We remark here that the problem of estimating normal directions from Voronoi cells is well-studied, and numerous efficient, robust algorithms exist <cit.>.



As in Definition <ref>, two points x,y ∈ X form a bottleneck if their normal lines coincide. This implies that the line connecting them contains both Vor_X(x) and Vor_X(y). 



Let A_ϵ be an ϵ-approximation of an algebraic curve X ⊂ℝ^2. We say a pair x_ϵ, y_ϵ∈ A_ϵ is an approximate bottleneck reach candidate if the line x_ϵ y_ϵ joining x_ϵ and y_ϵ meets each of Vor_A_ϵ(x_ϵ) and Vor_A_ϵ(y_ϵ) at short edges of those cells.




In Figure <ref> we show the approximate bottleneck reach candidates for 348 points sampled from the butterfly curve. The following result gives conditions for bottleneck pairs to be the limit of approximate bottleneck reach candidates.




Let {A_ϵ}_ϵ↘ 0 be a sequence of ϵ-approximations of a smooth algebraic curve X ⊂ℝ^2. If x,y is a bottleneck pair of X such that the midpoint of the line segment xy is on the medial axis of X and this is the only place where  xy meets the medial axis, then there are sequences x_ϵ, y_ϵ∈ A_ϵ of approximate bottleneck reach candidates converging to x and y. In particular, if x,y is a bottleneck pair of X that achieves the reach, then there are sequences x_ϵ, y_ϵ∈ A_ϵ of approximate bottleneck reach candidates converging to x and y.





Let x,y be a bottleneck pair and suppose that the midpoint of the line segment xy joining x and y is on the medial axis of X and this is the only place where  xy meets the medial axis. We note that this condition holds in particular if x and y are a bottleneck pair that achieves the reach. 
Then xy intersects the short edges of exactly two Voronoi cells Vor_A_ϵ(x_ϵ) and Vor_A_ϵ(y_ϵ) in a point v_ϵ. Then, x_ϵ and y_ϵ form an approximate bottleneck reach candidate by definition. We must then show that the sequence x_ϵ converges to x and the sequence y_ϵ converges to y.

Since v_ϵ is in the normal space of x,
there exists a neighborhood of x such that
the nearest point of the intersection of this neighborhood and X to v_ϵ is x. So for ϵ smaller than the radius of this neighborhood, one of the two points in A_ϵ on either side of x as one moves along X must be the one whose Voronoi cell contains v_ϵ. Since x_ϵ is the point whose Voronoi cell contains v_ϵ, we have that x_ϵ is one of the two closest points in A_ϵ to x, meaning that d(x,x_ϵ) ≤ϵ. Hence, x_ϵ converges to x. Similarly, y_ϵ converges to y.









§ REACH

















We may find the reach of the butterfly curve by taking the minimum of half the narrowest bottleneck distance and the minimum radius of curvature. This is shown in Figure <ref>. From the computations in Example <ref>, we find that the narrowest bottleneck distance is approximately 0.251. Meanwhile, from Example <ref>, we find that the minimum radius of curvature is approximately 0.104. Therefore, the reach of the butterfly is approximately 0.104.


In previous sections, we describe how the reach is the minimum of the minimal radius of curvature and half of the narrowest bottleneck distance. We also give equations for the ideal of the bottlenecks and for the ideal of the critical points of curvature. We now give  <cit.> code to compute these ideals for smooth algebraic curves X ⊂^2.
Here, the expression for  comes from using Lagrange multipliers to find critical points of the affine radius of curvature subject to the constraint given by the curve f. 
Finding the points in these ideals, using for example  <cit.>, and taking appropriate minimums gives the reach of X. See <cit.> for an alternate technique for computing the reach in Julia. 







Alternatively, one can estimate the reach from a point sample. The paper <cit.> provides a method to do so. We provide a substantially different method that relies upon computing Voronoi and Delaunay cells of points sampled from the curve. We have already discussed how to approximate bottlenecks and curvature using Voronoi cells. This gives the following Voronoi-based Algorithm <ref> for approximating the reach of a curve.



The reach is equivalently defined as the minimum distance to the medial axis, which suggests the following Delaunay-based Algorithm <ref> for estimating the reach. This algorithm is susceptible to sample error, and to give accurate results would require more sophisticated techniques.



The approximate methods can be used with curves of higher degree, while the symbolic methods are hard to compute for curves with degrees even as low as 4, but give a more accurate estimate for the reach. This suggests that more work can be done to develop fast and accurate methods to compute the reach of a variety. 






§ ACKNOWLEDGEMENTS


We thank Paul Breiding, Diego Cifuentes, Yuhan Jiang, Daniel Plaumann, Kristian Ranestad, Rainer Sinn, Bernd Sturmfels, and Sascha Timme for helpful discussions. We thank the referees for their thoughtful comments and suggestions for improving the paper. Research on this project was carried out while the authors were based at the Max Planck Institute for Mathematics in the Sciences (MPI-MiS) in Leipzig, Germany. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE 1752814. Any opinions, findings, and
conclusions or recommendations expressed in this material are those of the authors and do not
necessarily reflect the views of the National Science Foundation.




plain


