A boundary-partition-based Voronoi diagram of d-dimensional balls: definition, properties, and applications

In computational geometry, different ways of space partitioning have been developed, including the Voronoi diagram of points and the power diagram of balls. In this article, a generalized Voronoi partition of overlapping d-dimensional balls, called the boundary-partition-based diagram, is proposed. The definition, properties, and applications of this diagram are presented. Compared to the power diagram, this boundary-partition-based diagram is straightforward in the computation of the volume of overlapping balls, which avoids the possibly complicated construction of power cells. Furthermore, it can be applied to characterize singularities on molecular surfaces and to compute the medial axis that can potentially be used to classify molecular structures.


Introduction
The Voronoi diagram [1][2][3] is a partition of a Euclidean plane into regions, based on the distance to points in a specific subset of the plane. In any dimension d, given a finite set of points {p 1 , p 2 , . . . , p N } in E d (the d-dimensional Euclidean space), the corresponding Voronoi cell R V (p i ) consists of every point whose distance to p i is not greater than its distance to any other point p j . That is to say, the Voronoi region is given by where the notation | · | denotes the Euclidean norm. Each Voronoi cell is generated by the intersection of half-spaces, and hence, is a convex polygon (see Fig. 1 (left) for a graphical illustration). The power diagram [4,5], also called the Laguerre-Voronoi diagram, provides another partition of the plane into polygonal cells with respect to a finite set of circles. For a finite set of spheres {S 1 , S 2 , . . . , S N } in E d with d ≥ 2 (circles in E 2 ), the power region R P (S i ) consists of all points whose power distances to S i are not larger than their power distances to any other sphere S j . The power distance from a point x ∈ E d to a sphere S i with center c i and radius r i is defined as The power region R P (S i ) is then given by It is worth to mention that the spheres could be overlapping and the power distance can be negative. The power diagram can be seen as a generalized Voronoi diagram, in the sense that the Euclidean norm in Eq. (1) is replaced with the power distance dist P (·, ·). If r i = 0 for each i, the power diagram degenerates to be the classical Voronoi diagram. Figure 1 (middle) shows an example of the power diagram for three discs. The volume of the union of d-balls can be computed by summing up the volume of each power cell (see [6,7] for details). However, computing the power cells could be tricky which involves characterizing the planar bases of pyramids [7]. In the general case, given a Euclidean subspace X ⊂ E d endowed with a distance function and a finite number of nonempty subsets {A 1 , A 2 , . . . , A N } in X, the corresponding Voronoi region R V (A i ) is the set of each point in X whose distance to A i is not greater than their distance to any other set A j . That is to say, the Voronoi region R V (A i ) is given by where the distance function between a point x and a set A defined as Most commonly, each subset A i is taken as a point and its corresponding Voronoi region R V (A i ) is consequently a polyhedron in E d . In particular, the Voronoi diagram in E 3 can be computed and visualized by some softwares such as the Voro++ [8] and the CGAL [9]. In addition to the power diagram, there are some other generalized Voronoi diagrams, such as the Möbius diagram [10], the anisotropic diagrams [11], the Apollonius diagram [12,13], and the centroidal Voronoi tessellations [14].
In this article, we propose a special Voronoi diagram of d-balls, called the boundary-partition-based (BPB) diagram, which is based on the boundary partition of the union of d-balls. The initial idea of this diagram started from the characterization of molecular surfaces in E 3 (see [15]). Figure 1 (right) provides a simple example of this diagram of three discs in the plane, whose boundary is divided into three open circular arcs {γ 1 , γ 2 , γ 3 } and three intersection points {x 1 , x 2 , x 3 }. We take {γ 1 , γ 2 , γ 3 , x 1 , x 2 , x 3 } as the sets {A i } in Eq. (4) for the generic Voronoi diagram, to obtain six BPB cells. This gives a partition of the union of discs, consisting of circular sectors and polygons. To compute the area of these discs, one can sum up simply the areas of three circular sectors (respectively in red, yellow, and blue) and one big polygon (in gray). Note that the area of a polygon is easy to compute using the Gauss-Green theorem, as its boundary is composed of line segments. Comparing to the power diagram, the BPB diagram is convenient for volume (area) computation, which will be presented later. In addition, it can be used to characterize the singularities of molecular surfaces.
The outline of this paper is as follows. In Section 2, we divide the boundary into patches of different dimensions and study the properties of the BPB cells in E d , by analyzing the signed distance from any point to the boundary of a union of d-balls. Then, in Section 3, we present the volume formula of a union of d-balls based on the boundary components. In Section 4, we introduce the application of the BPB diagram to characterize molecular surfaces. Finally, we draw some conclusions in Section 5.

Boundary-partition-based Voronoi diagram
We consider a finite set of (d −1)-dimensional spheres {S 1 , S 2 , . . . , S N } in E d , where the (d − 1)-sphere S i has center c i ∈ E d and radius r i for each 1 ≤ i ≤ N. The corresponding d-balls are consequently denoted by {B 1 , B 2 , . . . , B N } where B i = B(c i , r i ). In this article, the notation B(c, r) denotes the open ball with center c and radius r. The union of these d-balls is denoted by Ω, that is, (6) and the boundary of this union is denoted by Γ := ∂Ω. The open exterior region outside Ω is denoted by As the boundary Γ is a closed and compact set, for any point p ∈ E d , there exists at least one closest point on Γ to p. The signed distance function f Γ with respect to Γ is then given as follows As a consequence, the sets Ω, Γ , and Ω c can be mathematically written as Ω = Remark 1 Another way to characterize Ω, Γ , and Ω c is to use the signed distance By defining the following function

Boundary partition
As mentioned previously, the boundary Γ is composed of spherical patches of different dimensions from 0 to d − 1. We will provide a rigorous characterization in this subsection.
To do this, we define an index mapping I as follows: where I(x) collects all indices of the spheres containing x and m is the number of these spheres. Then, we define the intersection set of these spheres as follows: where the subscript i = I(x) = {i 1 , i 2 , . . . , i m } depends on x (see an example in Fig. 2). Further, the affine space generated by the associated centers is denoted by With the above notations, we propose the following lemma:

Lemma 1 In the case when S i is non-empty and dim(Λ
or a point (in the tangent case), where dim(·) denotes the dimension. Note that the 0-dimensional sphere is a pair of points. Furthermore, in the case when S i is non-empty and dim(Λ i ) = d, S i is a point.
Proof Let us take an arbitrary point x 0 ∈ S i . We rewrite S i in the following form: where Here, P i is the intersection of m − 1 hyperplanes that contain x 0 . It is an affine space In the case of dim(Λ i ) ≤ d − 1, we then have dim (P i ) ≥ 1. Note that S i is the intersection of the sphere S i 1 and the affine space P i . If S i is non-empty, then it is either a sphere with dim(S i ) = d − dim(Λ i ) − 1, or a single point in the degenerate case when P i is tangent to S i . We mention that in this degenerate case, S i is a point contained in Λ i with dim(Λ i ) ≤ d − 1. Furthermore, in the case of dim(Λ i ) = d, we have dim(P i ) = 0 and P i is consequently the point x 0 . If S i is nonempty, then S i is just the point x 0 .
For a given set of indices i ⊆ {1, 2, . . . , N}, if dim(S i ) ≥ 1, we can define the following set: which is open in S i . In this case, we can further divide Γ i as follows: where γ i,j ) denotes the dimension. This means that Γ i has been divided into different k-patches. For the sake of simplicity, we can reorder all patches {γ (k) i,j } based on the dimension k, by replacing the subscripts (i, j) with only one subscript i. As a consequence, the whole boundary Γ is classified into a set of patches {γ where n k is the number of k-patches and γ i , we know that I(x) remains the same for any point x ∈ γ (k) i . Therefore, we can generalize the definition of I as follows: is an arbitrary point on γ (k) i . Next, we analyze the signed distance f Γ to the boundary Γ , which involves in finding one closest point to any given point p (note that the uniqueness of the closest point is not guaranteed).

Analysis of the signed distance function
We want to analyze the signed distance f Γ from an arbitrary point to the boundary Γ of the union of d-balls. We first consider the case when the point lies outside the union and give the following lemma. The proof of this lemma is trivial and therefore skipped here.

Lemma 2 For any point
We now focus our attention to the case when p ∈ Ω. In fact, this case is studied from another perspective, in the sense that for any point x on Γ , we study the set of all points in Ω that have x as a closest point on Γ . We therefore define a mapping R such that ∀x ∈ Γ , which represents the region consisting of the points having x as a closest point. It therefore holds for any x ∈ Γ and p ∈ R(x), f Γ (p) = −dist (p, Γ ) = −|p − x|.
Further, the convexity of the set R(x) is ensured according to the following lemma.
. . , i m }. Here, the notation conv denotes the convex hull of a set of points.
Proof We first need to prove that ∀p 1 , To do this, we construct a function as follows: As a consequence, x is a closest point on Γ to p if and only if h(p) = 0. As an obvious fact, we have h(p) ≤ 0, ∀p ∈ E d . Since x is a closest point on Γ to p 1 and to p 2 , we have h( This means that x is also a closest point to p 0 . To prove that p 0 ∈ R(x), we should show further that p 0 ∈ Ω. As x is a closest point to p 1 , we know that the ball centered at p 1 with radius |p 1 − x| is covered by Ω, i.e., Similarly, we have B(p 2 , |p 2 − x|) ⊆ Ω. Notice that the line segment p 1 p 2 is covered by the union of these two balls (intersecting at x), which implies that p 0 ∈ p 1 p 2 ⊆ Ω.
Since p 0 has x as a closest point and p 0 ∈ Ω, we therefore have p 0 ∈ R(x). So far, we have proved that R(x) is convex. Further, if there exists a set of spheres {S i t } t=1,...,m each containing x ∈ Γ , it is obvious that x is a closest point of c i t . In addition, we know that x ∈ Γ is a closest point to itself. Due to the convexity of R(x), we then have conv(x, Remark 2 It is well-known that the Voronoi cells of a finite number of points are convex. In Lemma 3, we have actually proved that this convexity also holds for an infinite number of points (forming the boundary surface).
represents a convex cone (as in linear algebra) with apex x.
(2) If there exists another point x ∈ Γ satisfying x = x and x ∈ S i , then R(x) = conv(x, c i ).
Proof (1) According to Theorem 19.1 in the book [16], we can write where is a real number, N 0 denotes the number of half-spaces, and (·, ·) denotes the Euclidean scalar product in E d equivalent to the · notation. Each half-space H s corresponds to a hyperplane P s defined by Without loss of generality, we suppose that |n s | = 1 for each half-space H s . Furthermore, each hyperplane P s is supposed to contain at least one ray starting from x in a direction among {v 1 , v 2 , . . . , v m }. In fact, the half-space H s can be removed in Eq. (23) if P s does not contain any such ray. To prove the first statement of the theorem, it suffices to show that ∀p ∈ R(x) and ∀1 ≤ s ≤ N 0 , it follows that p ∈ H s . Let p ∈ R(x) and s be fixed. As mentioned above, P s contains a ray starting from We now construct a small curve ζ s (x) starting from x in the direction of n s lying on Γ , of the form where α(x) ≥ 0 is a function with respect to x satisfying α(0) = 0, v is a nonzero vector in E d , and ε is a sufficiently small positive number. In order to choose v, we define the following nonempty set: As a consequence, the vector v and the function α(x) can be constructed as follows: and where The aboveconstructed v has the following two properties: and With the above construction of ζ s (x) and the properties of v, we can now state the following claim.
The proof of this claim is presented in Appendix 1. Let us focus on the proof of the first statement in Theorem 1. For any point p ∈ R(x), we know that x is a closest point on Γ to p and ζ s (x) ∈ Γ, ∀x ∈ [0, ε]. Therefore, we have which yields that x − p, ζ s (0) ≥ 0. It is not difficult to find that ζ s (0) = n s and we then obtain (p, n s ) ≤ (x, n s ) = b s . This implies that p ∈ H s for each subscript s. The proof of the first statement in the theorem is complete.
(2) We now prove the second statement in Theorem 1. Suppose that there exists another point which yields that Furthermore, since x is a closest point on Γ to p, we have Since x = x, we consequently have the inequality m t=1 λ t ≤ 1, which means that p ∈ conv(x, c i ). Then, we have R(x) ⊆ conv(x, c i ).

Corollary 1 For each boundary component
Proof In the case of k ≥ 1, γ (k) i contains infinitely many points, each of which is contained and only contained by the set of spheres {S i 1 , S i 2 , . . . , S i m }. Therefore, the second statement in Theorem 1 can be applied.

Partition of the union of d -balls
In this part, we introduce the general concept of the BPB diagram, which gives a partition of E d . According to Lemma 1, the partition of Ω c is directly based on the simple function F (·) given by Eq. (8). Alternatively, one can also decompose Ω c using the power distance or even simply treat Ω c as one entire cell. Thus, we are only interested in the partition of Ω.
The mapping R maps any point x ∈ Γ to a subregion of Ω that collects all points having x as a closest point. We now generalize the definition 20 of R as follows: which maps any subset γ of Γ to a subregion of Ω such that each point in the subregion has a closest point in γ . Further, we define the BPB cell corresponding to γ (k) i as follows: where the distance function dist(·, ·) is given by Eq. (5). As a consequence, each BPB cell R BP (γ (k) i ) satisfies the following relationship: According to Theorem 1 and Corollary 1, in the case where i = I(γ where i = I(γ To better understand the cell R(γ (0) i ), we define the following subregion of Ω by R 0 will be further analyzed later in Section 3.
Here, we provide two examples of the BPB diagram. Figure 3 provides a partition of some discs in E 2 , respectively, obtained from the power diagram (computed by F. McCollum's package [17]) and the proposed BPB diagram. In the BPB diagram, the union of these discs is divided into circular sectors and polygons (constituting R 0 ). Note that the boundary of R 0 is composed of line segments with the disc centers and the intersection points as endpoints. As a consequence, the area of R 0 can be obtained directly using the Gauss-Green theorem, while the power cells could be complicated to compute. Figure 4 provides an example of the power diagram and the BPB diagram of three intersected balls in E 3 . The union of these 3-balls is divided into spherical sectors (in red), double-cone cells (in yellow), and tetrahedrons (in blue), respectively, corresponding to the 2-patches, 1-patches (circular arcs) and 0-patches (intersection points) on Γ . The volumes of BPB cells are can be computed conveniently, similarly to the 2D case, which will be explained with details in Section 3.
In summary, the group of subregions {R(γ i . Therefore, this newly proposed diagram allows to find efficiently a closest point on Γ to p and consequently allows to compute the signed distance f Γ (p). Since Lemma 3, Theorem 1, and Corollary 1 provide an accurate description of any cell R γ (k) i , we emphasize that the volume (resp. area) of the union of balls (resp. discs) can be computed based on the BPB diagram, which is explained in the next section.

Volume of the union of d-balls
The BPB diagram can be used to calculate the volume of Ω, given all components of Γ . In E 2 and E 3 , the boundary components are easy to compute, while this could be difficult in higher dimension (so as the power diagram). In E d with d ≥ 4, the BPB diagram builds a relationship between the boundary and the volume through the mapping R.

General formula
Consider an arbitrary k-patch γ is part of a k-sphere with center c (k) i and radius r (k) i (see Lemma 1). According to Eq. 33, we can compute where i = I(γ

Lemma 4 Let γ (k) i
be an arbitrary k-patch with 1 ≤ k ≤ d − 1 contained on a k-sphere. The d-volume of R(γ (k) i ) can be characterized as follows: where Vol (k) (γ ) denotes the k-volume of a k-dimensional surface γ , r i denotes the radius of the k-sphere containing γ (k) i and B(·, ·) is the Beta function.
Proof We define two sets is the center of the k-sphere containing γ (k) i . From Lemma 1, we then conclude that ⊥ . Further, and are respectively of dimension k and d − k − 1. According to Eq. 37, we can write R(γ (k) i ) as follows: (39) Since ⊥ , the volume infinitesimal dy can be written as We can consequently compute where B(·, ·) is the Beta function [18].
We now consider an arbitrary intersection point γ is an endpoint of some 1-patches (circular arcs) denoted by {γ i ) denotes the affine space defined in Eq. 12. In the degenerate case, as presented in the proof of Lemma 1. In this case, γ i is actually generated by a (d − 1)-sphere and a tangent affine space. According to Theorem 1, we have R(γ As a consequence, the volume of the BPB cell corresponding to a degenerate intersection point is This implies that the degenerate intersection points can be ignored in the computation of Vol (d) (R 0 ). In the nondegenerate case, γ (0) i is generated by the intersection of a (d − 1)-sphere and some line passing through the sphere. In this case, there exists some 1-patches having γ (0) i as an endpoint, implying that K i > 0. Further, it holds that according to Lemma 1. Then, we can denote the (d − 1)-dimensional face corresponding to γ (1) ij by where c I(γ (1) ij ) is a set of spherical centers given by Eq. 11 and F ij is actually a tetrahedron in some (d − 1)-hyperplane. We define the set of all nondegenerate intersection points as P 0 . Then, we propose the following lemma for computing Vol (d) (R 0 ) where R 0 is defined in Eq. 35.

Lemma 5
The volume of R 0 can be computed as where n ij denotes the outward-pointing normal vector of F ij . We make a convention that n ij = 0 when γ (0) i is a degenerate intersection point.
Proof Denote by R 0 the union of all BPB cells associated with nondegenerate intersection points in P 0 , that is, According to the definition Eq. 35 of R 0 and Eq. 43, we know that Claim The boundary of R 0 can be characterized as where F 0 is some subset on The proof of Claim 3.1 is presented in Appendix 2. According to the Gauss-Green theorem, we can further compute where dσ y denotes the surface measure. In the last equality, we use the fact that F ij lies in a hyperplane and that n ij · y is constant.
In summary, given the components of its boundary Γ , we obtain an explicit expression of the volume of the union of balls Ω according to Lemmas 4 and 5 as follows:

Analytical volume in 3D
We have given an explicit formula of the volume of Ω, which is based on the BPB diagram. In the cases of E 2 and E 3 , the different components of Γ are not difficult to compute. In this subsection, we consider the case of E 3 as an illustration. The boundary of the union of 3-balls is constituted by the intersection points, the circular arcs (or circles), and the spherical 2-patches. The length of a circular arc is easy to compute and the area of a spherical 2-patch can be computed by the Gauss-Bonnet theorem [19]. For the sake of completeness, we present here the explicit formula of the area of a spherical 2-patch γ with the notations in Fig. 5 as follows: (see [15] for details) where α j is the angle at vertex v i between two neighboring circular arcs e j −1 and e j , k e j is the geodesic curvature of e j , |e j | is the length of e j , and Area (γ ) is the area of γ . In addition, χ is the Euler characteristic of γ , which equals to 2 minus the number of loops forming the boundary of γ . From the volume formulation Eq. 51 in where r (2) i denotes the radius of the 2-patch γ (2) i , r (1) i denotes the radius of the circular arc γ (1) i , d (1) i denotes the distance between the centers of the two spheres generating γ (1) i , and Area(γ (2) i ) is computed according to Eq. 52. The faces {F ij } are 2D polygons enclosed by specific line segments (connecting the sphere centers and the intersection points) and in the common case are triangles. The above formula Eq. 53 is convenient to be computed.
The computation of the volume of molecular cavities is an elementary problem in biology and chemistry. There are plenty of works on it using the power diagram, such as [20][21][22][23][24]. Based on the proposed BPB diagram, we test 18 molecules in Matlab, with the geometry data derived from the Protein Data Bank (PDB), including caffeine, 1yjo, 1etn, 1b17, 101m, 2k4c, 3wpe, 1kju, 1a0t, 1a0c, 4xbg, 4cql, 5any, 4wht, 4qy1, 4by9, 4u8u, and 4y5z. The numerical results are presented in Table 1 and Fig. 6. We should mention that these volumes are exact if the machine errors are not taken into account. In addition, the run time appears to scale roughly linearly with respect to N.

Remark 4
The volume formula Eq. 53 of the BPB diagram only involves the computation of boundary components. In other words, any buried ball does not contribute to the volume. While the power diagram requires to compute more, because one has to compute not only the boundary components but also the boundaries of all power cells, which could be tricky and nonrobust due to the possible degenerate cases as pointed out in [7]. From this point of view, the BPB diagram is more convenient to be implemented and can save computational cost for the volume computation. N represents the number of balls (atoms). n 0 , n 1 , and n 2 represent the number of intersection points, the number of circular arcs (or circles), and the number of spherical patches on the boundary Fig. 6 The run time for volume computation with respect to the number of atoms, implemented in Matlab on a mac with processor 2.5-GHz Intel Core i7 Remark 5 In the high-dimensional space E d with d > 3, it is difficult to compute the surface volumes of the boundary components on Γ in a general case. As a consequence, one cannot compute the analytical volume easily based on the BPB diagram, unless the boundary components are known. Nevertheless, the same difficulty exists for the power diagram.
In the next section, we will discuss about another application of the BPB diagram for characterizing molecular surfaces.

Characterization of molecules
For the sake of completeness, we introduce briefly how the BPB diagram can be used to characterize molecules in E 3 .

Singularity problem
In computational chemistry, the "smooth" molecular surface [25] is actually defined as the level set of a signed distance function, which strictly speaking is not always smooth. The related singularities have caused trouble when meshing or visualizing molecular surfaces. Some meshing algorithms of molecular surfaces even fail when the surface contains singularities. However, the BPB diagram allows us to compute all surface singularities a priori, which was first introduced in our previous work [15,26]. In the following content, we will briefly present the basic idea.
In an implicit solvation model, the solute molecule is commonly regarded as a union of atomic balls, such as the van der Waals balls with radii r v,i . Meanwhile, the solvent molecules are simply idealized as spherical probes with radius r p . The solvent-excluded surface (SES), denoted by Γ ses , is the boundary of the cavity where solvent probes can not touch (the solute atoms and the solvent probes can not overlap). We take Γ as the boundary of the union of balls with atomic centers and radii r i = r v,i + r p . Then, the SES can be represented mathematically as the following level set: where f Γ is the signed distance defined in Eq. 7. Figure 7 illustrates the signed distance of the benzene molecule in the XY plane. The SES is the so-called smooth molecular surface, which however might have plenty of singularities. The BPB diagram gives a partition of the union of balls, based on the signed distance to different boundary components. As a consequence, given any point x in the union, one can find conveniently its closest point(s) on Γ by determining which BPB cell this point lies in. Note that any point x ∈ Γ ses is a singularity (in the sense that f Γ is nondifferentiable at x) if and only if x has more than one closest points on Γ . However, the number of closest points can be counted directly as soon as the BPB diagram is constructed, which means that all surface singularities can be computed a priori.  Figure 8 illustrates the SES of a protein (generated by our MolSurfComp package [27]), where the green curves highlight singularities. It is well known that Γ ses is composed of three types of patches: convex spherical patches (in red), toroidal patches (in yellow), and concave spherical patches (in blue). These patches correspond respectively to the 2-patches, the 1-patches (circular arcs in yellow), and the The above formula is computable as soon as the BPB cell is given (see [15] for details). In particular, the singularities on each SES-patch lie on the boundary of the BPB cell, because any interior point of the BPB cell has only one closest point on Γ .
Remark 6 Generally speaking, given an arbitrary surface in 3D, the signed distance from one point to the surface is expensive to compute. However, in the special case for balls, the BPB diagram allows to compute this value directly.

Medial axis of molecule
The medial axis of an object is the set of all points having more than one closest point on the object's boundary. According to the definition of BPB cells, we can claim that the medial axis of a molecule is part of the boundaries of BPB cells. The concept of medial axis was first introduced by Blum [28] and was originally referred to as the topological skeleton. It has been shown that the medial axis of an object is always homotopy equivalent to the object itself [29] and the medial axis is useful for shape descriptions. As a consequence, it is potentially a good idea to use the medial axes of molecules for classifying their structures.
In fact, the theory on the BPB diagram allows us to compute the medial axis of 3-balls. As mentioned above, the boundary of the union of 3-balls is composed of kpatches with k = 0, 1, 2. Given a k-patch γ (k) i , let i = I(γ (k) i ) be the indices of the 3-balls generating this k-patch. In the case when k = 1, 2, according to Corollary 1, we can claim that conv(c i ) is part of the medial axis. The reason is that ∀x ∈ conv(c i ), any point on γ (k) i is a closest point of x. In other words, x has infinite number of closest points. In the case when k = 0, given an arbitrary intersection point γ For the sake of illustration, we provide a simple example of the medial axis of molecule 1yjo in Fig. 9. As mentioned above, the molecule and its medial axis are homotopy equivalent.

Conclusion
In this article, we have introduced the so-called boundary-partition-based (BPB) Voronoi diagram of d-balls in the Euclidean space E d , which can be seen as an alternative partition to the power diagram. We have studied the properties of this diagram and presented its applications. Compared to the power diagram, the BPB diagram is more convenient to be implemented for volume computations of balls, avoiding the possibly complicated computation of power cells. In addition, this BPB diagram can be used to compute singularities on molecular surfaces and to compute the molecular medial axis for possible structure classification in the context of protein-ligand binding affinity prediction. At this moment, these applications are restricted to E 2 and E 3 just like the power diagram. We expect nevertheless possible applications in higher dimensions.
Funding Information Open Access funding provided by Projekt DEAL. This work is partially supported by the National Natural Science Foundation of China (Grant No. 11901281).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/.
Since (x + λv t , n s ) ≤ b s , ∀λ ≥ 0, we know that (v t , n s ) ≤ 0. In the case of (v t , n s ) < 0, according to Eq. 56, there exists a small enough number ε i t > 0 such that ∀x ∈ [0, ε i t ], |ζ s (x) − c i t | ≥ |x − c i t | = r i t .
Besides, in the case of (v t , n s ) = 0 implying v t ∈ A s , according to Eq. 28 and (56), we have ∀0 ≤ x ≤ |v|, In particular, when t = t 0 , we have both (v t , n s ) = 0 and (v, v − v t ) = 0. As a consequence, which means that ζ s (x) ∈ S i t 0 . So far, we have proved that ∀x ∈ [0, ε i t ], ζ s (x) does not lie in the interior of S i t and lies on the sphere S i t 0 . Second, we consider a sphere S j ∈ {S i 1 , S i 2 , . . . , S i m } that does not contain x. In this case, we have |x − c j | − r j > 0. Therefore, there exists a small number ε j > 0 such that ∀x ∈ [0, ε j ], This yields that ∀x ∈ [0, ε j ], that is to say, ζ s (x) lies outside the sphere S j . In summary, there exists a possibly small number ε > 0 such that ∀x ∈ [0, ε], ζ s (x) does not cross any sphere S i and lies on the sphere S i t 0 , which implies that ζ s (x) ∈ Γ .
It is sufficient to prove that any point x ∈ ∂ R 0 belongs to either some face F ij or some set F 0 with dim(F 0 ) ≤ d − 2. ∂ R 0 can be divided into two sets U 1 := {x ∈ ∂ R 0 | all closest points of x belong to P 0 }, and U 2 := {x ∈ ∂ R 0 | there exists a closest point of x not contained in P 0 }.
In the following content, we prove that if x ∈ U 1 , then x belongs to a certain face F ij , while if x ∈ U 2 , then x belongs to F 0 which will be defined later.
Step 1: In the case of x ∈ U 1 , we can find a sequence of points {x n } in Ω\R 0 such that x n tends to x. Correspondingly, there exists a sequence of points {a n } on Γ , where a n is one closest point of x n . Since x has finitely many closest points in P 0 and the total number of k-patches is finite, we can extract a subsequence of a n such that this subsequence lies on some k-patch γ (k) with k ≥ 1 and converges to some nondegenerate intersection point γ i . Without loss of generality, we can therefore suppose that a n tends to γ (0) i and a n ∈ γ (k) , ∀n. As a consequence, γ (0) i is on the boundary of γ (k) and further, there exists a 1-patch γ (1) ij on γ (k) , satisfying c I(γ (k) ) ⊆ c I(γ (1) ij ) .
Due to the fact that x n ∈ conv a n , c I(γ (k) ) , we then have by taking n → ∞.
Step 2: In the case of x ∈ U 2 , we want to prove that x belongs to some F 0 .
According to the definition of U 2 , x has at least one closest point a that is not a nondegenerate intersection point. Here, we mention the fact that for any point y belonging to the open line segment ax with endpoints a and x, a is the unique closest point of y on Γ , which can be easily proven by contradiction.
On the one hand, if a is not an intersection point, then a lies on some k-patch γ (k) with k ≥ 1. According to Theorem 1, we know that R(a) = conv a, c I(γ (k) ) . (64) Considering that the latter convex hull, we obtain that x ∈ conv(c I(γ (k) ) ) of dimension dim(conv(c I(γ (k) ) )) ≤ d −2, since otherwise, x will has a unique closest point on Γ . On the other hand, if a is a degenerate intersection point, then we have a ∈ Λ I(a) with dim Λ I(a) ≤ d − 1. Since R(a) ⊆ Λ I(a) according to Theorem 1, it holds that dim (R(a)) ≤ d − 1. Here, we actually have a ∈ R(a) and R(a) is a convex set from Lemma 3. Due to the fact mentioned above, we obtain that x ∈ R(a) only lies on ∂R(a) of dimension dim (∂R(a)) ≤ d − 2.