An efficient algorithm for approximate Voronoi diagram construction on triangulated surfaces

Voronoi diagrams on triangulated surfaces based on the geodesic metric play a key role in many applications of computer graphics. Previous methods of constructing such Voronoi diagrams generally depended on having an exact geodesic metric. However, exact geodesic computation is time-consuming and has high memory usage, limiting wider application of geodesic Voronoi diagrams (GVDs). In order to overcome this issue, instead of using exact methods, we reformulate a graph method based on Steiner point insertion, as an effective way to obtain geodesic distances. Further, since a bisector comprises hyperbolic and line segments, we utilize Apollonius diagrams to encode complicated structures, enabling Voronoi diagrams to encode a medial-axis surface for a dense set of boundary samples. Based on these strategies, we present an approximation algorithm for efficient Voronoi diagram construction on triangulated surfaces. We also suggest a measure for evaluating similarity of our results to the exact GVD. Although our GVD results are constructed using approximate geodesic distances, we can get GVD results similar to exact results by inserting Steiner points on triangle edges. Experimental results on many 3D models indicate the improved speed and memory requirements compared to previous leading methods.


Introduction
Voronoi diagrams are fundamental data structures that have been widely investigated in computational geometry, and applied in various science and engineering disciplines, e.g., molecular biology [1], image processing [2,3], computer graphics [4][5][6], machining and manufacturing [7], and mobile communication [8], etc. A Voronoi diagram can be viewed as the minimization diagram of a finite ensemble of continuous functions [9]. A Voronoi diagram divides the embedding space into sub-regions, each being the domain closer to a given seed point than the others. We may define many variants of Voronoi diagrams depending on the embedding space. In Euclidean space, the Voronoi diagram is well understood and studied for rasterization data [10,11]. For spaces with non-Euclidean metrics, people have developed Voronoi diagrams on parametric surfaces [12], polyhedral surfaces [13,14], spherical spaces [15], hyperbolic spaces [16], etc. However, there are distinct differences between Euclidean and non-Euclidean metrics, so many established attributes based on the Euclidean metric do not hold for non-Euclidean metrics. For instance, a 2D Euclidean Voronoi component is always convex, while a geodesic Voronoi component is often non-convex. In this paper, we specifically consider Voronoi diagrams embedded in triangular meshes, based on the geodesic metric (geodesic Voronoi diagrams, GVDs).
Owing to its simplicity and flexibility, the triangular mesh is the most frequently used 3D representation in computer graphics [17]. Voronoi diagrams on triangular meshes have been used in a range of applications, such as 3D mesh reconstruction [18,19], tree skeleton extraction [13], trajectory planning [20], shape segmentation [21], and so on. However, making use of geodesic Voronoi diagrams on triangular surfaces, rather than diagrams based on the Euclidean metric, presents many unresolved challenges. Often, in order to construct geodesic Voronoi diagrams, exact geodesic algorithms are used, such as the MMP algorithm [22] and the VTP algorithm [23], to obtain the geodesic metric [13,14,24]. Due to the excellent features of these algorithms, they can supply adequate geodesic information about edges for GVD construction. However, since these exact algorithms are slow and memory hungry, they present a stumbling block to constructing GVDs and further practical applications. Furthermore, the boundaries of GVD regions, bisectors, consist of hyperbolic and line segments that are have particular, more complex characteristics than those of Euclidean Voronoi diagrams.
In this paper, we focus on avoiding superfluous computation of geodesic distances to improve speed and memory usage. As computing geodesic distances using exact algorithms is expensive, we construct graphs using Steiner points placed on mesh edges to report distances. In the propagation process, we regard each Steiner point as a simplified window and adapt existing window filtering rules to remove many redundant events. Theoretically, a bisector in a Voronoi diagram on a mesh may include linear, hyperbolic, and parabolic segments. We make use of the Apollonius diagram, also called the additively weighted Voronoi diagram, to depict these structures locally defined in the plane of a single mesh triangle.
Based on these ideas, we present an efficient method for constructing an approximate GVD for point-source generators. Once distance propagation terminates, geodesic distances and the GVD become concurrently available. Thanks to the highly selective filtering rules, our algorithm runs in O(mnlogn) time and takes O(mn) space on an n-face mesh with t generators, where m is the number of Steiner points inserted on each edge. Figure 1 is a representative of experimental results on a large number of 3D models, indicating that our algorithm has improved speed and memory efficiency compared over previous leading methods.
In summary, our key ideas include the following: instead of using exact methods, to obtain geodesic distances, increasing speed and reducing memory usage. (2) We utilize local Apollonius diagrams with weighted Steiner points to partition mesh triangles and extract hyperbolic and linear segments, to encode the complicated GVD bisector structures. (3) We suggest a method of approximate GVD evaluation based on discrete Fréchet distance, and use it to show that our approximate GVD results are close to the exact results. Our work is presented in detail as follows. Section 2 briefly recalls related work. Section 3 provides background on the Apollonius diagram and geodesic distance computation based on Steiner points. Section 4 explains how to construct a GVD using point-source sites. We explain our algorithm for building GVD on triangle meshes in detail in Section 5. Section 6 experimentally compares our method to previous methods. We give some applications of the proposed GVD algorithm in Section 7.

Discrete geodesics
The discrete geodesic problem (DGP) focuses on computing geodesic distances on triangle meshes [22]. There are three main kinds of methods for solving it: discrete wavefront propagation methods, PDE methods, and graph-based methods. We refer to Refs. [25,26] for a comprehensive survey. Discrete wavefront propagation methods use a window to encode geodesic paths sharing the same edge sequence. They generally propagate wavefronts across triangle faces in a Dijkstra-like sweep. The MMP algorithm [22,27] and the CH algorithm [28] are two representative methods. The MMP algorithm maintains a partitioning scheme for each directed edge while the CH algorithm uses a one angle one split rule so that at most O(n) windows are allowed to have two children. Many variants and improvements [23,[29][30][31][32] propose different rules to filter redundant windows in the propagation process. Although discrete wavefront propagation methods can provide exact geodesic information and adequate information for bisectors, these methods are computationally expensive.
PDE methods [33][34][35][36] compute geodesic distances using partial differential equations (PDEs) on a smooth manifold. In detail, they solve the Eikonal equation ∇u(t) = 1 with boundary condition u(s) = 0 on discrete domains, where s is the source. PDE methods are easy to implement and efficient. However, they provide only a first-order approximation, which is extremely sensitive to the mesh tessellation.
Graph-based algorithms [37][38][39][40][41][42][43][44] are approximation algorithms; they can balance accuracy and speed very well. These algorithms are characterized by an approximation ratio . A typical goal is to compute a (1 + )-approximate geodesic path as quickly as possible. Graph-based algorithms usually deal with the discrete geodesic problem by constructing a dense undirected graph G on the target manifold mesh and use graph-based search techniques to answer the distance-and-path query. In general, approaches add Steiner points on mesh edges or even on triangle faces. For example, Lanthier et al. [39,40] constructed a dense graph by inserting Steiner points on each mesh edge to achieve accurate results. Aleksandrov et al. [41,45] proposed adding Steiner points more densely near vertices than near mid-edges to ensure that the length of the approximate path is within a factor (1 + ) of the shortest geodesic path. Later, they [42] found that adding Steiner points along the bisectors of triangles also works. A survey on graph-based methods is presented in Ref. [25]. Generally, the Steiner-point insertion algorithms need to construct a graph G in which the complexity of G grows quadratically with respect to the number m of Steiner points inserted on each mesh edge. Consequently, the computational cost grows rapidly as m increases. Since efficiently computing geodesic is crucial in GVD computation, we reformulate a graph method [46] based on Steiner point insertion; we regard each Steiner point as a simplified window and use existing window filtering rules for efficiently propagating geodesic distances.

Voronoi diagram
Mathematically, a Voronoi diagram is a partition of the plane into regions, each closer to one of a given set of objects. Using the same definition, the construction of Voronoi diagrams can be extended to metric spaces other than Euclidean spaces [47,48]. For non-Euclidean spaces, there are many works dedicated to Voronoi diagrams on smooth surfaces, e.g., Riemannian manifolds M [49,50], parametric surfaces [12], spheres S 2 [15,51], and hyperbolic spaces H 2 [16]. The reader is referred to Refs. [10,52] for detailed surveys.
Owing to their simplicity and a high degree of freedom, non-differentiable polyhedral surfaces, especially triangular surfaces, are the most popular representation for 3D objects. Therefore, constructing a Voronoi diagram on triangle mesh surfaces using the geodesic metric is an important topic in computer graphics. Since many properties of a smooth manifold no longer hold, transferring Voronoi diagrams from Euclidean metric to geodesic metric still presents many unresolved challenges. Kimmel and Sethian [53] computed Voronoi diagrams on triangle meshes using the fast marching method [33]. However, since that method computes only a first-order approximation to geodesic distance, it may produce bad results on poor triangulations. Liu et al. [13] studied the analytic structure of iso-contours, bisectors, and Voronoi diagrams on a triangular mesh M . They proposed practical algorithms for constructing GVDs using the MMP algorithm for exact geodesic computation. Later, they extended their work to construct intrinsic Delaunay triangulations from dual geodesic Voronoi diagrams [54,55]. Xu et al. [24] proposed an algorithm for constructing GVDs with polyline generators. They introduced a new concept, the local Voronoi diagram (LVD), to represent the GVD bisectors. Recently, Qin et al. [14] proposed constructing Voronoi diagrams using the window-VTP algorithm which runs 3-8 times faster than Liu et al.'s method [13]. Although exact geodesic algorithms like MMP and VTP can compute Voronoi diagrams accurately, obtaining exact geodesic distances is time-consuming, presenting a stumbling block to constructing GVDs and further practical applications. In this paper, we employ an improved Steiner-point graph method to balance time and accuracy of GVD construction.

Steiner points
In computational geometry, a Steiner point is a point that is not part of the input to a geometric problem but is added during the solution of the problem. It usually creates a better solution than using the original points alone [56]. In geodesic computation, Steiner points are regularly applied in graph-based methods [39,40,42] to balance accuracy and time. By inserting m Steiner points into each mesh edge, and connecting any pair of Steiner points in the same triangle, we can construct an augmented graph G. In the original Steiner-point graph algorithms [39][40][41][42], when a Steiner point p is assigned a shorter distance, it must broadcast an update event to all of the neighboring Steiner points and vertices: see In order to seek a balance between accuracy and speed in such graph algorithms, Meng et al. [46] proposed an improved Steiner-point graph method. In it, the inserted Steiner points play much like the same role as windows in exact geodesic methods. We can filter useless broadcast events by adapting existing window pruning rules. In order to efficiently construct Voronoi diagrams based on the geodesic metric, we adopt filtering strategies and revise Meng et al.'s original method for computation of the multisource geodesic distance field.

Apollonius diagram
The Voronoi diagram is a fundamental data structure defined as the minimization diagram of a finite set of continuous functions [10,52,57]. It subdivides the embedding space Ω ⊂ R d into certain regions based on the Euclidean distance to a set of points , called a site or generator, dominates the subregion: This is called the Voronoi cell of the site x i . We may define many variants of Voronoi diagrams using different distance functions and embedding spaces. Power diagrams and Apollonius diagrams are two classical variants of Voronoi diagrams (see Fig. 3).
A power diagram is a type of weighted Voronoi diagram. Instead of each region comprising the points closest to a site, it comprises those points with smallest power distance for a particular circle [58,59].
Each site x i with weight w i dominates the subregion or power cell: A power diagram degenerates to a Voronoi diagram when the weights are identical. The Apollonius diagram, or additively weighted Voronoi diagram, is another variant [60][61][62]. Differing from the power diagram which uses squared distance, each site x i with weight w i dominates the subregion, or Apollonius cell: The sets of points simultaneously belonging to two Apollonius cells are called Apollonius edges. The Apollonius edges between two adjacent cells usually contain linear and hyperbolic segments, and is the motivation for applying Apollonius diagrams rather than other Voronoi diagram variants to construct the GVD.

GVD with point-source sites
Given a set of sites S = {s 1 , · · · , s t } on mesh M, let G s i (p) denote the geodesic distance from site s i to point p on M. The geodesic Voronoi cell associated with site s i can be defined as Furthermore [14], the boundaries of geodesic Voronoi cells are given by the sets of points q satisfying: The common boundary of sites s i and s j can also be regarded as their bisector.  Proof. The key is to analyze the bisectors as they appear in a triangle F. In the Steiner-point insertion method, m 0, Steiner points are inserted into each triangle edge. There are two cases to consider according to whether the Steiner points used as relays to compute distance to p have equal or unequal distances to their respective sites: • Case 1: when G s 1 (q 1 ) = G s 2 (q 2 ), point p lies on the line segment bisecting sites s 1 and s 2 . • Case 2: when G s 1 (q 1 ) = G s 2 (q 2 ), point p satisfies G s 1 (q 1 ) + q 1 p = G s 2 (q 2 ) + q 2 p , which implies that p is on a hyperbolic segment with foci q 1 and q 2 .

Bisectors
If three or more Steiner points give the same distance to point p, it is a branch point. This completes the proof. Figure 4 illustrates the above situations.
In order to find the boundaries of GVDs, we employ Apollonius diagrams to guide the capture of these complicated structures, as we next explain.

Local Apollonius diagrams
Unlike traditional Voronoi diagrams, Apollonius diagrams consist of linear and hyperbolic segments.
Taking (q 1 , −G s 1 (q 1 )) and (q 2 , −G s 2 (q 2 )) as two weighted points, their additively weighted distances to p are respectively Thus, point p has an equal-weighted distance to sites s 1 and s 2 via Steiner points q 1 and q 2 , respectively. More generally, bisectors B (may be with a treelike structure) can be reported by computing the Apollonius diagram with respect to a set of Steiner points on the boundary of triangle F. We may represent a triangle by its vertices: the set of Steiner points inserted into the edges of F. When the Steiner-point graph method terminates, each Steiner point has a geodesic distance G s i (q k ) measured from generator s i , i = 1, · · · , t. Taking the set of weighted points where M cannot exceed the total number of inserted Steiner points. The boundaries of regions are enclosed by line segments and hyperbolic segments. Each segment lies on the border of two distinct regions belonging to different generators, and these form the boundaries of the GVD. See Fig. 5 for an example. Therefore, using the strategy above, by constructing the local Apollonius diagram we can extract the GVD structure triangle by triangle on mesh M.

Redundant primitives
While this approach to constructing the geodesic Voronoi diagram on a mesh M could be used in a direct way by first computing the Apollonius diagram on each mesh face, and then finding diagram edges that belong to the GVD, this is time-consuming, since typically only a handful of mesh triangles include GVD boundaries. To develop an efficient method to construct the GVD of mesh M, it is critical to distinguish those triangles which include the GVD boundaries [14], instead of conducting a brute force search over all triangles. Figure 6 shows a point p which is the intersection between a triangle edge and a GVD boundary. Point p must satisfy the constraint in Eq. (5) and is shared by two adjacent Steiner points connected to two different sources. We call points like p critical points; triangles including GVD boundaries always contain such points. In other words, any triangle incorporating GVD structures has Steiner points propagated from different sources, giving criteria for determining the redundant primitives on a mesh: Given a triangle mesh M = (V, E) and distinct     Using these definitions to discard redundant primitives, we need only process a few remaining triangles of mesh M to efficiently extract the geodesic Voronoi diagram.

Basis
Based on the ideas in Section 4, we present an efficient and practical algorithm to construct geodesic Voronoi diagrams for point source sites. The input is a triangular mesh M and a set of sites S = {s 1 , · · · , s t }. The output is the approximate GVD (AGVD) for M. To provide the multi-source geodesic distance field used in constructing the GVD, we modify Meng et al.'s method [46]. We use a priority queue Q to organize the Steiner points, which represent wavefronts in the propagation process.

GVD edge extraction
We discard redundant primitives during distance propagation in the Steiner-point graph method. During each distance iteration, the geodesic computation method propagates the distances of Steiner points using the least geodesic distance. Therefore, the distance of the Steiner point at the head of the priority queue Q does not decrease. If the edges of a non-discarded triangle ABC have both distances shorter than the distance of the point at the head of the queue, we can conclude that ABC is behind the geodesic wavefront; all Steiner points behind the wavefront are no longer updated. We say such triangles are inactive. Thus, we check whether or not ABC contains GVD edges.
In exact geodesic algorithms, there are O(n) windows on each edge of ABC. Computation of the exact distance from the source to each edge is time-consuming. In our practical implementation, we use an upper bound of the distance from the source s to each edge e [24]. The upper bound distance d(s, e) is (d(a)+d(b)+ e )/2, where d(a), d(b) are the geodesic distances from the source to the endpoints of edge e and e is its length.
Frequently checking for inactive triangles can extract GVD edges in early stages and consume lower memory. However, it has a significant influence on time. In order to balance memory use and time, we check for inactive triangles every λn (an integer) distance iterations, where n is the number of triangles in M: smaller λ results in more frequent checks. After every λn distance iterations, we construct the local Apollonius diagram and extract GVD edges from non-discarded triangles in the inactive region.

GVD construction
For the Steiner points on opposite edges of site s i in its 1-ring neighbourhood, we initially assign distances and push them into a priority queue Q. Then, our algorithm progressively computes geodesic distances and extracts GVD structures by performing the following steps in each distance iteration.
(1) Take the Steiner point p from the head of Q and propagate from p to Steiner points on its adjacent triangles. (2) Find the adjacent Steiner points which have shorter distances via p and insert them into the priority queue Q. (3) If the number of distance iterations is an integral multiple of λn or the priority queue Q is empty, find non-discarded triangles to construct the local Apollonius diagram, and extract GVD edges. In order to robustly extract GVD edges, besides Steiner points on each retained (non-discarded) triangle, the surrounding inserted Steiner points of target triangles are also added into the construction process to build joint local Apollonius diagrams, i.e., Steiner points in the 1-ring neighbourhood of retained triangle F are employed as relay points to partition and extract GVD edges from F. Algorithm 1 gives further details of GVD construction.

Complexity analysis
We now analyze the time and space complexity of our algorithm to demonstrate its theoretical efficiency. Previous GVD methods [13,14,24] based on the exact geodesic metric require O(n 2 log n) time and O(n 2 ) space to construct the GVD. Our method has an obvious theoretical advantage in both time and space. Although our GVD is constructed using approximate geodesic distances, we can get a result close to the exact result by inserting Steiner points into triangle edges. We evaluate our algorithm in Section 6.

Setting
We implemented our algorithm in Microsoft Visual C++ 2015 without using additional numeric packages. All experiments were conducted on a computer with a 3.60 GHz Intel i7-9700K CPU and 8 GB memory. We evaluated our method on both 3D models widely used in the graphics community, and the Thingi10k 3D shape repository [64], which contains a large number of man-made anisotropic meshes and some have an extremely high degree of anisotropy.

Error evaluation
Our GVD results are based on approximate geodesic distances, and we can balance accuracy and time required by varying the number of Steiner points inserted into triangle edges. Here, we employ Fréchet distance to assess error in the constructed GVD; it is widely used to measure the similarity of polylines [65,66]. This distance is also known as dog-leash distance, as it can be intuitively explained in terms of the shortest leash needed when a person and their dog walk at varying speeds along the respective curves.
Converting the curves into polylines is a typical approach [67] for computing the Fréchet distance between arbitrary curves. The method is based on all pairings between the endpoints of line segments on each of the polygonal curves, and uses the paired distances to approximate Fréchet distance, as follows. Let polygonal curve C in R 3 be a continuous function: [0, n] → R 3 . Let σ(C) denote the sequence (C(0), · · · , C(n)) of endpoints on polygonal curve C. Let P and Q be polygonal curves, with corresponding sequences σ(P ) = (u 1 , · · · , u p ) and σ(Q) = (v 1 , · · · , v q ) . A coupling L between P and Q is a sequence (u a 1 , v b 1 ), · · · , (u a m , v b m ) of distinct pairs from σ(P ) × σ(Q). The length L of the coupling L is the length of the longest link in L: The discrete Fréchet distance between curves P and Q is then defined as F(P, Q) = min{ L | L is a coupling between P and Q} Let us consider the GVD on a genus-r (r 0) mesh M with sites S = {s 1 , · · · , s t }, t 2, and assume each cell GC(s i ) consists of a set of closed curves C = {c 1 , · · · , c k }, we define the Fréchet distance of cell GC(s i ) to be k j=1 F(c j , c j ) where c j is the curve for the exact GVD Voronoi cell and c j is the curve of our GVD Voronoi cell. We further define the Fréchet distance between the exact GV D and our GV D on mesh M with t sites as In order to evaluate the approximate GVD results more intuitively, we define the error function δ( GV D) as This error measure depends on the choice of the constant c. If it is chosen too small or too large, the function poorly describes GVD error. Our extensive experiments show that the best choice of c ≈ 10.
Obviously, δ( GV D) is in [0, 1] where 0 indicates that our GVD and the exact GVD result are identical and 1 means they are completely different. In order to construct GVD, we use the Steiner point insertion method [46] to acquire the geodesic distance field. Errors in geodesic distances are controlled by inserted Steiner points on triangle edges which indirectly influence the accuracy of the AGVD construction. Intuitively, the more inserted Steiner points on each edge, the higher expected accuracy of the AGVD result. Our method can balance the time, memory, and accuracy well by varying the number of inserted Steiner points.
Taking the 50k-face Cow model as an example in Fig. 9, we place 30 generators on the model and insert 5, 8, and 12 Steiner points on each edge to construct the AGVD. As Fig. 9 shows, the AGVD errors rapidly decrease with increasing Steiner points. With 12 Steiner points, our algorithm only takes 0.24 s and 1.32 MB memory to construct an AGVD with 0.07% error, while the exact GVD method [14] takes 0.59 s and 2.79 MB memory. Generally, we set m = 8 as the default number of Steiner points on each edge in our experimental tests.
Our algorithm uses relatively little time and memory to construct an approximate GVD with an extremely low error close to the exact GVD result. This has an obvious advantage in error insensitive applications, e.g., geodesic remeshing, which we consider in Section 7. Without a doubt, higher accuracy must require more computational cost. However, our algorithm can achieve much higher accuracy at the cost of a small extra computational amount. As Fig. 10 shows, when m is set to 10, the time for GVD construction is 0.51 s, and when m increases to 30, the time used is only 2.26 s. Table 1 gives further experimental results to emphasize the balance between time, memory, and accuracy.

Assessment
We use running time and peak memory to compare performance between our method and other leading methods. There are many factors which influence the performance of GVD construction methods. Here, we select the three most significant quantities: number m of Steiner points per triangle edge, number t of source sites, control parameter λ balancing time and memory. Since choice of number m of Steiner points was discussed in Section 6.2, we concentrate on the other two factors here.

Number of source sites
We first consider the effect of increasing the number of source sites, using various models. We used 2- 2000 source sites and measured AVGD time and memory consumption. Results are shown in Fig. 11 and Table 2.
As Table 2 shows, the time needed for GVD construction drops with an increasing number of generators, since each geodesic wavefront just covers a small region of the model surface. Furthermore, more time is needed because a relatively large t implies that  the local Apollonius diagram is more complicated, so takes longer to compute.

Balancing time and memory
Our algorithm considers inactive triangles after every λn distance iterations where λ controls the balance between time and memory. Smaller λ results in more frequently checking the inactive triangles, so extracting GVD edges earlier and consuming less memory. However, this has a significant effect on time taken. Using a larger λ results in less frequent checking of the inactive triangles, which is faster but consumes more memory. We verify these claims and show the effect of changing λ in Fig. 12. We set λ = 1.0 by default value in our experimental tests.

Fig. 12
Balancing time and memory, using λ. We employ the 28k-face Kids model to construct a GVD with 50 sites and vary λ to control the balance between time taken and memory usage.

Qualitative comparison
In this section, we make a qualitative comparison between our method and the four typical methods most in common with our method.

Comparison to Ref. [13]
Liu et al. studied the analytic structure of bisectors and Voronoi diagrams on triangulated surfaces, and proposed a practical algorithm to construct the Voronoi diagram on a triangulated mesh M. They used the MMP method [22] to obtain geodesic information for GVD construction. Since the distance field on an edge e of M can have O(n) extrema, they partition e into sub-edges such that the distance field value on each sub-edge is monotone and linear. Compared to Liu et al.'s method, our method has two significant advantages. On one hand, their method needs to subdivide mesh triangles until the bisectors cross each face at most once. As a result, the method is sensitive to the mesh resolution while our method is robust on meshes of the same shape with different resolutions. See Section 6.5 for further details. On the other hand, using the Steiner-point graph method, we can get GVD results which are close to exact, using less time and memory. See Table 1.

Comparison to Ref. [24]
Unlike most previous GVD algorithms, this approach uses polyline generators. They revealed that a typical GVD bisector contains linear, hyperbolic, and parabolic segments. To overcome the challenge this presents, a new concept called the local Voronoi diagram (LVD) was introduced, which is a combination of additively weighted Voronoi diagrams and line-segment Voronoi diagrams defined locally in the plane. Our algorithm extracts GVD edges in a similar manner to this method. The noticeable difference is that they use the MMP framework to compute exact geodesic distances, while our algorithm relies on approximate geodesics computed by the Steiner point insertion method. Our relative advantage is that we can flexibly balance the time, memory, and error during GVD construction; this is an obvious advantage in error insensitive applications, as shown in Section 7.

Comparison to Ref. [14]
This paper proposed an efficient algorithm to construct the GVD, aiming to reduce redundant computation to save time and memory. A redundant window removal process is used during GVD construction: in the window-VTP algorithm it selectively retains windows on edges. The key issue is how to detect and remove redundant windows simultaneously with geodesic wavefront propagation. Compared with this method, our method is easier to implement and more computationally efficient. Using the 50k-face Monster model in a test, our method is about three times faster for an error of less than 0.3% (see Fig. 1). Meanwhile, our peak memory usage is about one third of that of Qin et al.'s method. 6.4.4 Comparison to Ref. [37] In this paper, Xin et al. presented an approximate algorithm to construct the GVD using point generators. They first computed multisource geodesic distances using generators on the target mesh. Upon termination of this phase, each vertex is labelled with its nearest source. They then construct bisectors by checking triangles containing at least two intersection points. If there are two intersection points, a segment is used to connect points. If there are three intersection points, a point is found inside the triangle and connected to all three intersections. Finally, the GVD is constructed by tracing the bisectors. The method is easy to implement and works well for high quality and resolution meshes. However, it may produce low quality results when the input meshes are poor, e.g., have irregular triangulations or sharp corners. The natural solution is to improve the triangulation quality or increase mesh resolution. However, remeshing and subdivision operations are time-consuming and may have their own issues. Compared to this method, our algorithm is robust to mesh resolution and triangulation quality, as well as providing an efficient way to construct an approximate GVD. We next discuss robustness of our method further.

Robustness
Our algorithm is robust to changes in mesh resolution and triangulation quality. In order to demonstrate robustness, we constructed the AGVD on a triangular mesh with different resolutions and varying quality. The Lo value [68] provides a general way to evaluate mesh quality; higher quality meshes have values closer to 1. AGVD results on the Lucy model at three different resolutions and three different qualities are shown in Fig. 13, indicating that our algorithm can construct reliable and relatively consistent AGVD results. We also examined our method under varying mesh anisotropy. We use the anisotropy measure defined in Xu et al. [31] for a triangular face f : where P , H, and S are respectively its half-perimeter, longest edge length, and area. The mesh anisotropy measure is then computed as the average of τ (f ) over all faces of the mesh. Larger τ indicates a higher degree of anisotropy. We selected two models from the Thingi10k repository with average anisotropy τ avg of 31.1 and 9.9, and maximal anisotropy τ max of 1062.3 and 25,183.7, respectively. Figure 14 shows our AGVD results for 50 sites on these highly anisotropic meshes. Fig. 14 Our algorithm produces accurate GVD results for meshes with highly anisotropic triangles (average anisotropy τavg = 31.1, 9.9, respectively). Using 12 Steiner points per edge, AGVD errors are 0.42% and 0.38%, demonstrating that our method copes well with anisotropic triangles.

Application to remeshing
That the Delaunay triangulation of a point set S is the dual of its Voronoi diagram is well known. Leibon and Letscher [49] showed that, if the sampling points are sufficiently dense, the dual triangulation of the Voronoi diagram on M exists. This observation has an application to remeshing the compact triangular models reconstructed from laser scanning data. In detail, we first compute a geodesic Voronoi diagram (GVD) by selecting sources as generators and then find its dual graph as the remeshed modelM. In this process, we can save time and memory by constructing an AGVD close to the exact result, making a new mesh more efficiently than related methods [13,14]. As a demonstration, we remesh the 200k-face Centaur model using 2000 randomly selected sources. Liu

Conclusions and future work
In this paper, we have presented an efficient and practical algorithm to construct geodesic Voronoi diagrams with point source sites. Unlike previous methods that depend on the exact geodesic metric to construct GVD, we reformulate a Steiner point insertion method [46], to obtain geodesic distances effectively. This substitution reduces GVD computation time and memory usage. We utilize local Apollonius diagrams with weighted Steiner points to partition mesh triangles and extract hyperbolic and linear segments for encoding the complicated GVD bisector structures. In order to balance memory and time requirements, we only process significant triangles instead of using a brute force search over all triangles. Every λn distance iterations, where n is the number of triangles in M, we construct the Apollonius diagram and extract GVD edges. We also suggest an evaluation measure based on discrete Fréchet distance to assess similarity between our result and the exact GVD result. Although our GVD results are based on approximate geodesic distances, we get GVD results close to exact results by inserting Steiner points on triangle edges.
In the future, we hope to enhance our algorithm so that it iteratively inserts Steiner points until a pre-specified error bound is reached. We also intend to investigate other applications that would benefit from approximate geodesic Voronoi diagrams, such as tree skeleton extraction and classification, point pattern analysis on a mesh [13], and so on.
R&D Program of China (2021YFB1715900), and the Joint Funds of the National Natural Science Foundation of China (U22A2033).

Declaration of competing interest
The authors have no competing interests to declare that are relevant to the content of this article.