Abstract
We present a novel framework to compute geodesics on implicit surfaces and point clouds. Our framework consists of three parts, particle based approximate geodesics on implicit surfaces, Cartesian grid based approximate geodesics on point clouds, and geodesic correction. The first two parts can effectively generate approximate geodesics on implicit surfaces and point clouds, respectively. By introducing the geodesic curvature flow, the third part produces smooth and accurate geodesic solutions. Differing from most of the existing methods, our algorithms can converge to a given tolerance. The presented computational framework is suitable for arbitrary implicit hypersurfaces or point clouds with high genus or high curvature.
1. Introduction
Geodesics play an important role in differential geometry. The applications range from finite element computation to computer aided geometric design, from computer animation to robotic navigation, and from brain flattening and warping in computational neuroscience to machine learning on manifolds. However, little research has been undertaken on point clouds so far, despite the fact that affordable 3D acquisition devices are becoming increasingly available. One of the main challenges is that all computation must be carried out on a set of scattered points rather than parameterized surfaces or meshes.
In this paper, we will tackle this challenge and concentrate on two scenarios: implicit surfaces and point clouds, which are the common forms used for representing point set surfaces. The presented geodesic computation framework works well for implicit surfaces and point clouds and can be easily applied to high dimensional datasets. Our framework consists of three parts, approximate geodesics on implicit surfaces, approximate geodesics on point clouds, and geodesic correction. Although point set surfaces might be given in an implicit form in many applications [1–6], it is impractical to generate an implicit surface for every point cloud.
Our work is motivated by the results presented in [7], which computes intrinsic distance functions and geodesics on implicit hypersurfaces by embedding surfaces into a Cartesian grid. The computation is performed with the well-known fast marching method [5]. This approach can handle surfaces of a high curvature, such as that shown in Figure 1. However how to generate a Cartesian grid for an implicit function remains a challenging issue, since the size of the Cartesian grids grows exponentially when they are subdivided iteratively, rendering the computing process unbearably slow. Moreover, the resulting geodesic paths lie on a grid offsetting from the surface; that is, they do not lie on the surface. To the best of our knowledge, other algorithms, such as variational curve design schemes [8–10], are unable to handle surfaces like that shown in Figure 1.
Our research is built on the work of [4], which presents a particle based approach of modelling implicit surfaces, and [11], which presents rendering point set surfaces by moving least squares. Our framework differs from the previous methods in that the geodesic computation is directly performed on the implicit surfaces or point clouds, not on the proximity mesh or grid. Our main contributions include the following.
(1) Accurate Geodesic Computation. When an approximate geodesic path is available, which may be an approximate path lying on the grid offsetting from the point set surface or a zigzag path on the implicit surface, our method can effectively reach a smooth and accurate geodesic solution on the point set surface. The distinct advantage is that the precision of the geodesic solution is controllable.
(2) Arbitrary Implicit Surfaces or Point Clouds. Our algorithms can work well on arbitrary implicit surfaces (or point clouds) irrespective of genus or curvature (positive, nonpositive, high, or low).
Moreover, our presented algorithms can perform geodesic computation on both 3D and high dimensional manifolds.
The remainder of this paper is organized as follows. In Section 2, we briefly review the related work. Our computational framework is presented in Section 3. The details of the implementation and discussions are given in Section 4. Finally, our conclusions are given in Section 5.
2. Previous Work
For a triangle mesh, the pioneering work is the MMP algorithm proposed by Mitchell et al. in [12], which provides a solution for the “single source, all destination” shortest path problem. For the worst case, the running time is . Another geodesic algorithm was presented in [13], whose time complexity is . Almost all recent work focuses on producing approximate geodesics with a guaranteed error bound. For example, additional edges are introduced into a mesh for geodesic computation [14]. Speeding up the MMP algorithm with a merging operation is also of interest [15]. A well-known algorithm is the fast marching method (FMM) [5], which computes approximate geodesics in time. Working on a triangular mesh, the FMM’s distance functions are sometimes error-prone. The resulting geodesic path is not always accurate. Some correction methods are consequently proposed as a post procedure following the FMM, such as applying the “straightest geodesics strategy” [16] to correct the geodesic path in [17].
In contrast, there have been only a few methods relating to implicit surfaces. A parametric curve can be easily represented in an implicit form, but an implicit curve cannot always be represented parametrically in dimensions other than two or one. A practical treatment in dealing with an implicit surface is to embed it into a mesh. In [7], an implicit surface is embedded into a Cartesian grid for geodesic computation. However, the resulting geodesic paths lie on the Cartesian grid, not on the implicit surface. Witkin and Heckbert in [4] proposed an alternative method of directly sampling and controlling implicit surfaces by using a particle system. Hence, one of the objectives of this paper is to develop an approach to computing geodesics directly on an implicit surface to improve the computational accuracy.
For point clouds, several authors have presented their individual methods. Klein and Zachmann [18] approximated geodesics as shortest paths on a geometric proximity graph over point clouds, called the spheres-of-influence graphs (SIG). Ruggeri et al. in [19] further improved the accuracy of the approximate geodesics defined on the SIG. Hofer and Pottmann [8] proposed to compute geodesic curves as energy minimizing discrete curves constrained on a moving least square surface. However, a distinct defect is that these methods are usually sensitive to noise. Hence, it is hard to control the precision. Our proposed algorithms can converge to a desired solution with some specified tolerance.
3. Geodesic Computational Framework
A geodesic is calculated between two specified points on a surface. The proposed approach in this paper calculates the geodesics of all pairs of specific sample points on an implicit surface or a point cloud. Usually, this set of specific sample points is given in advance. However, for measuring an implicit surface or a point cloud, resampling the surface is a vital step. We discuss this issue on two scenarios as follows.
3.1. On Implicit Surfaces
The basic idea of our sampling method is to drive the particles to float from the sources (where the geodesics start) to the destinations (where the geodesics end, also called sinks) on an implicit surface. The resulting trajectories are used as the initial estimates of the geodesics between the sources and the destinations. In order to produce smooth geodesics on the surface, in this paper, we introduce two modifications to the - formulation proposed in [4]. (1) We modify the particle density model and add a curvature term. This makes particle distribution dependent on the surface curvature resulting in more optimal distribution of the particles. (2) We present a repulsion-attraction scheme. It can not only yield smooth particle trajectories but also further speed up the process. To compute geodesics, the resulting particles are to form a connection graph covering the surfaces. However, unlike a mesh, the edges of the connection graph are the trajectories of the particles, not straight lines. The nodes of the connection graph (i.e., the initial given samples) can be arbitrarily placed on the surface. To obtain a smooth and accurate solution, geodesic curvature flow will be employed as the geodesic correction procedure, which is addressed in Section 3.3.
Our algorithm consists of three steps: (1) particle based sampling on implicit surfaces, (2) Dijkstra’s algorithm, and (3) geodesic correction. We first illustrate how to manipulate particles directly on an implicit surface, followed by geodesic computation.
3.1.1. Density Model
We are concerned with only the particle position on stationary surfaces here. For simplicity, the surface deformation terms are ignored from the - formulation thereafter. Consider an initial collection of floating particles randomly lying on an implicit surface of the implicit function , , where is the trajectory of the th particle. We use boldface letters to denote vectors and matrices and italics for scalars thereafter. The particles then repel each other so as to reach equilibrium with a uniform distribution on . The particle motion can be described as follows: where is the derivative of th particle w.r.t. time , representing a tangent vector at which is orthogonal to the normal , is a velocity (or estimate) of , and is a constant. In our experiments, is set to 1. The first term of the right side of (1) is viewed as the tangent component of that pushes the particles away in the tangent plane at . This results in the particles offset to , especially in the areas of high curvature. Hence, the second term of (1) plays a role that pulls back the offset particles towards by using the Newton-Raphson approximation scheme. (For details, refer to [4].)
However, the original scheme of (1) has many deficiencies. The main problem is how to estimate the velocity in (1). Usually, the original - method can produce a homogeneous distribution of particles on an implicit surface. But, it is difficult to give an inhomogeneous distribution based on curvature. It has been observed that the new particles are only pushed out onto flat areas, which will become too crowded, while the particles are fleeing from the high curvature areas.
3.1.2. Curvature Term
A Hessian matrix is usually used to describe the differential structure of an implicit surface , and the eigenstructure of the Hessian matrix is used to gauge the curvature of . For an implicit surface in , the Hessian matrix at can be given as . It is well-known that the maximum and minimum eigenvalues of the Hessian at , respectively, correspond to the principal curvatures of at , while the corresponding eigenvectors are, respectively, the principal directions in . For an implicit surface in , the principal curvatures may be defined in a directly analogous fashion. Hence, the curvature vector of at a point can be described by the eigenvectors of the Hessian at as follows: where is the th eigenvalue of and is the corresponding eigenvector. The curvature term of (2) only depends on the differential structure of rather than others. It will be introduced into the estimate of in (1) as an extra force.
3.1.3. Repulsion-Attraction Scheme
For point to point geodesic computation, the particles from a single point (source) are expected to be attracted into a specified destination point (sink) quickly; therefore the trajectory of the particle which reaches the sink the quickest is the closest approximation of the shortest path from the source to the sink. Our repulsion scheme is thus extended to a repulsion-attraction scheme as follows.
The energy functional of attraction for each particle is defined as where denotes a sink on , is a constant, and denotes the sink number. The attraction force that each particle receives from is independent of each other. Minimizing by gradient descent with respect to the particle position yields the desired attraction force exerted to , that is, . Hence, there are three forces exerted to here, that is, the repulsion, the attraction forces, and curvature term. The resultant force is written as where denotes the number of nonsink particles. In the repulsion-attraction scheme, the sink number is a variable. The determination of the sinks is given later. It can be noted that there exist two balance constants in order to estimate . No parameters need to be tuned during iteration. This is desirable for the stability and effectiveness of the particle system. In our experiments, the particle system works well when setting both and to 1.
3.1.4. Approximation of Geodesics
The scheme of (1)–(4) describes the particle behaviour on an implicit surface. For the purpose of geodesic computation, the repulsion-attraction scheme is employed here. Suppose that the sources on surface are known. The goal is to compute the geodesics from some specified source to all destinations. Herein, we provide a procedure for particle generation. The particles are then distributed over by the scheme of (1)–(4). The basic idea is to utilize the rule of particle repulsion-attraction; that is, particles from the same sources repel each other; otherwise, they attract each other.
Figure 2 shows the particles’ trajectories form a network that covers the implicit surface. The given sources (i.e., A, B, C, D, and E) generate the new particles according to the distances from each source to their individual nearest neighboring particles. Let the derived particles from the same source be the sinks of the particles coming from other sources. The particles from the same sources repel each other; otherwise, they attract each other. Note that the sources are fixed, while the sinks are floating on here. When two particles derived from two different sources meet, these two sources are connected by the particles’ trajectories (see the dashed lines in Figure 2). At that time, the particles derived from these two sources are regarded as the homologous particles and repel each other. As a result, all the sources will be connected finally by the particles’ trajectories as shown in Figure 2. Our presented sampling approach is summarized as follows.
Procedure of Sampling an Implicit Surface(i)Input a given threshold , implicit function , and a set of sources.(ii)A loop is as follows:(a)each source generates new particle, if the shortest distance to its neighbours is more than ;(b)repulsion-attraction scheme equations (1)–(4);(c)if , where particles and are derived, respectively, from two different sources, then these two sources are connected by the recorded trajectories of and ;(d)recording the current location of each particle.(iii) This is done, until all the sources have been connected.
Moreover, we can apply Dijkstra’s algorithm [20] to the resulting connected graph for the approximate geodesic paths on . The whole algorithm for geodesic approximation is then described as follows.
Algorithm 1 (for approximate geodesic paths on implicit surface). Consider the following.(1)Input a set of sources.(2)Apply the procedure of sampling an implicit surface.(3)Apply Dijkstra’s algorithm to the resulting connected graph.
Now, let us consider the time complexity of the algorithm for sampling the implicit surface. Assume that there are given sources and drifting particles are produced on each time. Applying the repulsion-attraction scheme of (1)–(4), we can note that it involves an inner iteration, that is, to pull back new particles to by the Newton-Raphson iteration. In our experiments, this can be accomplished through a small number of iterations. For simplicity, we ignore this inner iteration here. There are particles on at the th iteration. It takes times for the repulsion-attraction computation between all pairs of points at that time. Hence, the upper bound of computational time is estimated as , where . When is large, that is, becomes small, only a small number of iterations are required to generate the connection graph.
Applying Dijkstra’s algorithm to the resulting graph, we can compute the approximate geodesic paths of all pairs of sources at the cost of . Hence, for all pairs of sources, the complexity of Algorithm 1 is of in total. This means that adding the process of implicit surface sampling does not cause the overall running time to increase exponentially.
In general, the implicit surfaces are utilized to fit discrete sampling points. It is straightforward to select the sources from the original sampling points manually. Moreover, the resulting geodesic paths are only approximations on the implicit surface. We will show that the approximate paths can converge to geodesic curves on the implicit surfaces through a geodesic correction procedure described in Section 3.3.
3.2. On Point Clouds
For a point cloud, an alternative is to approximate the point set surface by the Cartesian grid as used in [7]. The basic idea is conceptually simple, that is, to determine the Cartesian grid envelope which surrounds the implicit hypersurface and then apply “continuous Dijkstra”-like strategy to the grid for generating the geodesic estimations.
Herein, we set the mean of the points within a specified cell of the grid as the source node on the surface. This makes the source nodes evenly distributed over the point cloud surface. Moreover, the unorganized points are in general much more than the source nodes. However, the drawback is to build up a whole grid with request for an additional space, where denotes the dimensionality of manifolds and denotes the node number. It can be noted that this is a sparse multidimensional array. To avoid such large space complexity, the simplest approach is to store the nonempty cells into an array. Thanks to the ANN searching [21], we may then conveniently determine the neighboring cells of a query cell by it. The ANN searching usually needs spaces for the nonempty cell storage and an extra time of partitioning the point clouds into a data structure (e.g., - tree), that is, .
Based on the resulting grid, we can propagate the distance information from a specified node to all other nodes on it in a “continuous Dijkstra”-like manner. This is illustrated in Figure 3. The algorithm for geodesic approximation on point clouds is summarized as follows. As a result, computing all pairs of geodesics is of complexity.
Algorithm 2 (for approximate geodesics on a point cloud). Consider the following:(i)input: a point cloud , a source within the cell , and the intervals of a grid, that is, , , ;(ii)build up a searching tree for ANN running on the grid;(iii)loop:(a)determine the neighboring cells of the visited cell set by ANN searching on the grid, and denote them as ;(b)loop (within ):set the mean of points within the th neighboring cell as the destination node, and denote it as , which is not visited;obtain the “shortest” edge linking the and the visited cell set, that is, the “shortest” edge to the visited destination node set;combine this edge with the existing geodesic path to form the initial estimation of the geodesic from to ;(iv)until all the nonempty cells are visited.
3.3. Geodesic Correction
A geodesic curve has vanishing geodesic curvature everywhere. This suggests that an exact geodesic curve should have zero geodesic curvature at all of its points. For a discrete representation, we have to add another requirement that all samples must lie on the implicit surface. The presented correction approach below can guarantee that the sampling points of a curve converge to the implicit surface or point cloud surface and the curve on the surface converges to the geodesic. It is therefore more accurate than the existing mesh or grid based algorithms.
In order to make a curve converge to a geodesic curve, geodesic curvature flow is adopted here. This is because geodesic curvature flow eliminates the geodesic curvatures pointwise on a curve. For clarity, we first consider a regular surface and two endpoints , , between which we intend to compute a geodesic curve on . Let be an initial smooth parametric curve on with and , where denotes the arc-length parameter of the curve. If these two endpoints are fixed, we deform by the geodesic curvature flow on . The curve converges to a geodesic curve. The geodesic curvature flow is described as where denotes the geodesic curvature vector of . The geodesic curvature vector can be given as where is the 2nd-order derivative of w.r.t arc-length and is the normal to . This flow is also known as the Euclidean curve shortening flow [22] and has been widely applied in computational physics, image processing, and material sciences.
For a discrete case, the geodesic curve has to be approximated by a polyline containing a set of nodes. Correcting the approximate geodesic curve is thus carried out on the resulting polyline. Since we have no parametric forms, the second-order derivative of w.r.t is approximated using a vector triangle as shown in Figure 4. This is the increment of the tangent vector at ; that is, . If ’s projection onto the tangent plane at vanishes, the geodesic curvature of (6) is also bound to vanish at . When the geodesic curvature vanishes pointwise on , the deforming curve given by (5) reaches its stable state, that is, the geodesic curve.
Moreover, we can write the discrete version of (5) in a matrix form as follows:where denotes the number of sample points on , contains the coordinates of sample points, and is a coefficient matrix. Then, the geodesic curvature of (6) can be rewritten as where denotes the normal vector at the th sample point . As a result, the geodesic curvature flow of (5) can be rewritten as To speed up the iteration of (9), we apply the backward Euler method [23] to it, which yields where denotes the iterative step length, which is an arbitrary constant due to the backward Euler method. To solve the scheme equation (10) efficiently, we may use the preconditioned conjugate gradient (PCG) method here. Let . It can be noted that is sparse and symmetric. By having a good precondition, it can significantly improve the convergence rate of the conjugate gradient solver. The simpler the precondition, the faster the convergence. The is a block diagonal matrix, while the is a block tridiagonal matrix. All the eigenvalues of lie in the interval of . Matrix A is well conditioned for typical values of . In our implementation, we employ the usual diagonal preconditioner ; that is, , which speeds up the iteration with almost no overhead.
In this paper, we need to deal with two cases, implicit surfaces and point clouds. For the former, it can be noted that ; that is, we can regard as a velocity estimate at in (1). This is the constraint of the particles moving on the geodesic curve. Consequently, we can set and substitute it into (1) leading to the following evolution equation: It can be noted that the 2nd term pulls back the sample points onto the implicit surface . All the sampling points will lie on and converge to the geodesic accordingly. In our implementation, we split (11) into a two-step scheme. The 1st step employs (10) for updating the in (11), while the 2nd step is the Newton-Raphson iteration given by For point clouds, there is no implicit function available and we have to seek an alternative approach for estimating the normal vectors. We also have to consider how to restrict the geodesics on the point cloud surfaces similarly to the 2nd term in (11). Because of the above issues, we employ the moving least squares [11] to the point clouds and estimate the normal vectors for the curvature flow of (10).
For a sample point , we solve for the best tangent plane at that minimizes where denotes all the sample points. Herein, we employ the weight that is defined as where is a constant reflecting the anticipated spacing between neighboring points. The unknowns include the ’s coordinates and the parameters of affine here. For simplicity, let and , where denotes an initial clustering centre and denotes the normal vector. Such the minimization problem of (13) can be solved in an iterative manner; that is, we start with and first approximate as follows: and then, fixing the normal , we solve the following nonlinear minimization problem with the single unknown : In our implementation, we found that there existed a satisfied solution within . Note that we select the nodes of the approximate geodesic paths from Algorithm 2 in Section 3.2 as the initial . Once and are available, it is natural to update the nodes accordingly, resulting in the resulting geodesic paths being restricted on the point cloud surfaces.
There are scenarios where a specified approximate geodesic path is given. In this case, we first compute the normal vectors and update the location of each node equation (13) and then apply (10) to it for improving the geodesic accuracy.
In order to sample a geodesic curve evenly according to the geometric features, the samples are distributed depending on the curvature of the geodesic curve. This can be given by the proposed repulsion-attraction scheme equation (4), which is modified as Note that all the samples on a curve repel each other and the curvature term is replaced by the 2nd derivative of . Substituting it into (1), the particle motion equation is then rewritten as It can be noted that (18) only contains the repulsion motion. Once the equilibrium is reached, the samples will be distributed depending on the curvature of the curve.
Let us summarize the geodesic correction procedure below.
Procedure of Geodesic Correction(i)For the approximate geodesic paths from Algorithm 1,(a)apply Newton-Raphson iteration in (12) for pulling back the particles;(b)apply the scheme (10) for correcting geodesics;(c)apply the scheme equation (18) for a distribution depending on curvature;(d)check the interval of the successive samples and interpolate new samples accordingly;(e)repeat until the geodesics converge.(ii)For the approximate geodesic paths from Algorithm 2,(a)apply the scheme of (13) to the input geodesic paths for updating the nodes and the associated normal vectors;(b)apply the scheme of (10) for correcting geodesics;(c)apply the scheme equation (18) for a distribution depending on curvature;(d)check the interval of the successive samples and interpolate new samples accordingly;(e)repeat until the geodesics converge.
The time complexity of the geodesic correction procedure depends on the number of the sample points on a geodesic, that is, in (10). First, we need to pull back the samples to implicit surfaces one by one in both scenarios. Then, (10) is employed individually. Since matrix of (10) is sparse and block diagonal, the multiplication of matrix-vector only costs a linear time in the implementation of the preconditioned conjugate gradient method [23]. Assume that pulling back the samples to the surface costs , and solving (10) in time and solving (10) can be iterated times. Thus, one geodesic costs around . For one single source geodesic computation, the upper bound of the whole time complexity can be estimated as , where denotes the destination number. Due to the backward Euler scheme and the precondition technique, the convergence rate of (10) is drastically improved. In general, 2 iterations are sufficient to reach the solution with the relative error less than . Thus, the time complexity of one single source geodesic computation is of . For all pairs of geodesics, it costs . We also need to store each geodesic individually, which costs space in total.
3.4. Error Analysis
In this section, we will show that regardless of whether the approximate geodesic paths are obtained from Algorithm 1 (such as the zigzag curves on implicit surfaces) or from Algorithm 2 (such as the polylines lying on the Cartesian grid), with different intervals, the geodesic correction procedure can achieve the same precision automatically.
In order to estimate the approximation error, we assume that the number of nodes is fixed; that is, the correction procedure does not generate new nodes. Let us consider the local structure of a curve with an osculating circle as shown in Figure 5(a). The local approximate error is described as where denotes the arc-length from to and is the chord length . Applying the cosine law yields , where denotes the radius of curvature. Substituting into err gives . For a whole geodesic curve, the error is expressed as where denotes the number of samples on a curve and and depend on the local curvature of a curve. With regard to (10), we assume that is constant for a given geodesic curve. The error of the whole geodesic curve can thus be given by where is the estimate of the curve’s length. Because of the scheme equation (18), the samples are distributed depending on the curvature of the geodesic curve. It can guarantee that the Cst is constant. This means that our algorithms can converge to a desired solution with some specified tolerance.
(a)
(b)
Moreover, if the number of samples is variable, Figure 5(b) further shows that raising can improve numerical stability effectively.
4. Implementation
Our experiments consist of two parts, implicit surfaces and point clouds. All the codes are run on MatLab in a Pentium IV 3.2 GHz PC with 4 GB RAM. Implicit surfaces were generated from three point clouds (i.e., Stanford bunny, hand, and sculpture) by using the FastRBF toolbox (at http://www.farfieldtechnology.com/products/toolbox/).
4.1. Implicit Surfaces
To illustrate the curvature dependent distribution, we first performed the repulsion-attraction scheme equations (1)–(4) on the 3D model of the Stanford bunny. The reason why we used this model is because the ears have high curvatures (as shown in Figure 6, there are two images for front and back views), which can test the effectiveness of our method. The particles are generated from 39 fixed sources (which are marked as ). The sources disperse the particles over the implicit surface by continuously generating new particles. Due to the curvature term equation (2), it can be noted that many particles remain at the high curvature regions in Figure 6 allowing higher curvature regions to be more densely represented.
We further performed Algorithm 1 described in Section 3.1 on a synthetic 3D model of a roll that has two high curvature regions and a large surface area as shown in Figure 7. To visually demonstrate the trajectories of the particles, the particles are generated from only 5 fixed sources (which are marked as ). It can be noted that our algorithm does not produce the approximate geodesic paths for all pairs of sources on . This is because the particle system equations (1)–(4) terminate once all sources have been linked into a connected graph. Hence, applying Dijkstra’s algorithm to the resulting graph produces all pairs of geodesic approximations on . These approximate geodesic paths were further corrected by the geodesic correction procedure described in Section 3.3 as shown in Figure 7(b). For comparison, the approximate geodesics and the corrected ones are both depicted in Figure 7(c). One can see that some trajectories have some small zigzag lines around the sources. This indicates that the movement of the particles is very random in the neighborhood of the sources, since new particles are generated constantly around the sources. Moreover, the roll surface is a cylinder-like surface. Due to the random movement of particles, it is possible to generate two approximate geodesic paths sharing the same endpoints on the two sides of the surface. Performing the geodesic correction procedure on them results in two different geodesic curves on the surface (see yellow and black lines in Figure 7(d)). This can be easily understood by looking at the great circle on a sphere or a cylinder.
(a)
(b)
(c)
(d)
If the sources are evenly distributed over the surface, Algorithm 1 can generate the connected graph quickly. To illustrate it, we performed Algorithm 1 on a 3D human hand model, which contains 272 vertices. These vertices are regarded as the fixed sources and continue to generate new particles. However, the resulting network shown in Figure 8(a) is not a reasonable mesh for modelling the human hand. How to remesh it is beyond the remit of this paper. Herein, our goal is to carry out Dijkstra’s algorithm on a resulting connected graph for initial geodesic estimations. Each edge in Figure 8(a) corresponds to a particle trajectory, which is, namely, the initial geodesic path between two sources. Following Dijkstra’s algorithm, the geodesic correction procedure was further performed on the resulting initial geodesic estimations. We show the geodesics of “one source to all destinations” in Figure 8(b).
(a)
(b)
We also performed Algorithm 1 with the geodesic correction procedure on 2 implicit surfaces, that is, human hand and sculpture, for computing all pairs of geodesics. For graphical rendering purposes, only 3 sets of the geodesics of “one source to all destinations” are depicted on the surface of each model in Figure 9. Table 1 shows the running times. It can be noted that our method can handle the genus surfaces (e.g., sculpture surface with genus 2 in Figure 9).
(a)
(b)
4.2. Point Clouds
For a point cloud, we first utilize a Cartesian grid to cover it and then apply Algorithm 2 in Section 3.2 to the resulting grid. Independent of the cell size, the geodesic correction procedure in Section 3.3 can automatically interpolate new samples for the improvement of the geodesic quality. We performed Algorithm 2 on a semisphere that is covered by two Cartesian grids with different cell sizes as shown in Figure 10. It is evident that the corrected geodesic curves by the geodesic correction procedure are becoming the same if these two curves contain the same number of samples. It can also be noted that the error depends on the sample number rather than the original grid cell sizes.
(a) Cell size of 0.25
(b) Cell size of 1.5
For the comparison of the running times of Algorithms 1 and 2, we performed Algorithm 2 with the geodesic correction procedure on these 2 point clouds in Figure 9. The running times are shown in Table 1. It can be noted that the running time of Algorithm 2 is less than that of Algorithm 1. However, the subsequent application of the geodesic correction procedure takes nearly the same time. This is because, for the same specified precision, the procedure of geodesic correction automatically inserts new samples accordingly. This will result in the similar sample size regardless of Algorithms 1 or 2.
To further illustrate the effectiveness of Algorithm 2, we performed it on 3 point clouds with complex topological and geometric structures. The results are shown in Figure 11. It can be noted that Algorithm 2 can deal with point clouds with genus (e.g., point cloud of jar with genus 2).
5. Conclusions
In this paper, we have presented a novel framework for computing geodesics on implicit surfaces and point clouds. The main feature of this work is its ability to produce smooth and accurate geodesics on the surface. This is due to the introduction of geodesic curvature flow in our framework. In addition, since the geodesic correction procedure can distribute the sample points according to the curve curvature and further interpolate new samples automatically, the resulting geodesic solution can achieve the specified precision. Experiments have shown that our algorithm works well on arbitrary implicit surfaces and point clouds.
For a general case, if no sources are given, the implicit surface concerned must be modelled beforehand. The surface modelling step is not included here. Moreover, it can be noted that the geodesic correction procedure costs most of the running time compared to Algorithms 1 and 2. This is because geodesic curvature flow is introduced. Optimisation of the computational efficiency and accuracy of the geodesic correction procedure is left as future work.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.