- Email: [email protected]
A geometric method for computation of geodesic on parametric surfaces
Accepted Manuscript A geometric method for computation of geodesic on parametric surfaces
Peng Zhang, Ronglei Sun, Tao Huang
PII: DOI: Reference:
S0167-8396(15)00096-5 http://dx.doi.org/10.1016/j.cagd.2015.08.001 COMAID 1502
To appear in:
Computer Aided Geometric Design
Received date: Revised date: Accepted date:
17 November 2014 8 August 2015 8 August 2015
Please cite this article in press as: Zhang, P., et al. A geometric method for computation of geodesic on parametric surfaces. Comput. Aided Geom. Des. (2015), http://dx.doi.org/10.1016/j.cagd.2015.08.001
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Highlights • • • •
An efficient geometric algorithm aims to trace geodesic on parametric surfaces. Independent of the complex description of the geodesic equations. Simpler and faster than the existing geometric method. Step size of the algorithm adapts to the geometry of parametric surface.
A geometric method for computation of geodesic on parametric surfaces Peng Zhang, Ronglei Sun*, Tao Huang School of Mechanical Science and Engineering, State Key Laboratory of Digital Manufacturing Equipment and Technology, Huazhong University of Science and Technology, Wuhan, 430074, China *Corresponding author. Tel.: +86 27 87559416. E-mail address: [email protected]
A sketch of used notations Notation
Description
Used space
S C
Parametric surface Geodesic curve on S Preimage of C in parametric domain Ideal geodesic curve on S Starting point on curve C Unit tangent vector of C at Q0 Preimage of Q0 in parametric domain Direction of C ' at q0 Point lying in the neighborhood of q0 on C ' Three-dimensional image point of q1 on S
3D space 3D space Parameter space 3D space 3D space 3D space Parameter space Parameter space Parameter space 3D space 3D space 3D space 3D space 3D space Parameter space Parameter space Parameter space
C'
C Q0 T0 q0 L0
q1 Q1
T1
β0 n0 n1 L1 q2 qi qi p
Unit tangent vector of curve C at Q1 Principal normal vector of C at Q0 Unit normal vector of S at Q0 Unit normal vector of S at Q1 Direction of C ' at q1 Point lying in the neighborhood of q1 on C ' Point calculated by the proposed approach in step i The exact position of qi
Parameter space
Position calculation error of qi
Parameter space
ε1
Error caused by the direction calculation error of qi −1
Parameter space
ε2
The neglected remainder of the Taylor expansion
Parameter space
Qi
Three-dimensional image point of qi on S
3D space
Qi
Three-dimensional image point of qi on S
3D space
Ti −1
Unit tangent vector at Qi −1 .
3D space
Li −1
Direction of C ' at point qi −1 .
Parameter space
Li−1
The exact value of Li −1
Parameter space
Ε il−1
Calculation error of Li −1
Parameter space
Ε i−t 1
Calculation error of Ti−1
3D space
Ti−1
Unit tangent vector at Qi −1
3D space
ni−1
Unit normal vector of S at Qi-1
3D space
β i− 2
Principal normal vector of C at Qi-2.
3D space
β i− 2
Principal normal vector of C at Qi − 2
3D space
ni−1
Unit normal vector of S at Qi −1
3D space
Ti− 2
Unit tangent vector at Qi − 2
3D space
ε3
Local truncation error of qi
Parameter space
Εi
A geometric method for computation of geodesic on parametric surfaces
Abstract This paper proposes a geometric algorithm for computation of geodesic on surfaces. The geodesics on surfaces are traced in a simple way which is independent of the complex description of the geodesic equations. Through derivation process, the calculation error of this algorithm is obtained. A step size adjustment strategy which enables the step size adapt to the geometry of surface is introduced. The proposed method is also compared to some other well-known methods in this study. Many geodesics computed using these approaches on various B-spline surfaces or their equivalent tessellated surfaces have been presented. Experiments demonstrate that the proposed algorithm is efficient. Meanwhile, the results show that the step size adjustment strategy works well for most of the cases.
Key words: geometric; geodesic; geodesic equations; tessellated surface *Corresponding author at: State Key Lab of Digital Manufacturing Equipment and Technology. Huazhong University of Science and Technology, Wuhan 430074, PR China. Tel.: +86 18986129289 E-mail address: [email protected]
1. Introduction Geodesic on a surface is an intrinsic geometric feature that plays an important role in a diversity of applications. The importance of geodesic is well known in the computer simulation of design and manufacturing process of composite components. For tape laying process, geodesic allows minimizing the steering of the tape [1]. For filament winding process, it would be impossible to lay a filament in any way other than along a geodesic on a frictionless convex surface. Geodesic also finds its application in surface and shape processing, such as surface flattening [2, 3], segmentation, sampling, meshing and comparison of shapes [4]. A curve on a surface is called a geodesic if and only if the normal vector to the curve is everywhere parallel to the normal vector of the surface. A geodesic can also be defined as a curve with zero geodesic curvature. Available approaches for the computation of geodesic on surfaces can be classified broadly as analytical and numerical, and the later one is more widely used than the former one recently. The analytical approaches presented by Do Carmo are quite complex and closed form solutions cannot be found for geodesic on general surfaces [5]. As for the numerical approaches, Beck et al computed geodesic paths on a bicubic spline surface by using the fourth order Runge-Kutta method [6]. Patrikalakis and Badris examined geodesic curves on parametric surfaces when they constructed offset curves on Rational B-spline surfaces [7]. Sneyd and Peskin investigated the computation of geodesic paths on a generalized cylinder using a second order Runge-Kutta method [8]. Hotz and Hagen presented a geometric method for the construction of geodesics on arbitrary surface [9]. Their method is based on the fundamental property that geodesics are a generalization of straight lines on planes. The superiority of this method is that it makes us independent of the complex description of the surface. Ying and Candes proposed an approach for rapidly computing a very large number of geodesics on a smooth surface. Their approach is built upon the phase flow method and is especially well suited for the problem of computing creeping
rays [10]. Finding the shortest path between two points on surface is a classical problem and it has many important applications. Maekawa introduced an approach for finding the shortest path between two points on a free-form parametric surface based on a relaxation method relying on finite difference discretization [11]. Kasap et al presented a numerical method for the computation of geodesic between two points on surface [12]. They solve the nonlinear differential equations of geodesic via the finite-difference method and the iterative method. Chen proposed a geodesic-like curve that approaches to the geodesic on surface when the order of the curve reaches to infinity [13]. Chen’s geodesic-like curve has been proved to be accurate. Many accurate discrete methods approach geodesics or the shortest paths on tessellated surfaces [14, 15], polygonal surfaces [16, 17] and triangular meshes [18, 19]. Bose et al did a survey which gave an overview of theoretically and practically relevant algorithms to compute geodesic paths and distances on three-dimensional surfaces [20]. The algorithms are differentiated based on theoretical time complexity and approximation ratio. In general, conventional approaches for computation of geodesic can be classified broadly into three types. First, solve the nonlinear differential equations of geodesic by using a numerical integration method such as Runge-Kutta method. Second, compute discrete geodesics on discrete surfaces by following some rules. Last, construct geodesics on smooth surfaces by using geometric methods. The first type is elegant and accurate, but the differential equations of geodesic are very complicated and generally not easy to solve. Discrete geodesics have been gaining attention as computer becomes increasingly more powerful and discretized models become more prevalent in geometric modeling. However, the discrete geodesics can not be computed directly on the original smooth surface. There is little work on constructing geodesics directly on smooth surfaces by using geometric methods. The geometric method presented by Hotz and Hagen is complicated as it needs to compute projection of point of tangent plane in each step. Ravi Kumar’s method can be extended to estimate geodesics on regular surfaces [13]. However, this method does not perform on regular surfaces directly. The discrete geodesic is computed at first and then projected to the regular surface. In this work, we present a compact geometric method for constructing geodesics directly on smooth surfaces. Firstly, principle of the method is described in detail. Then, calculation error of the method is provided through derivation process. Finally, the method is tested on a group of different surfaces.
2. Principle of the algorithm This work describes an efficient method which aims to trace geodesic on parametric surfaces. The presented approach is independent of the complex description of the geodesic equations and is summarized as follows: We look at a surface S which is given by a parameterization r (u, v) = [ x(u, v), y (u, v), z (u, v)] , where x, y and z are differentiable functions of the parameters u and v. As this work focus on providing a computational method for the geodesic, we assume that S is regular. A curve C lying entirely on the surface S can be expressed in parametric form by: u = u ( s ), v = v( s ) , (1) where s is the arc length. Let Q0 = r (u0 , v0 ) denote a starting point on curve C in three-dimensional space, and T0 the unit tangent vector of curve C at Q0.
Then the three-dimensional directional vector T0 is pulled back to the parametric domain adopting a simple procedure [21], and we get the direction of C ' at q0 = (u0 , v0 ) which can be represented as L0 = (u '(q0 ), v '( q0 )) . Where C ' and q0 are the preimage of curve C and point Q0 in the parametric domain, respectively. Assume that q1 = (u1 , v1 ) is the point lying in the neighborhood of q0 on C ' , and it can be approximated as: (u1 , v1 ) = (u0 + u '( q0 )Δs, v0 + v '( q0 )Δs ) ,
(2)
where Δs is arch length increment. Let Q1 = r (u1 , v1 ) denote the three-dimensional image point of q1 on S. According to the Frenet-Serret formulas, the unit tangent vector of curve C at Q1 can be approximated as: T1 = T0 + k0 β 0 Δs , (3) where T1 is the unit tangent vector of curve C at Q1, k0 is the curvature of curve C at Q0 and β 0 is the principal normal vector of curve C at Q0. Suppose curve C is geodesic on surface S, and then Eq. (3) can be rewritten as: T1 = T0 ± k0 n0 Δs , where n0 is the unit normal vector of surface S at Q0 which is parallel to β 0 and is given by: n0 =
ru ( q0 ) × rv (q0 ) . |ru ( q0 ) × rv (q0 )|
(4)
(5)
ru and rv are the first order parametric derivatives of the surface patch with respect to u and v,
respectively. Then we consider the curvature of curve C at Q0. Let a denote the included angle made by T0 and T1, and n1 the unit normal vector of surface S at Q1 (see Fig. 1).
Figure. 1. The included angle a made by T0 and T1 From the definition of curvature, k0 can be obtained as: k0 = lim | Δs → 0
a |. Δs
(6)
Assume Δs is small, k0 is approached as follows: |sin a| |n1 .T0 | = . (7) Δs Δs To obtain the sign of n1 .T0 in Eq. (7), the neighborhood of Q0 on Curve C is approached by k0 =
using a circular-arc approximation, as illustrated in Fig. 2. In this work, the circular-arc approximation is deemed to be accurate enough as it is of second order accuracy.
n0 n1
n0 ( β 0 )
Q0
n1 C
T1
C
Q0
T0
Q1
T0
Q1
T1
β0
β 0 and n0 are in the same direction
β 0 and n0 are in the opposite direction
Figure. 2. The relationship between T0 and n1 in two cases In Fig. 2, Q0 and Q1 are both located on the circular-arc, and the direction in which the circular-arc is turning is determined by β 0 . As seen in Fig 2, when β 0 and n0 are in the same direction, k0 = −
n1 .T0 . Substitute it into Δs
Eq. (4) then it follows that: T1 = T0 −
(n1 .T0 ) n0 Δs . Δs
Similarly, when β 0 and n0 are in the opposite direction, k0 =
(8) n1 .T0 . Substitute it into Eq. Δs
(4), and then we get a same equation as Eq. (8). In summary, a universal equation is obtained to calculate T1. Next, T1 is pulled back to the parametric domain, thus we get L1, where L1 = (u '( q1 ), v '(q1 )) is the direction of C ' at q1. Next, we can get point q2 (u2 , v2 ) lying in the neighborhood of q1. Similarly, the subsequent geodesic points can be computed until the boundary of the surface is reached or the number of computed points exceeds a maximal point number. The sequence of geodesic points gives a unique parametric geodesic using Q0 as initial value. We are tracing a geodesic on surface S using Q0 and T0 as initial values. The main algorithm flow is described as follows: 1. The three-dimensional directional vector T0 is pulled back to the parametric domain of the surface S, and the direction of curve C ' at q0 is obtained. 2. Based on the direction of curve C ' at q0, the point lying in the neighborhood of q0 on C ' is derived as q1. 3. The unit normal vector of surface S at Q0 and Q1 is obtained as n0 and n1 separately through a simple procedure. 4. After the steps above, the three-dimensional directional vector T1 can be obtained through Eq. (8). 5. Then q1 and T1 are used to compute the subsequent geodesic point. The sequence of geodesic points gives a unique parametric geodesic using Q0 as initial value. The rest of the paper is organized as follows. In section 3, the calculation error of this algorithm is obtained. In section 4, we introduce a strategy which enables the step size adapt to the geometry of surface during the computation process. Section 5 shows the experimental results and section 6 concludes the paper.
3. Error analysis Let qi = (ui , vi ) denote the point calculated by the proposed method in step i (i = 2....m, m ≥ 2) , the exact position of qi is represented as qi = (ui , vi ) . The three-dimensional image point of qi and qi on S is represented as Qi and Qi , respectively. The main focus in this section is to acquire the difference between qi and qi . Q0 is the starting point of the geodesic, so q0 = q0 . As presented in section 2, the position of qi is obtained as: qi = qi −1 + Li −1Δs , where Li −1 = (u '( qi −1 ), v '(qi −1 )) is the direction of C ' at point qi −1 .
(9)
From the Taylor expansion, it follows that: qi = qi −1 + Li −1Δs +
1 M i −1Δs 2 + O ( Δs 3 ) , 2
(10)
where Li −1 = (u '(qi −1 ), v '(qi −1 )) , M i −1 = (u ''(qi −1 ), v ''(qi −1 )) and O () is big O notation. u ' and v ' are the first-order partial derivatives of u(s) and v(s) with respect to s, while u '' and v '' are
the second-order partial derivatives of u(s) and v(s). Combine Eq. (9) with Eq. (10), the position calculation error of qi is obtained as: 1 2
Ε ip = qi −1 − qi −1 + ( Li −1 − Li −1 )Δs + M i −1Δs 2 + O (Δs 3 ) ,
(11)
where the superscript p of Ei specifies the position calculation error. By referring to Eq. (11), Ε ip consists of three parts: 1. The position calculation error of qi −1 :
Ε ip−1 = qi −1 − qi −1 . 2. The error caused by the direction calculation error of qi −1 :
ε1 = Ε il−1Δs , where Ε il−1 = Li −1 − Li −1 . The superscript l of Ei−1 specifies the direction calculation error. 3. The neglected remainder of the Taylor expansion: 1 M i −1Δs 2 + O (Δs 3 ) . 2 In Eq. (11), ε1 + ε 2 is the local truncation error in step i. As can be seen from above, ε 2 is on
ε2 =
the order of O (Δs 2 ) . However, it is not easy to determine the order of ε1 as we do not know the expression of Ε i−l 1 . Ε i−l 1 will affect Ε ip while will be affected by Ε i−p1 , which reveals that there is a coupling relationship between ‘position calculation error’ and ‘direction calculation error’. In the subsequent part, in order to decoupling and acquire Ε ip , we will perform an overall
analysis. First of all, we make an assumption that this algorithm is of kth-order accuracy, where k = 0 or k = 1 . Then Ε ip can be represented as:
Ε ip = ω1Δs k +1 ,
(12)
where ω1 is constant depends on variable i. Next, we are going to show that k = 1 . By referring to section 2, Li-1 is derived by solving the following system of linear equations: ALi −1T = BTi −1T ,
(13)
ª |ru (qi −1 )|2 ru ( qi −1 ).rv (qi −1 ) º ª ru (qi −1 ) º where A = « », B=« » and Ti−1 is the unit tangent vector 2 |rv (qi −1 )| «¬ ru (qi −1 ).rv ( qi −1 ) »¼ ¬ rv (qi −1 ) ¼
at Qi −1 . At the same time, Li−1 can be got by solving a system of linear equations quite similar to Eq.(13): T
T
ALi −1 = BTi −1 ,
ª
ª ru (qi −1 ) º ru ( qi −1 ).rv (qi −1 ) º » , B=« » and 2 |rv (qi −1 )| «¬ ru (qi −1 ).rv ( qi −1 ) »¼ «¬ rv (qi −1 ) »¼
where A = «
|ru (qi −1 )|2
(14) Ti−1
is the unit tangent
vector at Qi −1 . From the Taylor expansion, it follows that: ru ( qi −1 ) = ru ( qi −1 ) + O ( Δs k +1 ) ,
(15)
rv (qi −1 ) = rv (qi −1 ) + O ( Δs k +1 ) .
(16)
From Eq. (13), Eq. (14), Eq. (15) and Eq. (16), through simplification, the following equation is expressed: AΕ il−1T = BΕ it−1T + O (Δs k +1 ) ,
(17)
where Ε it−1 = Ti −1 − Ti −1 . The superscript t of Ei−1 specifies the unit tangent vector calculation error. With Eq.(17), in order to obtain Ε i−l 1 , Ε i−t 1 is investigated. By referring to the Taylor expansion, it follows that: 1 Ti −1 = Ti − 2 +|T '(Qi − 2 )|β i − 2 Δs + T ''(Qi − 2 )Δs 2 + O ( Δs 3 ) . 2
(18)
where Ti− 2 is the unit tangent vector at Qi − 2 . |T '(Qi − 2 )| and β i− 2 are the curvature and
principal normal vector of C at Qi − 2 , respectively. C represents the ideal geodesic with Q0 and T0 as initial values. The unit tangent vector calculation error of Qi −1 is obtained as: |ni −1 .Ti − 2 | 1 β i − 2 )Δs + T ''(Qi − 2 )Δs 2 + O ( Δs 3 ) Δs 2 . |ni −1 .Ti − 2 | 1 =Ti − 2 − Ti − 2 + (|T '(Qi − 2 )| − ) β i − 2 Δs + T ''(Qi − 2 )Δs 2 +O ( Δs k + 2 ) Δs 2
Ε it−1 = Ti − 2 − Ti − 2 + (|T '(Qi − 2 )|β i − 2 −
(19)
where ni−1 is the unit normal vector of S at Qi-1 and β i− 2 is the principal normal vector of C at Qi-2. Move Ti− 2 and Ti−1 to a same starting point O, and draw a unit circle with O as center (see Fig 3). This unit circle will pass through the end points of Ti− 2 and Ti−1 , and the included angle a made by Ti− 2 and Ti−1 is equal to the length of arc MN .
Ti− 2
a Ti−1 Figure. 3. Move Ti− 2 and Ti−1 to a same starting point O As illustrated in Fig. 3, the following equation is obtained: JJJJJG JJJJJG |ni −1 .Ti − 2 | |sin a| |MN | |sin a| |MN | a JJJJJG = = = |cos | , Δs Δs Δs |MN | Δs 2
(20)
a 2 JJJJG Now we consider | MN | :
where we assume |cos | = 1 when ‘ Δs ’ is small.
JJJJG 1 |MN | = |Ti −1 − Ti − 2 | = |T '(Qi − 2 )Δs + T ''(Qi − 2 )Δs 2 + O ( Δs 3 )| . 2
(21)
Substitute Eq. (21) into Eq. (20), the following equation is expressed: |n .T | 1 1 −( |T ''(Qi − 2 )|Δs + O (Δs 2 )) ≤ i −1 i − 2 −|T '(Qi − 2 )| ≤ |T ''(Qi − 2 )|Δs +O( Δs 2 ) . 2 2 Δs
Substitute ni −1 = ni −1 + O (Δs k +1 ) into Eq. (22), the following equation is derived:
(22)
|n .T | 1 1 −( |T ''(Qi − 2 )|Δs + O (Δs k )) ≤ i −1 i − 2 −|T '(Qi − 2 )| ≤ |T ''(Qi − 2 )|Δs + O (Δs k ) . 2 2 Δs
As k = 0 or k = 1 , let |T '(Qi − 2 )|=
(23)
|ni −1 .Ti − 2 | + O (Δs k ) . Substitute it into Eq. (19), we arrive at a Δs
recursion formula:
Ε it−1 = Ti − 2 − Ti − 2 + (|ni −1 .Ti − 2 |−|ni −1 .Ti − 2 |) β i − 2 + O (Δs k +1 ) .
(24)
By referring to the recursion formula, the unit tangent vector calculation error of Qi −1 is obtained as:
Ε it−1 = Ti −1 − Ti −1 = ω 2 Δs k +1 ,
(25)
where ω 2 is constant depends on variable i. With Eq. (25) and Eq. (17), Ε i−l 1 can be derived as follows:
Ε il−1 = ω3 Δs k +1 ,
(26)
where ω3 is constant. Substitute Ε i−l 1 into Eq. (11), Ε ip is obtained as: 1 2
Ε ip = qi −1 − qi −1 + M i −1Δs 2 + ω3 Δs k + 2 + O (Δs 3 ) .
(27)
With Eq.(27), it is acquired that k = 1 . At the same time, Ε ip can be simplified as: 1 (28) 2 By using the recursion formula as Eq. (28), the position calculation error of qi is obtained:
Ε ip = qi −1 − qi −1 + M i −1Δs 2 + O (Δs3 ) .
Ε ip =
1 i ¦ M j −1Δs 2 + O (Δs3 ) . 2 j =1
(29)
Although Eq. (29) is obtained when i ≥ 2 , it can also be applied to the position calculation error of q1.
4. Step size adjustment strategy In this section, a strategy which enables the step size adapt to the geometry of surface is introduced. By referring to the mean value theorem, the local truncation error of qi can be rewritten as: 1 2
1 2
ε 3 = ε1 + ε 2 = ( u ''(δ i −1 )Δs 2 , v ''(ϕi −1 )Δs 2 ) ,
(30)
where both of δ i −1 and ϕi −1 are located between qi −1 and qi on curve C . Let λ denote the upper limit of the local truncation error of qi , we get: (|u ''(δ i −1 )|+|v ''(ϕi −1 )|) Δs 2 ≤ 2λ .
(31)
In each calculation step, we choose Δs =
2λ as the step size. (|u ''(δ i −1 )|+|v ''(ϕi −1 )|)
Then the position calculation error of qi is estimated as: 1 i −1 1 i −1 |Ε ip | = |¦ u ''(δ j )|Δs 2 + |¦ v ''(ϕ j )|Δs 2 ≤ iλ . 2 j =0 2 j =0
(32)
In the subsequent part, the approximate values of u ''( s) and v ''( s) are considered. From the derivation formula of composite function, it follows that: ªu '( s ) º r '( s ) = [ ru rv ] « ». ¬ v '( s ) ¼
(33)
Take the derivative of Eq. (33) on both sides with respect to s, and it follows: ªu ''( s ) º [ ru rv ] « v ''( s) » = M , ¬ ¼
(34)
where M = r ''( s ) − ruu u '( s ) 2 − 2ruv u '( s )v '( s ) − rvv v '( s ) 2 . By left multiply [ ru
rv ]
T
on both sides of Eq. (34), we get: ª| ru |2 « ¬« ru .rv
ru .rv º ªu ''( s ) º ª ru .M º »« ». »=« | rv |2 ¼» ¬ v ''( s ) ¼ ¬ rv .M ¼
(35)
By solving this system of linear equations, we can get u ''( s ) and v ''( s ) . Here is the detail of the step size adjustment strategy: We start with a starting point q0 (u0 , v0 ) and a starting direction T0. Further we have an initial step size h0 and λ . 1. By using the starting direction T0 and initial step size h0, q1 is obtained. ( n .T )n 2. u ''(q0 ) and v ''(q0 ) are obtained by solving Eq. (35), where − 1 0 0 h0
is used to
approach r ''(q0 ) . 3. Change the step size into h1 =
2λ , where h1 is the step size to calculate (| u ''(q0 ) | + | v ''( q0 ) |)
the position of q2. 4. The position of q2 is acquired as (u1 + u '(q1 )h1 , v1 + v '(q1 )h1 ) . Similarly, the subsequent geodesic points (q3, q4…) can be computed. In each step, we need to compute the position of qi and the step size hi for the next step. Although the time costs increases in each step, the overall time costs decreases generally. Remark : The quality of result is influence by step size, curve curvature and parameterization. Therefore, for different parameterization or curve curvature, we adjust the step size to maintain result with good quality.
5. Validation of the algorithm In this section, we perform comparisons between our method and some existing elegant methods. At first, our approach is compared with Beck’s approach on accuracy and efficiency as
both approaches work directly on the original surface. Then, both of the two methods mentioned above are compared with Ravi Kumar’s discrete method. Ravi Kumar’s discrete method had been extended to estimate geodesics on regular surfaces in Chen’s work. This method does not perform on regular surfaces directly. The discrete geodesic is computed at first and then projected to the regular surface. In this work, Liu’s point projection method is adopted [22]. Next, various surfaces are adopted to validate the proposed algorithms. All computer programs of the proposed algorithms are written in MATLAB 7.0 and run on a personal computer [Intel Core (TM) i5 Duo Processor, 2.8GHz, 2GB Memory] to test the accuracy and the efficiency. In the present work, all the geodesics are computed until they reach the boundary of the surface and the computation time is recorded. In order to measure the accuracy of the three kinds of algorithms, we record the distance dist (qco , qid ) between last point qco (u (qco ), v(qco )) of the computed geodesic and last point qid (u (qid ), v(qid )) of the ideal geodesic, where dist (qco , qid ) equals to | u (qco ) − u (qid ) | + | v(qco ) − v( qid ) | . The ideal geodesic on B-spline surface is not readily available, so we use the geodesic computed by using the fourth order Runge-Kutta method with extremely small step size as the ideal geodesic. Now, we will introduce the first two examples, respectively. In both of the two examples, our approach is compared with Beck’s approach on accuracy and efficiency. For the two approaches , we record their computation errors at a group of equal step size and their computation time required to achieve a given level of accuracy, respectively. For the purpose of equity, step size of the two algorithms in both of the examples is fixed. Meanwhile, the accuracy and efficiency of Ravi Kumar’s discrete method are also presented. In the first example, we consider a 2 × 2 B-spline surface. As illustrated in Fig. 4, the geodesics are computed on this surface in three different directions. The results of the comparison are demonstrated in Fig. 5-7. The results of Ravi Kumar’s discrete method are shown in Fig. 8. In this work, we simply record the efficiency of Ravi Kumar’s discrete method in V0(du:dv=1:1) as the results in V1(du:dv=1:2) and V2(du:dv=2:1) share a quite similar tendency as V0.
Figure. 4. Geodesics computed on a 2 × 2 B-spline surface and its equivalent tessellated surface in three different directions
0.1
1E-3 1E-4
Ours Beck's
1E-5
1E-4
1E-3
0.01
Computation time (s)
Error
0.01
1E-6
Ours Beck's
1000
100
10
1
0.1
0.1
1E-6
1E-5
1E-4
Step Size
1E-3
0.01
0.1
Error
(a) Error versus step size (b) Computation time versus error Figure. 5. Results of comparison between our approach and Beck’s approach in V0 for Example 1.
1E-4
1E-5
Ours Beck's
Computation time (s)
1E-3
Accuracy
Our Beck's
1000
0.01
1E-3
0.01
10
1
0.1 1E-7
1E-6 1E-4
100
0.1
1E-6
1E-5
Step size
1E-4
1E-3
0.01
0.1
Error
(a) Error versus step size (b) Computation time versus error Figure. 6. Results of comparison between our approach and Beck’s approach in V1 for Example 1.
Our Beck's
1000
1E-3
1E-4
Our Beck's 1E-5
Computation time (s)
Accuracy
0.01
100
10
1
0.1 1E-4
1E-3
0.01
Step size
0.1
1E-6
1E-5
1E-4
1E-3
0.01
0.1
Error
(a) Error versus step size (b) Computation time versus Error Figure. 7. Results of comparison between our approach and Beck’s approach in V2 for Example 1.
du:dv=1:1 du:dv=1:2 du:dv=2:1
10
Time (s)
Error
0.01
1E-3
R P C
100
1 0.1 0.01
1E-4
0
5000
10000
15000
20000
25000
0
5000
Grid quantity
10000
15000
20000
25000
Grid quantity
(a) Error versus grid quantity (b) Efficiency versus grid quantity in V0 in three directions Figure. 8. Results of Ravi Kumar’s discrete approach for Example 1. Grid quantity: the number of triangular meshes on tessellated surface. R: Reorganization time. The time needed to reorganize data structure of tessellated surface. It does not include the time needed to generate the tessellated surface. P: Projection time. The time needed to project discrete geodesic to original surface. C: Computation time. The time needed to compute discrete geodesic on tessellated surface.
For the 2 × 2 B-spline surface, as can be seen in Fig. 5(b), Fig. 6(b) and Fig. 7(b), although Beck’s approach is more accurate, our method are always able to achieve a given level of accuracy with less computation time, which indicates that our method is more efficient. This is mainly due to that the divergence of accuracy between the two approaches is not significant on the 2 × 2 B-spline surface. In the second example, we consider a 6 × 4 B-spline surface. The geodesics are computed on this surface in three different directions, as illustrated in Fig. 9. The ideal geodesic on B-spline surface is calculated in a same way as the first example. The results of the comparison are presented in Fig. 10-12. The results of Ravi Kumar’s discrete method are shown in Fig. 13.
Figure. 9. Geodesics computed on a 6 × 4 B-spline surface and its equivalent tessellated surface in three different directions
0.01
1E-4
1E-5
Our Beck's
Computation time (s)
1E-3
Error
Our Beck's
1000
1E-6
100
10
1 1E-3
0.01
0.1
1E-6
1
1E-5
Step size
1E-4
1E-3
Error
(a) Error versus step size (b) Computation time versus error Figure. 10. Results of comparison between our approach and Beck’s approach in V0 for Example 2.
1000
Error
1E-3 1E-4 1E-5
Our Beck's
1E-6 1E-3
0.01
0.1
Computation time (s)
0.01
Ours Beck's
100
10
1
1
1E-6
1E-5
Step size
1E-4
1E-3
0.01
Error
(a) Error versus step size (b) Computation time versus error Figure. 11. Results of comparison between our approach and Beck’s approach in V1 for Example 2.
0.01
1E-4
1E-5
Our Beck's 1E-3
0.01
0.1
Step length
1
Computation time (s)
Error
1E-3
1E-6
Our Beck's
1000
100
10
1
0.1
1E-6
1E-5
1E-4
1E-3
0.01
Accuracy
(a) Error versus step size (b) Computation time versus error Figure. 12. Results of comparison between our approach and Beck’s approach in V2 for Example 2.
du:dv=1:1 du:dv=1:2 du:dv=2:1
10
Error
Time (s)
1E-3
R P C
100
1 0.1 0.01
1E-4
0
5000
10000
15000
Grid quantity
20000
25000
0
5000
10000
15000
20000
25000
Grid quantity
(a) Error versus grid quantity (b) Efficiency versus grid quantity in V0 in three directions Figure. 13. Results of Ravi Kumar’s discrete approach for Example 2. The difference on efficiency of the two kinds of algorithm for the 6 × 4 B-spline surface is clearly illustrated in Fig. 10(b), Fig. 11(b) and Fig. 12(b). As shown in these diagrams, the comparison results of the two algorithms on efficiency in V0 and V1 depend on the level of the specified accuracy. To be more specifically, our approach is more efficient when the specified computation error is greater than 5 × 10-5 while Beck’s approach is more efficient when the computation error is less than 1× 10-5 . This is because Beck’s method is much more accurate than our approach with the equal step size when the specified computation error is less than 1× 10-5 . In this case, to obtain equivalent accuracy as Beck’s method, our approach generally have to adopt very small step size, which is at the expense of computational efficiency. The comparison results in V2 demonstrate that the divergence of accuracy between the two kinds of approaches is not significant and our method is more efficient, which is quite different from the results in V0 and V1. This is mainly due to the improvement of accuracy of our approach in V2, which can be seen when we compare Fig. 12(a) with Fig. 10(a) and Fig. 11(a). As shown Fig. 5-7 and Fig. 10-12, when the step size decreases, the errors generally decrease. Meanwhile, the time costs increase. For each step size, the errors of Beck’s algorithm are always smaller than those of our method, and this is more obvious in example 2, which indicates that Beck’s method is more accurate. However, as for the efficiency of the two kinds of algorithms about achieving a given level of accuracy, the results are different. Since our approach is less accurate, to acquire equivalent accuracy as Beck’s method, our approach have to adopt smaller step size. The step size of our approach depends on the level of divergence of accuracy between the two approaches. If the accuracy of Beck’s method is significantly better than our approach, the latter method might have to use very small step size, which will increase the time cost largely. It is well known that Beck’s method is computationally intensive. In this work, we have found that the time costs of Beck’s method is approximately 7 times larger than that of our approach with the same step size. Then, we have an estimation that our algorithm will be more efficient when the step size of our approach is bigger than 1/7 of the step size of Beck’s approach. The computational results of Ravi Kumar’s discrete method are presented in Fig. 8 and Fig. 13. As shown in these figures, the accuracy and time costs of Ravi Kumar’s discrete method are influenced by the grid quantity instead of the step size. With the increase of the grid quantity, the tessellated surface gets closer to the original surface, and the computation errors generally
decrease. However, the reorganization time increase obviously. The results show that the time needed to compute geodesic on tessellated surface is much less than that needed to reorganize data structure of the tessellated surface. The reorganization process is computationally expensive, so Ravi Kumar’s discrete method is more suitable for computation of massive geodesics on one surface. In this case, the average reorganization time in each computation is reduced largely. As shown in Fig. 5-8 and Fig. 10-13 , all of the three algorithms are able to compute geodesic on surface. Both of our algorithm and Beck’s method have better accuracy than Ravi Kumar’s discrete method as the former two directly use information of the original surface. Experimental results show that our algorithm is efficient, and the geodesics which are completely on the surfaces can be obtained. Meanwhile, we do not need to solve the differential equations of geodesic which are quite complicated. Next, in order to validate the effectiveness of our step size adjustment strategy, we will introduce another two examples. For convenience, our approach with step size adjustment strategy is called variable step algorithm, and with fixed step size is called fixed step algorithm. Then, variable step algorithm is compared with fixed step algorithm on accuracy and efficiency. The step size of the fixed step algorithm is determined by the minimum step size of the variable step algorithm. In the third example, the geodesics are computed on the surface given by the first example in three different directions. The results of the comparison are recorded in Table 1. Table 1 Results of comparison for Example 3. q0=(0, 0), V0(du:dv=1:1), λ = 1e − 5 Type
Fixed step algorithm
Variable step algorithm
Step size ( Δs )
0.0104 1.41557e-03 1.417 1242
max=0.0538; min=0.0104 2.07882e-03 0.875 419
Type
Fixed step algorithm
Variable step algorithm
Step size ( Δs )
0.0106 1.15169e-03 1.306 1211
max= 0.0618; min=0.0106 4.27124e-03 0.797 390
Type
Fixed step algorithm
Variable step algorithm
Step size ( Δs )
0.0112 2.36930e-03 1.391 1199
max= 0.0486; min=0.0112 5.17535e-03 0.953 467
Errors Time (s) Number of discrete points q0=(0, 0), V1(du:dv=1:2), λ = 1e − 5
Errors Time (s) Number of discrete points q0=(0, 0), V2(du:dv=2:1), λ = 1e − 5
Errors Time (s) Number of discrete points
As shown in Table 1, the accuracy of the fixed step algorithm is better than the variable step algorithm. This is due to the fact that the step size of the former algorithm is always not greater than the latter. However, as for the computational efficiency, the variable step algorithm has superiority. Furthermore, the number of discrete points on the geodesics computed by the variable step algorithm is smaller than half of that of the fixed step algorithm, which means the geodesics
computed by the former algorithm are more compact. In the fourth example, we consider another 2 × 2 B-spline. The geodesics are computed on this surface in three different directions, as illustrated in Fig. 14. The results of the comparison are recorded in Table 2.
Figure. 14. Geodesics computed on a 2 × 2 B-spline surface in three different directions Table 2 Results of comparison for Example 4. q0=(0, 0), V0(du:dv=1:1), λ = 1e − 5 Type
Fixed step algorithm
Variable step algorithm
Step size ( Δs )
0.0294 3.85031e-03 0.641 507
max=0.0800; min=0.0294 5.66168e-03 0.765 366
Type
Fixed step algorithm
Variable step algorithm
Step size ( Δs )
0.0162 2.47817e-03 0.828 690
max= 0.0843; min=0.0162 2.97750e-03 0.623 297
Type
Fixed step algorithm
Variable step algorithm
Step size ( Δs )
0.0196 2.27483e-03 0.797 677
max= 0.133; min=0.0196 3.39647e-03 0.625 297
Errors Time (s) Number of discrete points q0=(0, 0), V1(du:dv=1:2), λ = 1e − 5
Errors Time (s) Number of discrete points q0=(0, 0), V2(du:dv=2:1), λ = 1e − 5
Errors Time (s) Number of discrete points
As shown in Table 2, the comparison results of Example 4 in direction V1 and direction V2 are quite similar to the results of Example 3. However, as for direction V0, it takes the variable step algorithm more time on computational process even though the number of steps decreases, which means the reduced time costs due to the decrease of step numbers does not make up for the sum of the increased time costs in each step. As illustrated in Table 1 and Table 2, the average time costs in each step of the fixed step algorithm and the variable step algorithm are 1.2ms and 2.1ms, respectively. According to this, we have an estimate that the variable step algorithm is more efficient than the fixed step algorithm when the number of discrete points computed by the former
algorithm is smaller than 4/7 of that of the latter. Generally speaking, our variable step algorithm works well for most of the cases. The location of discrete points obtained by this algorithm distributes reasonably according to the geometry of the surface, and this could contribute to outstanding reduction of the number of discrete points on the geodesic, which means the geodesic computed by the variable step algorithm is more compact.
6. Conclusion This paper presents a geometric method for computation of geodesic on parametric surfaces. The geodesic which is completely on the surface can be obtained by this algorithm. Furthermore, we do not need to solve the differential equations of geodesic which are quite complicated and generally not easy to solve. Experimental results indicate that our approach is more efficient than Beck’s method on a 2 × 2 B-spline surface while the efficiency of our approach on a 6 × 4 B-spline surface depends on the tracing direction and the level of the specified accuracy. Generally speaking, our approach is more suitable in two cases. First, the parametric surface is relative flat. Second, the requirement for the accuracy is modest. In these cases, the accuracy of our algorithm is guaranteed. At the same time, our approach will trace the geodesic efficiently as we do not need to use quite small step size. One possible application of our method is path planning in the tape laying process. In this process, the composite tapes generally follow geodesic path on a mould surface which is nearly flat or nearly conventional cylinderical. Comparison results show that both of our algorithm and Beck’s method have better accuracy than Ravi Kumar’s discrete method. At the same time, Ravi Kumar’s discrete method is more suitable for computation of massive geodesics on one surface. The results also demonstrate that our step size adjustment strategy works well for most of the cases. The step size adjustment strategy has contributed to outstanding reduction of the number of discrete points on the geodesic, which means the geodesic computed with the step size adjustment strategy is more compact. However, there are also limitations in our algorithm. The function of the parameterized surface must be given in advance, which means this algorithm does not apply to general surfaces. In addition, the surface should contain no singular point. This is due to that we need to pull the direction of the curve back to the parametric domain in each step.
Acknowledgements This work was supported by the National Basic Research Program of China (2013CB035805) and the National Science and Technology Major Project of China (2010ZX04016-013).
References [1] Debout, P., Chanal, H., Duc, E., 2011. Tool path smoothing of a redundant machine: application to automated fiber placement. Computer-Aided Design 43 (2), 122–132. [2]
Azariadis, P., Aspragathos, N.A., 1997. Design of plane developments of doubly curved surfaces. Computer-Aided Design. 29 (10), 675–685.
[3]
Azariadis, P., Aspragathos, N.A., 2001. Geodesic curvature preservation in surface flattening through
constrained global optimization. Computer-Aided Design 33, 581–591. [4] Peyré, G. Péchaud, M. Keriven, R. Cohen, L., 2010. Geodesic Methods in Computer Vision and Graphics. In Found. Trends. Comput. Graph. Vis, 197-397. [5] Do Carmo, M.P., 1976. Differential Geometry of Curves and Surfaces. Prentice-Hall, Englewood Cliffs, NJ. [6] Beck, J.M., Farouki, R.T., Hinds, J.K., 1986. Surface analysis methods. IEEE Computer Graphics and Applications 6, 18–36. [7] Patrikalakis, N.W., Badris, L., 1989. Offsets of curves on rational B-spline surfaces. Engineering with computers 5, 39–46. [8] Sneyd, J., Peskin, C.S., 1990. Computation of geodesic trajectories on tubular surfaces. SIAM Journal of Scientific Statistical Computing 11, 230–241 [9] Hotz, I., Hagen, H., 2000. Visualizing geodesics. In: Proceedings IEEE Visualization, Salt Lake City, UT, PP. 311–318 [10] Maekawa, T., 1996. Computation of shortest paths on free-from parametric surfaces. Journal of Mechanical Design, ASME Transactions 118 (4), 499–508. [11]
Ying, L.X., Candes, E.j., 2006. Fast geodesics computation with the phase flow method. Journal of Computational Physicsm 220 (1), 6-18.
[12] Kasap, E., Yapici, M., Akyildiz, F.T., 2005. A numerical study for computation of geodesic curves. Applied Mathematics and Computation 171 (2), 1206–1213. [13] Chen, S.-G., 2010. Geodesic-like curves on parametric surfaces. Computer Aided Geometric Design 27, 106–117. [14] Tucker, C.L., 1997. Forming of advanced composites. In: Gutowski, T.G. (Ed.), Advanced Composites Manufacturing. Wiley, New York. [15] Ravi Kumar, G.V.V., Prabha Srinivasan., Devaraja Holla, V., Shastry, K.G., Prakash, B.G., 2003. Geodesic curve computations on surfaces. Computer Aided Geometric Design 20, 119–133. [16] Polthier, K., Schmies, M., 1998. In: Hege, H.C., Polthier, H.K. (EDS.), Straightest Geodesics On Polyhedral Surfaces in Mathematical Visualization. Springer-Verlag, Berlin. [17]
Kanai, T., Suzuki, H., 2001. Approximate shortest path on a polyhedral surface and its applications, Computer Aided Design 33 (11), 801–811.
[18] Martinez, D., Velho, L., Carvalho, P.C., 2005. Computing geodesics on triangular meshes. Computer & Graphics 29, 667–675. [19] Surazhsky, V., Surazhsky, T., Kirsanov, D., Gortler, S., Hoppe, H., 2005. Fast exact and approximate geodesics on meshes. In: Proc. of SIGGRAPH 2005. ACM Transactions on Graphics 24 (3), 553–560. [20]
Bose, P., Maheshwari, A., Shu, C., Wuhrer, S., 2011. A survey of geodesic paths on 3D surfaces. Computational Geometry: Theory and Applications 44 (9), 486–498.
[21] Piegl, L., Tiller, W., 1997. The NURBS book, 2nd Edition. Springer, New York. [22] Liu, X,-M., Yang, L., Yong, J.-H., Gu, H,-J., Sun, J,-G., 2009. A torus patch approximation approach for point projection on surfaces. Computer Aided Geometric Design 26 (5), 593–8.
Peng Zhang, Ronglei Sun, Tao Huang
PII: DOI: Reference:
S0167-8396(15)00096-5 http://dx.doi.org/10.1016/j.cagd.2015.08.001 COMAID 1502
To appear in:
Computer Aided Geometric Design
Received date: Revised date: Accepted date:
17 November 2014 8 August 2015 8 August 2015
Please cite this article in press as: Zhang, P., et al. A geometric method for computation of geodesic on parametric surfaces. Comput. Aided Geom. Des. (2015), http://dx.doi.org/10.1016/j.cagd.2015.08.001
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Highlights • • • •
An efficient geometric algorithm aims to trace geodesic on parametric surfaces. Independent of the complex description of the geodesic equations. Simpler and faster than the existing geometric method. Step size of the algorithm adapts to the geometry of parametric surface.
A geometric method for computation of geodesic on parametric surfaces Peng Zhang, Ronglei Sun*, Tao Huang School of Mechanical Science and Engineering, State Key Laboratory of Digital Manufacturing Equipment and Technology, Huazhong University of Science and Technology, Wuhan, 430074, China *Corresponding author. Tel.: +86 27 87559416. E-mail address: [email protected]
A sketch of used notations Notation
Description
Used space
S C
Parametric surface Geodesic curve on S Preimage of C in parametric domain Ideal geodesic curve on S Starting point on curve C Unit tangent vector of C at Q0 Preimage of Q0 in parametric domain Direction of C ' at q0 Point lying in the neighborhood of q0 on C ' Three-dimensional image point of q1 on S
3D space 3D space Parameter space 3D space 3D space 3D space Parameter space Parameter space Parameter space 3D space 3D space 3D space 3D space 3D space Parameter space Parameter space Parameter space
C'
C Q0 T0 q0 L0
q1 Q1
T1
β0 n0 n1 L1 q2 qi qi p
Unit tangent vector of curve C at Q1 Principal normal vector of C at Q0 Unit normal vector of S at Q0 Unit normal vector of S at Q1 Direction of C ' at q1 Point lying in the neighborhood of q1 on C ' Point calculated by the proposed approach in step i The exact position of qi
Parameter space
Position calculation error of qi
Parameter space
ε1
Error caused by the direction calculation error of qi −1
Parameter space
ε2
The neglected remainder of the Taylor expansion
Parameter space
Qi
Three-dimensional image point of qi on S
3D space
Qi
Three-dimensional image point of qi on S
3D space
Ti −1
Unit tangent vector at Qi −1 .
3D space
Li −1
Direction of C ' at point qi −1 .
Parameter space
Li−1
The exact value of Li −1
Parameter space
Ε il−1
Calculation error of Li −1
Parameter space
Ε i−t 1
Calculation error of Ti−1
3D space
Ti−1
Unit tangent vector at Qi −1
3D space
ni−1
Unit normal vector of S at Qi-1
3D space
β i− 2
Principal normal vector of C at Qi-2.
3D space
β i− 2
Principal normal vector of C at Qi − 2
3D space
ni−1
Unit normal vector of S at Qi −1
3D space
Ti− 2
Unit tangent vector at Qi − 2
3D space
ε3
Local truncation error of qi
Parameter space
Εi
A geometric method for computation of geodesic on parametric surfaces
Abstract This paper proposes a geometric algorithm for computation of geodesic on surfaces. The geodesics on surfaces are traced in a simple way which is independent of the complex description of the geodesic equations. Through derivation process, the calculation error of this algorithm is obtained. A step size adjustment strategy which enables the step size adapt to the geometry of surface is introduced. The proposed method is also compared to some other well-known methods in this study. Many geodesics computed using these approaches on various B-spline surfaces or their equivalent tessellated surfaces have been presented. Experiments demonstrate that the proposed algorithm is efficient. Meanwhile, the results show that the step size adjustment strategy works well for most of the cases.
Key words: geometric; geodesic; geodesic equations; tessellated surface *Corresponding author at: State Key Lab of Digital Manufacturing Equipment and Technology. Huazhong University of Science and Technology, Wuhan 430074, PR China. Tel.: +86 18986129289 E-mail address: [email protected]
1. Introduction Geodesic on a surface is an intrinsic geometric feature that plays an important role in a diversity of applications. The importance of geodesic is well known in the computer simulation of design and manufacturing process of composite components. For tape laying process, geodesic allows minimizing the steering of the tape [1]. For filament winding process, it would be impossible to lay a filament in any way other than along a geodesic on a frictionless convex surface. Geodesic also finds its application in surface and shape processing, such as surface flattening [2, 3], segmentation, sampling, meshing and comparison of shapes [4]. A curve on a surface is called a geodesic if and only if the normal vector to the curve is everywhere parallel to the normal vector of the surface. A geodesic can also be defined as a curve with zero geodesic curvature. Available approaches for the computation of geodesic on surfaces can be classified broadly as analytical and numerical, and the later one is more widely used than the former one recently. The analytical approaches presented by Do Carmo are quite complex and closed form solutions cannot be found for geodesic on general surfaces [5]. As for the numerical approaches, Beck et al computed geodesic paths on a bicubic spline surface by using the fourth order Runge-Kutta method [6]. Patrikalakis and Badris examined geodesic curves on parametric surfaces when they constructed offset curves on Rational B-spline surfaces [7]. Sneyd and Peskin investigated the computation of geodesic paths on a generalized cylinder using a second order Runge-Kutta method [8]. Hotz and Hagen presented a geometric method for the construction of geodesics on arbitrary surface [9]. Their method is based on the fundamental property that geodesics are a generalization of straight lines on planes. The superiority of this method is that it makes us independent of the complex description of the surface. Ying and Candes proposed an approach for rapidly computing a very large number of geodesics on a smooth surface. Their approach is built upon the phase flow method and is especially well suited for the problem of computing creeping
rays [10]. Finding the shortest path between two points on surface is a classical problem and it has many important applications. Maekawa introduced an approach for finding the shortest path between two points on a free-form parametric surface based on a relaxation method relying on finite difference discretization [11]. Kasap et al presented a numerical method for the computation of geodesic between two points on surface [12]. They solve the nonlinear differential equations of geodesic via the finite-difference method and the iterative method. Chen proposed a geodesic-like curve that approaches to the geodesic on surface when the order of the curve reaches to infinity [13]. Chen’s geodesic-like curve has been proved to be accurate. Many accurate discrete methods approach geodesics or the shortest paths on tessellated surfaces [14, 15], polygonal surfaces [16, 17] and triangular meshes [18, 19]. Bose et al did a survey which gave an overview of theoretically and practically relevant algorithms to compute geodesic paths and distances on three-dimensional surfaces [20]. The algorithms are differentiated based on theoretical time complexity and approximation ratio. In general, conventional approaches for computation of geodesic can be classified broadly into three types. First, solve the nonlinear differential equations of geodesic by using a numerical integration method such as Runge-Kutta method. Second, compute discrete geodesics on discrete surfaces by following some rules. Last, construct geodesics on smooth surfaces by using geometric methods. The first type is elegant and accurate, but the differential equations of geodesic are very complicated and generally not easy to solve. Discrete geodesics have been gaining attention as computer becomes increasingly more powerful and discretized models become more prevalent in geometric modeling. However, the discrete geodesics can not be computed directly on the original smooth surface. There is little work on constructing geodesics directly on smooth surfaces by using geometric methods. The geometric method presented by Hotz and Hagen is complicated as it needs to compute projection of point of tangent plane in each step. Ravi Kumar’s method can be extended to estimate geodesics on regular surfaces [13]. However, this method does not perform on regular surfaces directly. The discrete geodesic is computed at first and then projected to the regular surface. In this work, we present a compact geometric method for constructing geodesics directly on smooth surfaces. Firstly, principle of the method is described in detail. Then, calculation error of the method is provided through derivation process. Finally, the method is tested on a group of different surfaces.
2. Principle of the algorithm This work describes an efficient method which aims to trace geodesic on parametric surfaces. The presented approach is independent of the complex description of the geodesic equations and is summarized as follows: We look at a surface S which is given by a parameterization r (u, v) = [ x(u, v), y (u, v), z (u, v)] , where x, y and z are differentiable functions of the parameters u and v. As this work focus on providing a computational method for the geodesic, we assume that S is regular. A curve C lying entirely on the surface S can be expressed in parametric form by: u = u ( s ), v = v( s ) , (1) where s is the arc length. Let Q0 = r (u0 , v0 ) denote a starting point on curve C in three-dimensional space, and T0 the unit tangent vector of curve C at Q0.
Then the three-dimensional directional vector T0 is pulled back to the parametric domain adopting a simple procedure [21], and we get the direction of C ' at q0 = (u0 , v0 ) which can be represented as L0 = (u '(q0 ), v '( q0 )) . Where C ' and q0 are the preimage of curve C and point Q0 in the parametric domain, respectively. Assume that q1 = (u1 , v1 ) is the point lying in the neighborhood of q0 on C ' , and it can be approximated as: (u1 , v1 ) = (u0 + u '( q0 )Δs, v0 + v '( q0 )Δs ) ,
(2)
where Δs is arch length increment. Let Q1 = r (u1 , v1 ) denote the three-dimensional image point of q1 on S. According to the Frenet-Serret formulas, the unit tangent vector of curve C at Q1 can be approximated as: T1 = T0 + k0 β 0 Δs , (3) where T1 is the unit tangent vector of curve C at Q1, k0 is the curvature of curve C at Q0 and β 0 is the principal normal vector of curve C at Q0. Suppose curve C is geodesic on surface S, and then Eq. (3) can be rewritten as: T1 = T0 ± k0 n0 Δs , where n0 is the unit normal vector of surface S at Q0 which is parallel to β 0 and is given by: n0 =
ru ( q0 ) × rv (q0 ) . |ru ( q0 ) × rv (q0 )|
(4)
(5)
ru and rv are the first order parametric derivatives of the surface patch with respect to u and v,
respectively. Then we consider the curvature of curve C at Q0. Let a denote the included angle made by T0 and T1, and n1 the unit normal vector of surface S at Q1 (see Fig. 1).
Figure. 1. The included angle a made by T0 and T1 From the definition of curvature, k0 can be obtained as: k0 = lim | Δs → 0
a |. Δs
(6)
Assume Δs is small, k0 is approached as follows: |sin a| |n1 .T0 | = . (7) Δs Δs To obtain the sign of n1 .T0 in Eq. (7), the neighborhood of Q0 on Curve C is approached by k0 =
using a circular-arc approximation, as illustrated in Fig. 2. In this work, the circular-arc approximation is deemed to be accurate enough as it is of second order accuracy.
n0 n1
n0 ( β 0 )
Q0
n1 C
T1
C
Q0
T0
Q1
T0
Q1
T1
β0
β 0 and n0 are in the same direction
β 0 and n0 are in the opposite direction
Figure. 2. The relationship between T0 and n1 in two cases In Fig. 2, Q0 and Q1 are both located on the circular-arc, and the direction in which the circular-arc is turning is determined by β 0 . As seen in Fig 2, when β 0 and n0 are in the same direction, k0 = −
n1 .T0 . Substitute it into Δs
Eq. (4) then it follows that: T1 = T0 −
(n1 .T0 ) n0 Δs . Δs
Similarly, when β 0 and n0 are in the opposite direction, k0 =
(8) n1 .T0 . Substitute it into Eq. Δs
(4), and then we get a same equation as Eq. (8). In summary, a universal equation is obtained to calculate T1. Next, T1 is pulled back to the parametric domain, thus we get L1, where L1 = (u '( q1 ), v '(q1 )) is the direction of C ' at q1. Next, we can get point q2 (u2 , v2 ) lying in the neighborhood of q1. Similarly, the subsequent geodesic points can be computed until the boundary of the surface is reached or the number of computed points exceeds a maximal point number. The sequence of geodesic points gives a unique parametric geodesic using Q0 as initial value. We are tracing a geodesic on surface S using Q0 and T0 as initial values. The main algorithm flow is described as follows: 1. The three-dimensional directional vector T0 is pulled back to the parametric domain of the surface S, and the direction of curve C ' at q0 is obtained. 2. Based on the direction of curve C ' at q0, the point lying in the neighborhood of q0 on C ' is derived as q1. 3. The unit normal vector of surface S at Q0 and Q1 is obtained as n0 and n1 separately through a simple procedure. 4. After the steps above, the three-dimensional directional vector T1 can be obtained through Eq. (8). 5. Then q1 and T1 are used to compute the subsequent geodesic point. The sequence of geodesic points gives a unique parametric geodesic using Q0 as initial value. The rest of the paper is organized as follows. In section 3, the calculation error of this algorithm is obtained. In section 4, we introduce a strategy which enables the step size adapt to the geometry of surface during the computation process. Section 5 shows the experimental results and section 6 concludes the paper.
3. Error analysis Let qi = (ui , vi ) denote the point calculated by the proposed method in step i (i = 2....m, m ≥ 2) , the exact position of qi is represented as qi = (ui , vi ) . The three-dimensional image point of qi and qi on S is represented as Qi and Qi , respectively. The main focus in this section is to acquire the difference between qi and qi . Q0 is the starting point of the geodesic, so q0 = q0 . As presented in section 2, the position of qi is obtained as: qi = qi −1 + Li −1Δs , where Li −1 = (u '( qi −1 ), v '(qi −1 )) is the direction of C ' at point qi −1 .
(9)
From the Taylor expansion, it follows that: qi = qi −1 + Li −1Δs +
1 M i −1Δs 2 + O ( Δs 3 ) , 2
(10)
where Li −1 = (u '(qi −1 ), v '(qi −1 )) , M i −1 = (u ''(qi −1 ), v ''(qi −1 )) and O () is big O notation. u ' and v ' are the first-order partial derivatives of u(s) and v(s) with respect to s, while u '' and v '' are
the second-order partial derivatives of u(s) and v(s). Combine Eq. (9) with Eq. (10), the position calculation error of qi is obtained as: 1 2
Ε ip = qi −1 − qi −1 + ( Li −1 − Li −1 )Δs + M i −1Δs 2 + O (Δs 3 ) ,
(11)
where the superscript p of Ei specifies the position calculation error. By referring to Eq. (11), Ε ip consists of three parts: 1. The position calculation error of qi −1 :
Ε ip−1 = qi −1 − qi −1 . 2. The error caused by the direction calculation error of qi −1 :
ε1 = Ε il−1Δs , where Ε il−1 = Li −1 − Li −1 . The superscript l of Ei−1 specifies the direction calculation error. 3. The neglected remainder of the Taylor expansion: 1 M i −1Δs 2 + O (Δs 3 ) . 2 In Eq. (11), ε1 + ε 2 is the local truncation error in step i. As can be seen from above, ε 2 is on
ε2 =
the order of O (Δs 2 ) . However, it is not easy to determine the order of ε1 as we do not know the expression of Ε i−l 1 . Ε i−l 1 will affect Ε ip while will be affected by Ε i−p1 , which reveals that there is a coupling relationship between ‘position calculation error’ and ‘direction calculation error’. In the subsequent part, in order to decoupling and acquire Ε ip , we will perform an overall
analysis. First of all, we make an assumption that this algorithm is of kth-order accuracy, where k = 0 or k = 1 . Then Ε ip can be represented as:
Ε ip = ω1Δs k +1 ,
(12)
where ω1 is constant depends on variable i. Next, we are going to show that k = 1 . By referring to section 2, Li-1 is derived by solving the following system of linear equations: ALi −1T = BTi −1T ,
(13)
ª |ru (qi −1 )|2 ru ( qi −1 ).rv (qi −1 ) º ª ru (qi −1 ) º where A = « », B=« » and Ti−1 is the unit tangent vector 2 |rv (qi −1 )| «¬ ru (qi −1 ).rv ( qi −1 ) »¼ ¬ rv (qi −1 ) ¼
at Qi −1 . At the same time, Li−1 can be got by solving a system of linear equations quite similar to Eq.(13): T
T
ALi −1 = BTi −1 ,
ª
ª ru (qi −1 ) º ru ( qi −1 ).rv (qi −1 ) º » , B=« » and 2 |rv (qi −1 )| «¬ ru (qi −1 ).rv ( qi −1 ) »¼ «¬ rv (qi −1 ) »¼
where A = «
|ru (qi −1 )|2
(14) Ti−1
is the unit tangent
vector at Qi −1 . From the Taylor expansion, it follows that: ru ( qi −1 ) = ru ( qi −1 ) + O ( Δs k +1 ) ,
(15)
rv (qi −1 ) = rv (qi −1 ) + O ( Δs k +1 ) .
(16)
From Eq. (13), Eq. (14), Eq. (15) and Eq. (16), through simplification, the following equation is expressed: AΕ il−1T = BΕ it−1T + O (Δs k +1 ) ,
(17)
where Ε it−1 = Ti −1 − Ti −1 . The superscript t of Ei−1 specifies the unit tangent vector calculation error. With Eq.(17), in order to obtain Ε i−l 1 , Ε i−t 1 is investigated. By referring to the Taylor expansion, it follows that: 1 Ti −1 = Ti − 2 +|T '(Qi − 2 )|β i − 2 Δs + T ''(Qi − 2 )Δs 2 + O ( Δs 3 ) . 2
(18)
where Ti− 2 is the unit tangent vector at Qi − 2 . |T '(Qi − 2 )| and β i− 2 are the curvature and
principal normal vector of C at Qi − 2 , respectively. C represents the ideal geodesic with Q0 and T0 as initial values. The unit tangent vector calculation error of Qi −1 is obtained as: |ni −1 .Ti − 2 | 1 β i − 2 )Δs + T ''(Qi − 2 )Δs 2 + O ( Δs 3 ) Δs 2 . |ni −1 .Ti − 2 | 1 =Ti − 2 − Ti − 2 + (|T '(Qi − 2 )| − ) β i − 2 Δs + T ''(Qi − 2 )Δs 2 +O ( Δs k + 2 ) Δs 2
Ε it−1 = Ti − 2 − Ti − 2 + (|T '(Qi − 2 )|β i − 2 −
(19)
where ni−1 is the unit normal vector of S at Qi-1 and β i− 2 is the principal normal vector of C at Qi-2. Move Ti− 2 and Ti−1 to a same starting point O, and draw a unit circle with O as center (see Fig 3). This unit circle will pass through the end points of Ti− 2 and Ti−1 , and the included angle a made by Ti− 2 and Ti−1 is equal to the length of arc MN .
Ti− 2
a Ti−1 Figure. 3. Move Ti− 2 and Ti−1 to a same starting point O As illustrated in Fig. 3, the following equation is obtained: JJJJJG JJJJJG |ni −1 .Ti − 2 | |sin a| |MN | |sin a| |MN | a JJJJJG = = = |cos | , Δs Δs Δs |MN | Δs 2
(20)
a 2 JJJJG Now we consider | MN | :
where we assume |cos | = 1 when ‘ Δs ’ is small.
JJJJG 1 |MN | = |Ti −1 − Ti − 2 | = |T '(Qi − 2 )Δs + T ''(Qi − 2 )Δs 2 + O ( Δs 3 )| . 2
(21)
Substitute Eq. (21) into Eq. (20), the following equation is expressed: |n .T | 1 1 −( |T ''(Qi − 2 )|Δs + O (Δs 2 )) ≤ i −1 i − 2 −|T '(Qi − 2 )| ≤ |T ''(Qi − 2 )|Δs +O( Δs 2 ) . 2 2 Δs
Substitute ni −1 = ni −1 + O (Δs k +1 ) into Eq. (22), the following equation is derived:
(22)
|n .T | 1 1 −( |T ''(Qi − 2 )|Δs + O (Δs k )) ≤ i −1 i − 2 −|T '(Qi − 2 )| ≤ |T ''(Qi − 2 )|Δs + O (Δs k ) . 2 2 Δs
As k = 0 or k = 1 , let |T '(Qi − 2 )|=
(23)
|ni −1 .Ti − 2 | + O (Δs k ) . Substitute it into Eq. (19), we arrive at a Δs
recursion formula:
Ε it−1 = Ti − 2 − Ti − 2 + (|ni −1 .Ti − 2 |−|ni −1 .Ti − 2 |) β i − 2 + O (Δs k +1 ) .
(24)
By referring to the recursion formula, the unit tangent vector calculation error of Qi −1 is obtained as:
Ε it−1 = Ti −1 − Ti −1 = ω 2 Δs k +1 ,
(25)
where ω 2 is constant depends on variable i. With Eq. (25) and Eq. (17), Ε i−l 1 can be derived as follows:
Ε il−1 = ω3 Δs k +1 ,
(26)
where ω3 is constant. Substitute Ε i−l 1 into Eq. (11), Ε ip is obtained as: 1 2
Ε ip = qi −1 − qi −1 + M i −1Δs 2 + ω3 Δs k + 2 + O (Δs 3 ) .
(27)
With Eq.(27), it is acquired that k = 1 . At the same time, Ε ip can be simplified as: 1 (28) 2 By using the recursion formula as Eq. (28), the position calculation error of qi is obtained:
Ε ip = qi −1 − qi −1 + M i −1Δs 2 + O (Δs3 ) .
Ε ip =
1 i ¦ M j −1Δs 2 + O (Δs3 ) . 2 j =1
(29)
Although Eq. (29) is obtained when i ≥ 2 , it can also be applied to the position calculation error of q1.
4. Step size adjustment strategy In this section, a strategy which enables the step size adapt to the geometry of surface is introduced. By referring to the mean value theorem, the local truncation error of qi can be rewritten as: 1 2
1 2
ε 3 = ε1 + ε 2 = ( u ''(δ i −1 )Δs 2 , v ''(ϕi −1 )Δs 2 ) ,
(30)
where both of δ i −1 and ϕi −1 are located between qi −1 and qi on curve C . Let λ denote the upper limit of the local truncation error of qi , we get: (|u ''(δ i −1 )|+|v ''(ϕi −1 )|) Δs 2 ≤ 2λ .
(31)
In each calculation step, we choose Δs =
2λ as the step size. (|u ''(δ i −1 )|+|v ''(ϕi −1 )|)
Then the position calculation error of qi is estimated as: 1 i −1 1 i −1 |Ε ip | = |¦ u ''(δ j )|Δs 2 + |¦ v ''(ϕ j )|Δs 2 ≤ iλ . 2 j =0 2 j =0
(32)
In the subsequent part, the approximate values of u ''( s) and v ''( s) are considered. From the derivation formula of composite function, it follows that: ªu '( s ) º r '( s ) = [ ru rv ] « ». ¬ v '( s ) ¼
(33)
Take the derivative of Eq. (33) on both sides with respect to s, and it follows: ªu ''( s ) º [ ru rv ] « v ''( s) » = M , ¬ ¼
(34)
where M = r ''( s ) − ruu u '( s ) 2 − 2ruv u '( s )v '( s ) − rvv v '( s ) 2 . By left multiply [ ru
rv ]
T
on both sides of Eq. (34), we get: ª| ru |2 « ¬« ru .rv
ru .rv º ªu ''( s ) º ª ru .M º »« ». »=« | rv |2 ¼» ¬ v ''( s ) ¼ ¬ rv .M ¼
(35)
By solving this system of linear equations, we can get u ''( s ) and v ''( s ) . Here is the detail of the step size adjustment strategy: We start with a starting point q0 (u0 , v0 ) and a starting direction T0. Further we have an initial step size h0 and λ . 1. By using the starting direction T0 and initial step size h0, q1 is obtained. ( n .T )n 2. u ''(q0 ) and v ''(q0 ) are obtained by solving Eq. (35), where − 1 0 0 h0
is used to
approach r ''(q0 ) . 3. Change the step size into h1 =
2λ , where h1 is the step size to calculate (| u ''(q0 ) | + | v ''( q0 ) |)
the position of q2. 4. The position of q2 is acquired as (u1 + u '(q1 )h1 , v1 + v '(q1 )h1 ) . Similarly, the subsequent geodesic points (q3, q4…) can be computed. In each step, we need to compute the position of qi and the step size hi for the next step. Although the time costs increases in each step, the overall time costs decreases generally. Remark : The quality of result is influence by step size, curve curvature and parameterization. Therefore, for different parameterization or curve curvature, we adjust the step size to maintain result with good quality.
5. Validation of the algorithm In this section, we perform comparisons between our method and some existing elegant methods. At first, our approach is compared with Beck’s approach on accuracy and efficiency as
both approaches work directly on the original surface. Then, both of the two methods mentioned above are compared with Ravi Kumar’s discrete method. Ravi Kumar’s discrete method had been extended to estimate geodesics on regular surfaces in Chen’s work. This method does not perform on regular surfaces directly. The discrete geodesic is computed at first and then projected to the regular surface. In this work, Liu’s point projection method is adopted [22]. Next, various surfaces are adopted to validate the proposed algorithms. All computer programs of the proposed algorithms are written in MATLAB 7.0 and run on a personal computer [Intel Core (TM) i5 Duo Processor, 2.8GHz, 2GB Memory] to test the accuracy and the efficiency. In the present work, all the geodesics are computed until they reach the boundary of the surface and the computation time is recorded. In order to measure the accuracy of the three kinds of algorithms, we record the distance dist (qco , qid ) between last point qco (u (qco ), v(qco )) of the computed geodesic and last point qid (u (qid ), v(qid )) of the ideal geodesic, where dist (qco , qid ) equals to | u (qco ) − u (qid ) | + | v(qco ) − v( qid ) | . The ideal geodesic on B-spline surface is not readily available, so we use the geodesic computed by using the fourth order Runge-Kutta method with extremely small step size as the ideal geodesic. Now, we will introduce the first two examples, respectively. In both of the two examples, our approach is compared with Beck’s approach on accuracy and efficiency. For the two approaches , we record their computation errors at a group of equal step size and their computation time required to achieve a given level of accuracy, respectively. For the purpose of equity, step size of the two algorithms in both of the examples is fixed. Meanwhile, the accuracy and efficiency of Ravi Kumar’s discrete method are also presented. In the first example, we consider a 2 × 2 B-spline surface. As illustrated in Fig. 4, the geodesics are computed on this surface in three different directions. The results of the comparison are demonstrated in Fig. 5-7. The results of Ravi Kumar’s discrete method are shown in Fig. 8. In this work, we simply record the efficiency of Ravi Kumar’s discrete method in V0(du:dv=1:1) as the results in V1(du:dv=1:2) and V2(du:dv=2:1) share a quite similar tendency as V0.
Figure. 4. Geodesics computed on a 2 × 2 B-spline surface and its equivalent tessellated surface in three different directions
0.1
1E-3 1E-4
Ours Beck's
1E-5
1E-4
1E-3
0.01
Computation time (s)
Error
0.01
1E-6
Ours Beck's
1000
100
10
1
0.1
0.1
1E-6
1E-5
1E-4
Step Size
1E-3
0.01
0.1
Error
(a) Error versus step size (b) Computation time versus error Figure. 5. Results of comparison between our approach and Beck’s approach in V0 for Example 1.
1E-4
1E-5
Ours Beck's
Computation time (s)
1E-3
Accuracy
Our Beck's
1000
0.01
1E-3
0.01
10
1
0.1 1E-7
1E-6 1E-4
100
0.1
1E-6
1E-5
Step size
1E-4
1E-3
0.01
0.1
Error
(a) Error versus step size (b) Computation time versus error Figure. 6. Results of comparison between our approach and Beck’s approach in V1 for Example 1.
Our Beck's
1000
1E-3
1E-4
Our Beck's 1E-5
Computation time (s)
Accuracy
0.01
100
10
1
0.1 1E-4
1E-3
0.01
Step size
0.1
1E-6
1E-5
1E-4
1E-3
0.01
0.1
Error
(a) Error versus step size (b) Computation time versus Error Figure. 7. Results of comparison between our approach and Beck’s approach in V2 for Example 1.
du:dv=1:1 du:dv=1:2 du:dv=2:1
10
Time (s)
Error
0.01
1E-3
R P C
100
1 0.1 0.01
1E-4
0
5000
10000
15000
20000
25000
0
5000
Grid quantity
10000
15000
20000
25000
Grid quantity
(a) Error versus grid quantity (b) Efficiency versus grid quantity in V0 in three directions Figure. 8. Results of Ravi Kumar’s discrete approach for Example 1. Grid quantity: the number of triangular meshes on tessellated surface. R: Reorganization time. The time needed to reorganize data structure of tessellated surface. It does not include the time needed to generate the tessellated surface. P: Projection time. The time needed to project discrete geodesic to original surface. C: Computation time. The time needed to compute discrete geodesic on tessellated surface.
For the 2 × 2 B-spline surface, as can be seen in Fig. 5(b), Fig. 6(b) and Fig. 7(b), although Beck’s approach is more accurate, our method are always able to achieve a given level of accuracy with less computation time, which indicates that our method is more efficient. This is mainly due to that the divergence of accuracy between the two approaches is not significant on the 2 × 2 B-spline surface. In the second example, we consider a 6 × 4 B-spline surface. The geodesics are computed on this surface in three different directions, as illustrated in Fig. 9. The ideal geodesic on B-spline surface is calculated in a same way as the first example. The results of the comparison are presented in Fig. 10-12. The results of Ravi Kumar’s discrete method are shown in Fig. 13.
Figure. 9. Geodesics computed on a 6 × 4 B-spline surface and its equivalent tessellated surface in three different directions
0.01
1E-4
1E-5
Our Beck's
Computation time (s)
1E-3
Error
Our Beck's
1000
1E-6
100
10
1 1E-3
0.01
0.1
1E-6
1
1E-5
Step size
1E-4
1E-3
Error
(a) Error versus step size (b) Computation time versus error Figure. 10. Results of comparison between our approach and Beck’s approach in V0 for Example 2.
1000
Error
1E-3 1E-4 1E-5
Our Beck's
1E-6 1E-3
0.01
0.1
Computation time (s)
0.01
Ours Beck's
100
10
1
1
1E-6
1E-5
Step size
1E-4
1E-3
0.01
Error
(a) Error versus step size (b) Computation time versus error Figure. 11. Results of comparison between our approach and Beck’s approach in V1 for Example 2.
0.01
1E-4
1E-5
Our Beck's 1E-3
0.01
0.1
Step length
1
Computation time (s)
Error
1E-3
1E-6
Our Beck's
1000
100
10
1
0.1
1E-6
1E-5
1E-4
1E-3
0.01
Accuracy
(a) Error versus step size (b) Computation time versus error Figure. 12. Results of comparison between our approach and Beck’s approach in V2 for Example 2.
du:dv=1:1 du:dv=1:2 du:dv=2:1
10
Error
Time (s)
1E-3
R P C
100
1 0.1 0.01
1E-4
0
5000
10000
15000
Grid quantity
20000
25000
0
5000
10000
15000
20000
25000
Grid quantity
(a) Error versus grid quantity (b) Efficiency versus grid quantity in V0 in three directions Figure. 13. Results of Ravi Kumar’s discrete approach for Example 2. The difference on efficiency of the two kinds of algorithm for the 6 × 4 B-spline surface is clearly illustrated in Fig. 10(b), Fig. 11(b) and Fig. 12(b). As shown in these diagrams, the comparison results of the two algorithms on efficiency in V0 and V1 depend on the level of the specified accuracy. To be more specifically, our approach is more efficient when the specified computation error is greater than 5 × 10-5 while Beck’s approach is more efficient when the computation error is less than 1× 10-5 . This is because Beck’s method is much more accurate than our approach with the equal step size when the specified computation error is less than 1× 10-5 . In this case, to obtain equivalent accuracy as Beck’s method, our approach generally have to adopt very small step size, which is at the expense of computational efficiency. The comparison results in V2 demonstrate that the divergence of accuracy between the two kinds of approaches is not significant and our method is more efficient, which is quite different from the results in V0 and V1. This is mainly due to the improvement of accuracy of our approach in V2, which can be seen when we compare Fig. 12(a) with Fig. 10(a) and Fig. 11(a). As shown Fig. 5-7 and Fig. 10-12, when the step size decreases, the errors generally decrease. Meanwhile, the time costs increase. For each step size, the errors of Beck’s algorithm are always smaller than those of our method, and this is more obvious in example 2, which indicates that Beck’s method is more accurate. However, as for the efficiency of the two kinds of algorithms about achieving a given level of accuracy, the results are different. Since our approach is less accurate, to acquire equivalent accuracy as Beck’s method, our approach have to adopt smaller step size. The step size of our approach depends on the level of divergence of accuracy between the two approaches. If the accuracy of Beck’s method is significantly better than our approach, the latter method might have to use very small step size, which will increase the time cost largely. It is well known that Beck’s method is computationally intensive. In this work, we have found that the time costs of Beck’s method is approximately 7 times larger than that of our approach with the same step size. Then, we have an estimation that our algorithm will be more efficient when the step size of our approach is bigger than 1/7 of the step size of Beck’s approach. The computational results of Ravi Kumar’s discrete method are presented in Fig. 8 and Fig. 13. As shown in these figures, the accuracy and time costs of Ravi Kumar’s discrete method are influenced by the grid quantity instead of the step size. With the increase of the grid quantity, the tessellated surface gets closer to the original surface, and the computation errors generally
decrease. However, the reorganization time increase obviously. The results show that the time needed to compute geodesic on tessellated surface is much less than that needed to reorganize data structure of the tessellated surface. The reorganization process is computationally expensive, so Ravi Kumar’s discrete method is more suitable for computation of massive geodesics on one surface. In this case, the average reorganization time in each computation is reduced largely. As shown in Fig. 5-8 and Fig. 10-13 , all of the three algorithms are able to compute geodesic on surface. Both of our algorithm and Beck’s method have better accuracy than Ravi Kumar’s discrete method as the former two directly use information of the original surface. Experimental results show that our algorithm is efficient, and the geodesics which are completely on the surfaces can be obtained. Meanwhile, we do not need to solve the differential equations of geodesic which are quite complicated. Next, in order to validate the effectiveness of our step size adjustment strategy, we will introduce another two examples. For convenience, our approach with step size adjustment strategy is called variable step algorithm, and with fixed step size is called fixed step algorithm. Then, variable step algorithm is compared with fixed step algorithm on accuracy and efficiency. The step size of the fixed step algorithm is determined by the minimum step size of the variable step algorithm. In the third example, the geodesics are computed on the surface given by the first example in three different directions. The results of the comparison are recorded in Table 1. Table 1 Results of comparison for Example 3. q0=(0, 0), V0(du:dv=1:1), λ = 1e − 5 Type
Fixed step algorithm
Variable step algorithm
Step size ( Δs )
0.0104 1.41557e-03 1.417 1242
max=0.0538; min=0.0104 2.07882e-03 0.875 419
Type
Fixed step algorithm
Variable step algorithm
Step size ( Δs )
0.0106 1.15169e-03 1.306 1211
max= 0.0618; min=0.0106 4.27124e-03 0.797 390
Type
Fixed step algorithm
Variable step algorithm
Step size ( Δs )
0.0112 2.36930e-03 1.391 1199
max= 0.0486; min=0.0112 5.17535e-03 0.953 467
Errors Time (s) Number of discrete points q0=(0, 0), V1(du:dv=1:2), λ = 1e − 5
Errors Time (s) Number of discrete points q0=(0, 0), V2(du:dv=2:1), λ = 1e − 5
Errors Time (s) Number of discrete points
As shown in Table 1, the accuracy of the fixed step algorithm is better than the variable step algorithm. This is due to the fact that the step size of the former algorithm is always not greater than the latter. However, as for the computational efficiency, the variable step algorithm has superiority. Furthermore, the number of discrete points on the geodesics computed by the variable step algorithm is smaller than half of that of the fixed step algorithm, which means the geodesics
computed by the former algorithm are more compact. In the fourth example, we consider another 2 × 2 B-spline. The geodesics are computed on this surface in three different directions, as illustrated in Fig. 14. The results of the comparison are recorded in Table 2.
Figure. 14. Geodesics computed on a 2 × 2 B-spline surface in three different directions Table 2 Results of comparison for Example 4. q0=(0, 0), V0(du:dv=1:1), λ = 1e − 5 Type
Fixed step algorithm
Variable step algorithm
Step size ( Δs )
0.0294 3.85031e-03 0.641 507
max=0.0800; min=0.0294 5.66168e-03 0.765 366
Type
Fixed step algorithm
Variable step algorithm
Step size ( Δs )
0.0162 2.47817e-03 0.828 690
max= 0.0843; min=0.0162 2.97750e-03 0.623 297
Type
Fixed step algorithm
Variable step algorithm
Step size ( Δs )
0.0196 2.27483e-03 0.797 677
max= 0.133; min=0.0196 3.39647e-03 0.625 297
Errors Time (s) Number of discrete points q0=(0, 0), V1(du:dv=1:2), λ = 1e − 5
Errors Time (s) Number of discrete points q0=(0, 0), V2(du:dv=2:1), λ = 1e − 5
Errors Time (s) Number of discrete points
As shown in Table 2, the comparison results of Example 4 in direction V1 and direction V2 are quite similar to the results of Example 3. However, as for direction V0, it takes the variable step algorithm more time on computational process even though the number of steps decreases, which means the reduced time costs due to the decrease of step numbers does not make up for the sum of the increased time costs in each step. As illustrated in Table 1 and Table 2, the average time costs in each step of the fixed step algorithm and the variable step algorithm are 1.2ms and 2.1ms, respectively. According to this, we have an estimate that the variable step algorithm is more efficient than the fixed step algorithm when the number of discrete points computed by the former
algorithm is smaller than 4/7 of that of the latter. Generally speaking, our variable step algorithm works well for most of the cases. The location of discrete points obtained by this algorithm distributes reasonably according to the geometry of the surface, and this could contribute to outstanding reduction of the number of discrete points on the geodesic, which means the geodesic computed by the variable step algorithm is more compact.
6. Conclusion This paper presents a geometric method for computation of geodesic on parametric surfaces. The geodesic which is completely on the surface can be obtained by this algorithm. Furthermore, we do not need to solve the differential equations of geodesic which are quite complicated and generally not easy to solve. Experimental results indicate that our approach is more efficient than Beck’s method on a 2 × 2 B-spline surface while the efficiency of our approach on a 6 × 4 B-spline surface depends on the tracing direction and the level of the specified accuracy. Generally speaking, our approach is more suitable in two cases. First, the parametric surface is relative flat. Second, the requirement for the accuracy is modest. In these cases, the accuracy of our algorithm is guaranteed. At the same time, our approach will trace the geodesic efficiently as we do not need to use quite small step size. One possible application of our method is path planning in the tape laying process. In this process, the composite tapes generally follow geodesic path on a mould surface which is nearly flat or nearly conventional cylinderical. Comparison results show that both of our algorithm and Beck’s method have better accuracy than Ravi Kumar’s discrete method. At the same time, Ravi Kumar’s discrete method is more suitable for computation of massive geodesics on one surface. The results also demonstrate that our step size adjustment strategy works well for most of the cases. The step size adjustment strategy has contributed to outstanding reduction of the number of discrete points on the geodesic, which means the geodesic computed with the step size adjustment strategy is more compact. However, there are also limitations in our algorithm. The function of the parameterized surface must be given in advance, which means this algorithm does not apply to general surfaces. In addition, the surface should contain no singular point. This is due to that we need to pull the direction of the curve back to the parametric domain in each step.
Acknowledgements This work was supported by the National Basic Research Program of China (2013CB035805) and the National Science and Technology Major Project of China (2010ZX04016-013).
References [1] Debout, P., Chanal, H., Duc, E., 2011. Tool path smoothing of a redundant machine: application to automated fiber placement. Computer-Aided Design 43 (2), 122–132. [2]
Azariadis, P., Aspragathos, N.A., 1997. Design of plane developments of doubly curved surfaces. Computer-Aided Design. 29 (10), 675–685.
[3]
Azariadis, P., Aspragathos, N.A., 2001. Geodesic curvature preservation in surface flattening through
constrained global optimization. Computer-Aided Design 33, 581–591. [4] Peyré, G. Péchaud, M. Keriven, R. Cohen, L., 2010. Geodesic Methods in Computer Vision and Graphics. In Found. Trends. Comput. Graph. Vis, 197-397. [5] Do Carmo, M.P., 1976. Differential Geometry of Curves and Surfaces. Prentice-Hall, Englewood Cliffs, NJ. [6] Beck, J.M., Farouki, R.T., Hinds, J.K., 1986. Surface analysis methods. IEEE Computer Graphics and Applications 6, 18–36. [7] Patrikalakis, N.W., Badris, L., 1989. Offsets of curves on rational B-spline surfaces. Engineering with computers 5, 39–46. [8] Sneyd, J., Peskin, C.S., 1990. Computation of geodesic trajectories on tubular surfaces. SIAM Journal of Scientific Statistical Computing 11, 230–241 [9] Hotz, I., Hagen, H., 2000. Visualizing geodesics. In: Proceedings IEEE Visualization, Salt Lake City, UT, PP. 311–318 [10] Maekawa, T., 1996. Computation of shortest paths on free-from parametric surfaces. Journal of Mechanical Design, ASME Transactions 118 (4), 499–508. [11]
Ying, L.X., Candes, E.j., 2006. Fast geodesics computation with the phase flow method. Journal of Computational Physicsm 220 (1), 6-18.
[12] Kasap, E., Yapici, M., Akyildiz, F.T., 2005. A numerical study for computation of geodesic curves. Applied Mathematics and Computation 171 (2), 1206–1213. [13] Chen, S.-G., 2010. Geodesic-like curves on parametric surfaces. Computer Aided Geometric Design 27, 106–117. [14] Tucker, C.L., 1997. Forming of advanced composites. In: Gutowski, T.G. (Ed.), Advanced Composites Manufacturing. Wiley, New York. [15] Ravi Kumar, G.V.V., Prabha Srinivasan., Devaraja Holla, V., Shastry, K.G., Prakash, B.G., 2003. Geodesic curve computations on surfaces. Computer Aided Geometric Design 20, 119–133. [16] Polthier, K., Schmies, M., 1998. In: Hege, H.C., Polthier, H.K. (EDS.), Straightest Geodesics On Polyhedral Surfaces in Mathematical Visualization. Springer-Verlag, Berlin. [17]
Kanai, T., Suzuki, H., 2001. Approximate shortest path on a polyhedral surface and its applications, Computer Aided Design 33 (11), 801–811.
[18] Martinez, D., Velho, L., Carvalho, P.C., 2005. Computing geodesics on triangular meshes. Computer & Graphics 29, 667–675. [19] Surazhsky, V., Surazhsky, T., Kirsanov, D., Gortler, S., Hoppe, H., 2005. Fast exact and approximate geodesics on meshes. In: Proc. of SIGGRAPH 2005. ACM Transactions on Graphics 24 (3), 553–560. [20]
Bose, P., Maheshwari, A., Shu, C., Wuhrer, S., 2011. A survey of geodesic paths on 3D surfaces. Computational Geometry: Theory and Applications 44 (9), 486–498.
[21] Piegl, L., Tiller, W., 1997. The NURBS book, 2nd Edition. Springer, New York. [22] Liu, X,-M., Yang, L., Yong, J.-H., Gu, H,-J., Sun, J,-G., 2009. A torus patch approximation approach for point projection on surfaces. Computer Aided Geometric Design 26 (5), 593–8.