Geodesic Hermite Spline Curve on Triangular Meshes

Curves on a polygonal mesh are quite useful for geometric modeling and processing such as mesh-cutting and segmentation. In this paper, an effective method for constructing C1 piecewise cubic curves on a triangular mesh M while interpolating the given mesh points is presented. The conventional Hermite interpolation method is extended such that the generated curve lies on M. For this, a geodesic vector is defined as a straightest geodesic with symmetric property on edge intersections and mesh vertices, and the related geodesic operations between points and vectors on M are defined. By combining cubic Hermite interpolation and newly devised geodesic operations, a geodesic Hermite spline curve is constructed on a triangular mesh. The method follows the basic steps of the conventional Hermite interpolation process, except that the operations between the points and vectors are replaced with the geodesic. The effectiveness of the method is demonstrated by designing several sophisticated curves on triangular meshes and applying them to various applications, such as mesh-cutting, segmentation, and simulation.


Introduction
Curves on a polygonal mesh play a fundamental role in geometric modeling and mesh processing [1,2], such as remeshing [3], mesh-cutting [4], and mesh segmentation [5]. However, little attention has been paid to the curves on a polygonal mesh, compared with the freeform curves [6,7] in Euclidean space, and only a few methods [8,9] have been developed to generate these curves.
Martínez et al. [8] introduced a geodesic Bézier curve on a triangular mesh by applying the de Casteljau algorithm to the geodesic path between control points on a triangular mesh. The generated curve lies exactly on a triangular mesh and approximates the shape of the control polygon formed by the geodesics connecting consecutive control points. They also presented a method for joining two geodesic Bézier curves to obtain a piecewise Bézier spline curve while maintaining C 1 continuity.
Some applications, such as mesh-cutting and segmentation, require an interpolation curve for given mesh points rather than an approximation curve. For these applications, it is not suitable to use the geodesic Bézier curve introduced in a previous report [8]. In Euclidean space, it is quite easy to construct an interpolation curve, and a variety of methods [6,7] have been developed for their purposes. However, it would be challenging to define an interpolation curve on a triangular mesh consistently.
The Hermite spline is a fundamental technique for data interpolation. For two given points p 0 , p 1 and their tangent vectors v 0 , v 1 , the cubic Hermite curve C(t) satisfying C(0) = p 0 , C(1) = p 1 and C (0) = v 0 , C (1) = v 1 can be represented by the following Bézier curve: where B 3 i (t) = ( 3 i )(1 − t) 3−i t i is the ith Bernstein polynomial of degree 3, defined on t ∈ [0, 1]. There exist several considerations when applying the Hermite interpolation scheme to a triangular mesh. The above equation involves the addition of a tangent vector to a point, which determines the intermediate control points. In general, these control points are not on a triangular mesh; thus, the generated interpolation curve also is not. Figure 1 shows the Hermite spline curve (in black) that interpolates four mesh points (in green). The Hermite spline curve is constructed by (1) without considering the Bunny model and thus they do not lie on the model; the curve segment between second and third points penetrates the model. A simple approach based on projection can be employed to ensure that the generated curve is laid on a triangular mesh. In Figure 1, the pink and yellow curves show the results of projecting the Hermite spline curve (in black) onto the Bunny model in two different viewing directions. Depending on the projection direction, different curves are obtained, and an invalid projection may occur when there is no intersection between the curve and the model. Consequently, this approach is not suitable for constructing an interpolation curve on a triangular mesh.
To resolve this limitation and construct consistent interpolation curves, it is necessary to define appropriate operations between points and vectors on a triangular mesh such as the addition of a vector to a point in (1). For this, the concept of the straightest geodesic that is a straight line from a mesh point in a given direction is employed [10]. The shortest geodesic between two points is approximated by the straightest geodesic that defines a geodesic vector on a triangular mesh.
In this paper, a simple and effective method is presented for constructing a Hermite spline curve that interpolates given points on a triangle mesh by combining previous methods [7,8,10]. The blue curve in Figure 1 shows the geodesic Hermite spline curve generated using the proposed method. The resulting curve is smooth and guaranteed to be on a triangular mesh while interpolating the given mesh points (in green). The main contributions of the method are summarized as follows: • Geometric operations between points and vectors in Euclidean space are extended to those on a triangular mesh by leveraging the shortest geodesic with the straightest geodesic. These operations play an important role in constructing an interpolation curve on a triangular mesh. • The method is simple to implement and easy to use because the conventional Hermite interpolation scheme is followed, except that a tangent vector in (1) is replaced with a geodesic vector and the related operations are defined on a triangular mesh.
• The method generates a geodesic Hermite spline curve that lies exactly on a triangular mesh while interpolating the given mesh points. To the best of our knowledge, our method is a new approach to solving this problem and compared with common projection methods, the proposed method guarantees consistent results and robustness. Moreover, the generated curve exhibits a symmetry property, i.e., the same interpolation curve is obtained if the mesh points are specified in the reverse order.
The remainder of this paper is organized as follows. In Section 2, some related work on the shortest geodesic and the straightest geodesic on a triangular mesh is reviewed, and the appropriate geometric operations needed for constructing Hermite spline curve on a triangular mesh are introduced in Section 3. In Section 4, the construction of Hermite spline curve by combining the conventional Hermite interpolation scheme with new geodesic operations is explained. In Section 5, the effectiveness of the proposed method is demonstrated by showing several examples of interpolation curves on triangular meshes. Finally, the conclusions and suggested future research are provided in Section 6.

Related Work
The aim is to construct piecewise smooth curves on a triangular mesh while interpolating the given points on the mesh. The simplest ones are the connected polylines on the mesh passing through the given points. If one wants to find effective ones with minimum length, the polylines become the shortest geodesics. In this section, previous work directly related to the focus of this paper is reviewed, including geodesic computation and its extension to the geodesic Bézier curve on a triangular mesh.
The shortest path between a pair of points in a plane is the line segment connecting them, but the algorithm must be extended to a Riemann manifold to obtain the shortest path on the three-dimensional mesh [9]. Many geometric problems can be solved efficiently by determining the shortest geodesic between a pair of points on the mesh. Therefore, much research [11][12][13][14] has been conducted on this topic, increasing the efficiency and accuracy of the computation algorithms.
Kimmel et al. [12] computed the shortest path from a triangular mesh by extending fast marching method to the triangular domain. The Mitchell-Mount-Papadimitriou (MMP) [13] and Chen-Han (CH) [11] algorithms are representative methods for computing the shortest geodesic on a polygonal mesh. Both are based on the subdivision of polyhedral surfaces and the continuous Dijkstra algorithm, with the MMP algorithm having O(n 2 logn) and the CH algorithm having O(n 2 ) of time complexity. Surazhsky et al. [15] implemented the MMP algorithm and presented a faster way to find solutions by approximating the geodesic path. Xin et al. [16] significantly reduced the time and space complexity of the CH algorithm by eliminating unnecessary windows using current estimates of distances to the vertex and maintaining priority queues, as with the Dijkstra algorithm.
Polthier et al. [10,17] introduced a new type of geodesic on a polyhedral surface, called a straightest geodesic against the concept of the local shortest geodesic. The straightest geodesic is a straight line from a mesh point in a given direction. Computing the straightest geodesic can be considered an initial value problem with a unique solution, whereas the shortest geodesic solves the boundary value problem of finding the shortest path between two endpoints. If the straightest geodesic only passes through the edges, it is equivalent to the shortest geodesic. The straightest geodesic finds an optimal direction that minimizes the geodesic curvature when it meets a mesh vertex. The computation of the straightest geodesic is detailed in Section 3.
There have been many other studies on geodesic computation on a polygonal mesh [18][19][20][21]. Cheng et al. [18] constructed a smooth surface that approximates a polygonal mesh and computed a geodesic curve on the surface by solving a first-order ordinary differential equation of tangent vector. The discrete geodesic is then obtained by projecting the geodesic curve on the smooth surface onto the polygonal mesh. Lawonn et al. [19] proposed a method for smoothing polylines on a triangular mesh based on local reduction in geodesic curvature and made it possible for users to adjust their proximity to be the straightest geodesic. Qin et al. [20] proposed scenarios for effective window propagation and pruning, and developed triangle-oriented region growing techniques to reduce computational costs and memory usage significantly. Sharp et al. [21] presented a technique for calculating a geodesic on polyhedral surfaces by repeatedly performing edge flips based on the Delaunay flip algorithm [22,23]. As previous study showed [24], research on the shortest and straightest geodesic has continued until recently. In this paper, the straightest geodesic is used to represent a vector on a triangular mesh, and then the related operations between points and vectors on a triangular mesh are extended. Based on these operations, smooth interpolation curves on a triangular mesh can easily be obtained using the conventional Hermite interpolation scheme without significant modification.
Because the shortest geodesics show only C 0 continuity at junction points, they are not sufficient for smoothly interpolating points on a triangular mesh. To overcome this limitation and generate a smooth curve on a manifold or triangular mesh, Park et al. [9] generalized the de Casteljau algorithm in SE(3) (special Euclidean group) and interpolated keyframes using Bézier curves. Their method is theoretically sound and excellent, but it cannot be applied directly to a polygonal mesh because it depends on the exponential and logarithmic maps between a manifold and its tangent space at identity. For practical usage, Martínez et al. [8] introduced a geodesic Bézier curve defined on a triangular mesh. They extended the de Casteljau algorithm in a triangular mesh by employing geodesic linear interpolation and generated a Bézier curve lying exactly on the domain mesh. In this study, the geodesic Bézier curve is further extend to geodesic Hermite spline curves to interpolate given points on a triangular mesh smoothly.

Geodesic Operations on Mesh
Computing a geodesic Hermite spline curve on a triangular mesh requires extended operations between the points and vectors on the mesh. In this section, these operations which include (i) the definition of a geodesic vector, (ii) subtraction between a pair of points, and (iii) parallel translation of a geodesic vector over the mesh, are defined.

Geodesic Vector on Mesh
Given a triangular mesh M, let p be a point on a triangle f ∈ M and v be a vector in the tangent space T p M. A "geodesic vector", denoted by v p on M is defined as the straightest geodesic proposed by Polthier et al. [10]. A straightest geodesic is a directed polyline on M, which emanates from p in direction v. Computing the straightest geodesic, also known as "geodesic tracing", can be summarized as follows.

1.
Compute the intersection point p of line p + tv with the edge of the triangle f containing p.

2.
Determine the next straight direction v by unfolding the adjacent triangle that shares the intersecting edge into a common plane.

3.
Repeat steps (1) and (2) by setting p = p and v = v until the length of the geodesic is equal to v . Figure 2 shows the geodesic vector v p (in pink) generated by the geodesic tracing algorithm. Please note that geodesic tracing can be considered the discrete exponential map exp p : T p M → M defined by exp p (v) = v p and it is restricted in the local neighborhood of p to construct one-to-one correspondence (see the rightmost image in Figure 3).  The angle sum Θ i of a vertex V i is defined as: where θ i is the interior angle at V i of triangle f j sharing V i . Intuitively, Θ i captures the flatness of the vertex, and (2π − Θ i ) is often used as a discrete analog of the Gaussian curvature. During geodesic tracing, the ambiguity of the unfolding triangle occurs when the intersection point p coincides with the mesh vertex V i . In this case, a new direction v is chosen such that the geodesic curvature is minimized. More specifically, the direction that bisects the angle sum Θ i is chosen as the new direction v while preserving the symmetry property of the straightest geodesic. In Figure 2, the geodesic vector meets vertex V i ; thus, the direction with symmetric angles(θ l = θ r ) is chosen for the next direction. Figure 3 shows 30 straightest geodesics emanating from a single point on a torus model. As the figure goes to the right, it shows a longer and longer straightest geodesic with self-intersections possibly.
Scalar multiplication of a geodesic vector can be achieved by scaling its length, and the addition of two geodesic vectors, denoted by v p ⊕ u q , can be computed by successive geodesic tracing if the endpoint of v p coincides with q. Otherwise, it is necessary to move u q to u p over the mesh so that its start point becomes p, as detailed in the following subsection. Here, the symbols ⊕ and are used for addition and subtraction to distinguish them from those in Euclidean space. Please note that scalar multiplication combined with addition allows a linear combination of geodesic vectors, which produces a new geodesic vector on M.

Geodesic Subtraction of Points
Let p and q be points on a triangular mesh M. In Euclidean space, subtraction of p from q yields a vector that does not lie on the mesh in general. However, based on the definition of a geodesic vector on M, subtraction of p from q should result in a geodesic vector from p to q. In this section, the subtraction of points on a triangular mesh, equivalent to, the addition of a geodesic vector to a point, is defined.
For this, we employ the shortest geodesic from p to q, denoted by pq. As mentioned in Section 2, there exist many efficient algorithms for geodesic computation on a polygonal mesh. Here, the shortest geodesic is computed by propagating the geodesic distance field and then back tracing using the distance field [13,15]. Figure 4 shows the geodesic distance field propagated from a source point p and the shortest geodesic (in pink) to the target point q. Alternatively, other methods [11,12,[18][19][20][21] can be replaced to improve the computational efficiency and accuracy. In general, the shortest geodesic consists of a single straightest geodesic or a sequence of straightest geodesics connected at the saddle vertex [11,13] which is a vertex V i with Θ i > 2π. If the shortest geodesic pq consists of a single straightest geodesic, or equivalently passes through only edges, it can be represented by a single geodesic vector v p as: where d g (p, q) is the length of pq, andṽ p ∈ T p M is the normalized tangent vector from p to the intersection of pq with the edge of the triangle containing p-see Figure 5a. If the shortest geodesic pq passes through the saddle vertices V k (k = 1, · · · , n), it cannot be represented by a single geodesic vector, but the sum of the geodesic vectors as follows: Figure 5b shows an example of the shortest geodesic, consisting of two geodesic vectors connected at the saddle vertex V 0 . In such case, it is possible to retain all geodesic vectors to represent the shortest geodesic exactly. However, such case rarely occurs and thus approximation to pq is taken by selecting the first geodesic vector v p for a practical usage. Finally, the subtraction of p from q is defined as the geodesic vector v p and the following notation is used to distinguish it from one in Euclidean space:

Parallel Translation of Vector
Let q be a point and v p be a geodesic vector on a triangular mesh M. Because the start position of v p does not coincide with q, one cannot add v p to q directly. To make this operation allowable, it is necessary to move v p to v q so that it starts from q. In this section, the geodesic operations of points and vectors in M are completed by introducing the parallel translation of a geodesic vector over M.
To move v p to q in parallel, the shortest geodesic pq is used again. Let θ be the angle between v p and pq (see Figure 6). One extends pq by geodesic tracing and determines a unit vectorṽ ∈ T q M such that the extended pq andṽ make an angle θ. Finally, the parallel translation of v p to q, denoted by v q , can be computed as follows: Intuitively, parallel translation of a geodesic vector on a triangular mesh is achieved by moving the geodesic vector while maintaining the angle formed by the shortest geodesic. Figure 6 illustrates the translation of a geodesic vector over a sphere model. The geodesic operations for points and vectors make it possible to perform algebraic computations on a triangular mesh, such as a linear combination of vectors and a convex combination of points. For example, Figure 7 shows the centroid s = p+q+r 3 of triangle ∆pqr on a sphere, which can be computed by the geodesic operations as follows:

Geodesic Hermite Spline Curve on Mesh
Now, an interpolation curve on a triangular mesh can be created using the geodesic operations introduced in the previous section. First, the interpolation parameters are determined by the geodesic length between consecutive input points. Second, the tangent vector for each point is computed by subtracting of consecutive points and third, the intermediate control points for each Bézier curve are computed by appropriate geodesic operations. Finally, the geodesic Hermite spline curve consists of a set of cubic Bézier curves generated by this method.

Hermite Spline Curve
The method is based on Hermite interpolation. Therefore, first, the conventional Hermite spline curve is summarized. Given two points p 0 and p 1 , and their tangent vectors v 0 and v 1 , the Hermite interpolation curve C(u) defined on u ∈ [a, b] while satisfying C(a) = p 0 , C(b) = p 1 and C (a) = v 0 , C (b) = v 1 can explicitly be constructed as: where H 3 i (t), 0 ≤ t = u−a b−a ≤ 1 is the i th Hermite basis function defined as follows: Based on the decomposition of Hermite basis functions H 3 i (t) into Bernstein ones B 3 i (t), the interpolation curve C(u) can be represented by a cubic Bézier curve as: where 0 ≤ t = u−a b−a ≤ 1 and the control points are determined as follows: In this paper, a Bézier form (3) for a Hermite curve is employed because its intermediate control points b 1 and b 2 can be handled efficiently by the geodesic operations defined in the previous section. By applying the Hermite interpolation method to a pair of points p i and p i+1 , and their tangent vectors v i and v i+1 , one can construct a Hermite spline curve with C(u i ) = p i and C (u i ) = v i , where u i is an interpolation parameter for p i .

Interpolation Parameters
In general, interpolation parameters u i and tangent vectors v i are not given; thus, it is necessary to determine them from given points p 0 , p 1 , · · · , p n . Common choices for interpolation parameters of a Hermite spline curve include uniform, chord length (chordal), and centripetal methods [25]. The chordal or centripetal method can solve the problem of loops or self-intersections caused by uniform parameters [26]. Yuksel et al. [27] showed that the centripetal method has superior properties compared with other methods. Here, chordal and centripetal methods are both used, and their effects are compared in Section V.
Given the points p 0 , p 1 , · · · , p n , the interpolation parameters can be determined as follows: The chord length and centripetal parameters are obtained when α = 1 and α = 1 2 , respectively. To extend this method to the points on a triangular mesh, the Euclidean distance · should be replaced with the geodesic distance d g (·, ·) as follows: (4) Figure 8 shows an example of the chordal parameters for four points on a sphere model using the geodesic distances.

Geodesic Tangent Vectors
To construct a Hermite spline curve with C 1 continuity at C(u i ) in Euclidean space, it is necessary to assign a tangent vector v i to p i . Various methods for determining the tangent vectors of the Hermite spline curve exist, and the cardinal spline is a widely used technique. Let u i be the interpolation parameter for p i . The tangent vector v i is computed as follows: where c ∈ [0, 1] is a tension parameter that controls the magnitude of the tangent vector. However, the tangent vector in (5) does not lie on a triangular mesh and therefore cannot be used for constructing an interpolation curve on a mesh. In this paper, the Catmull-Rom spline is employed with c = 0.5 [28] and the subtraction between two points is replaced with the geodesic operation to produce the "geodesic tangent vector" as follows: where the geodesic vector p i+1 p i−1 moves to p i by parallel translation (see Figure 9). For an open Hermite spline curve, the geodesic tangent vectors v p 0 and v p n for the endpoints can be excluded or zero vectors. Here, zero vectors are assigned to the endpoints. A closed Hermite spline curve is quite useful for selecting a region bounded by the curve and, thus, is effectively used for such applications as mesh-cutting, texture mapping, and mesh hole filling. To construct a closed Hermite spline curve such that C(u 0 ) = C(u n ), two endpoints should be identical, p n = p 0 , and the common geodesic tangent vector should be assigned to both points as follows: . Figure 9 shows how to determine the geodesic tangent vector for each point p i for constructing the closed Hermite spline curve.

Geodesic Hermite Spline Curve
Once the interpolation parameters u 0 , u 1 , · · · , u n and geodesic tangent vectors v p 0 , v p 1 , · · · , v p n have been obtained, one can construct a Hermite spline curve that interpolates given mesh points p 0 , p 1 , · · · , p n . The Hermite spline curve consists of piecewise C 1 -continuous Bézier curves. Here, the focus is on a single geodesic Bézier curve C i (t) satisfying C i (u i ) = p i and C i (u i ) = v p i . As mentioned, C i (t) can be represented by a cubic Bézier curve C i (t) = ∑ 3 k=0 b k B 3 k (t), and the control points b i can be computed as follows: Because a cubic Hermite curve is used, both end control points b 0 and b 3 are simply set to be p i and p i+1 , respectively, and intermediate control points b 1 and b 2 are calculated by geodesic tangent vectors and constants related to the difference between interpolation parameters u i and u i+1 , as in (6). Equation (6) is equivalent to (3), except that the operations are replaced with geodesic ones-that is, the geodesic operations on a triangular mesh make it possible to reuse the conventional Hermite spline technique for constructing an interpolation curve on the mesh.

Experimental Results
The proposed method was implemented on a PC with an Intel i7-9700F CPU, 24.0-GB main memory, and an NVIDIA GeForce RTX 2060 graphics card, using the C++ language. The geodesic Hermite spline curve discussed can be used in various applications, including curve design on a polygonal mesh, as well as mesh-cutting and simulation of a moving object along the interpolation curves. The effectiveness of the method was demonstrated by designing several sophisticated curves on triangular meshes and applying them to various applications such as mesh-cutting, segmentation and simulation.

Hermite Spline Curves on Mesh
To create a geodesic Hermite spline curve, a user must select at least two points on a triangular mesh at runtime. Figure 11 shows four user-selected points (in green) on a sphere model and the shortest geodesics (in pink) between a pair of points. Although a sequence of shortest geodesics can be used for interpolating given points, they are inadequate for providing C 1 -continuous curves. The Hermite spline curve interpolating four points is shown in blue, where three Bézier curves are joined with C 1 continuity. In this example, chordal parameters and Catmull-Rom splines were employed to determine the geodesic tangent vectors. The intermediate control points determined by the geodesic operations are shown in black. Given n mesh points, this method generates (n − 1) cubic Bézier curves with C 1 continuity. Please note that the generated curve shows a symmetry property, i.e., the same interpolation curve is obtained if the mesh points are specified in the reverse order.   Table 1, and the time for constructing geodesic Hermite spline curves for each model is listed in the last row in Table 1, where sampling points on the interpolation curve is not included in the construction time.  The proposed technique shows a real-time performance for most of test examples presented in this paper. If the user moves an interpolation point to a new position, the method generates a new Hermite spline curve interpolating the new position in real time (see the Video S1).

Mesh-Cutting with Closed Curve
The selection of a region of interest, followed by cutting is one of the fundamental operations in geometric modeling and processing. In particular, a sophisticated boundary curve is required for mesh segmentation and partial texture mapping. The method provides an effective solution to these applications. By creating a closed Hermite spline curve, the user can select a region of interest on the model and cut the region with ease. Figure 13 shows the process of generating a closed Hermite spline curve near the eye of the Armadillo model. Four interpolation points near the left eye region were selected interactively and the resulting interpolation curve was constructed to separate the eye of the Armadillo model. Figure 14a shows the selection of the left eye of the model using the closed interpolation curve and Figure 14b shows the result of cutting the selected region.

Simulation
The Hermite spline curve on the mesh can also be used for the simulation of an object or camera moving along the generated path. Figure 15a shows the Bunny model moving along a curved path on a spherical model. In the figure, the yellow curve is the Hermite spline curve on the sphere, and the red curve indicates the path through which the Bunny model passed. Figure 15b shows the local axis of the curved paths on a sphere, Bunny and cylinder models, all of which are generated by the method presented. (a)

Control Parameters
Our technique provides the user with two parameters α in (4) and c in (5) for controlling the shape of the generated Hermite spline curve. In this section, the effects of these control parameters when creating a geodesic Hermite spline curve are discussed. Figure 16 shows the effects of changing the tension parameter c ∈ [0, 1] in (5). The tension parameter controls the length of the geodesic tangent vector, and therefore influences the sharpness of the generated curve at the interpolation points. In this example, 12 points are given for the open Hermite spline curve (in yellow) and 9 points are given for the closed Hermite spline curve (in pink). The interpolation parameters are computed by the chord length method for all tests. As the tension parameter increases, the generated Hermite spline curve approaches the shortest geodesic (see the Video S1).  Figure 17 shows the effects of changing the control parameter α in (4) which influences the interpolation parameters u i , where pink curves are generated by centripetal parameter (α = 0.5) and blue ones are generated by chordal parameter (α = 1). Depending on the interpolation parameters of the given points, the method produces slightly different results. As the tension parameter increases, the two curves become identical. These examples demonstrate that our method follows the same control mechanism as the conventional Hermite spline curves in Euclidean space.

Conclusions
An effective method for constructing a geodesic Hermite spline curve on a triangular mesh while interpolating the given mesh points was presented. To extend the conventional Hermite interpolation method to a triangular mesh, new geodesic operations between points and vectors on a triangular mesh were defined by using the concept of the shortest geodesic and straightest geodesic. The geodesic Hermite spline curve was constructed using the conventional Hermite interpolation method, except that the algebraic operations were replaced with the new geodesic operations.
The method is simple to implement and easy to use. A user simply selects the points on a triangular mesh, and then the interpolation curve is generated in real time. Moreover, the generated curve can be edited by simply moving the interpolation point to a new position interactively. In the experimental results, the effectiveness of the technique was demonstrated by showing several examples and applying the method to various applications.
In the current implementation, the subtraction of two points was approximated with a single geodesic vector, which can be improved by keeping all geodesic vectors. In future work, it is planned to apply the geodesic operations to other applications, such as computations of Voronoi diagram and medial axis bounded by geodesic curves on a polygonal model.