- Email: [email protected]
Geodesic-like curves on parametric surfaces
Computer Aided Geometric Design 27 (2010) 106–117
Contents lists available at ScienceDirect
Computer Aided Geometric Design www.elsevier.com/locate/cagd
Geodesic-like curves on parametric surfaces Sheng-Gwo Chen 1 Department of Applied Mathematics, National Chiayi University, Chia-Yi 600, Taiwan
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 5 March 2009 Received in revised form 31 August 2009 Accepted 1 October 2009 Available online 3 October 2009
Finding the geodesic on a surface is an important task in many fields. As the traditional discrete algorithms or numerical approach can only find a list of discrete points, this study aims to propose a new, elegant and accurate method for approximating geodesics on a surface. The proposed geodesic-like curve is a mathematical curve that approaches to the geodesic on surface when the order of the curve reaches to infinity. The geodesic-like curve is also compared to some other well-known methods in this study. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Geodesics play an important key of differential geometry and a powerful tool to improve many problems in CAGD, CAD/CAM, computational geometry, computer vision, etc. (Paluszny, 2008; Ravi Kumar et al., 2003a; Sánchez-Reyes and Dorado, 2008; Sprynski et al., 2008). A curve α on a regular surface is a geodesic if its osculating plane at each point on the curve contains the normal vector of the surface at that point. The osculating plane at α (t 0 ) is a plane generated by the tangent vector and the normal vector of α at α (t 0 ). A geodesic on surface can be also defined by a curve with zero geodesic curvature at each point or a solution of the system of geodesic equations (2). In the calculus of variation, the geodesic curve is a critical point of the length function or the energy function in Eqs. (6) and (7). Hence, the geodesics joining two points on surfaces have extreme (maximal or minimal) lengths and extreme energies. More details of geodesics can be found in Do Carmo (1976). Analytical approaches are quite complex and difficult to find on surfaces (Do Carmo, 1976). In 2000, Hotz and Hagen (2000) improved the geodesic from the first, second fundamental forms and numerical method. The numerical approach presented by Emin Kasap et al. (2005) for the computation of geodesics between two fixed points on a surface is elegant through the use of the system of geodesic equations in Eq. (2) and finite-difference. Many accurate discrete methods approach geodesics or the shortest paths on tessellated surfaces (Tucker, 1997), polygonal surfaces (Polthier and Schmies, 1998) and triangular meshes (Martinez et al., 2005; Surazhsky et al., 2005). These methods are listed in the references. A new method for estimating the geodesic is introduced in this paper. The geodesic-like curve which is a linear combination of base functions is a critical point of our energy function in Eq. (15). The base functions are considered as the Bernstein polynomials or B-spline functions in this paper. Obviously, the geodesic-like curve approaches the geodesic on surface when its order is large enough. The definition of geodesic-like curves on surfaces and some properties of geodesiclike curves are provided in Section 2. In Section 3, conditions to improve different kinds of geodesic problems are discussed and the method of geodesic-like curves with other recent researches.
1
E-mail address: [email protected]. Partially supported by NSC, Taiwan.
0167-8396/$ – see front matter doi:10.1016/j.cagd.2009.10.001
© 2009 Elsevier B.V.
All rights reserved.
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
Fig. 1. A proper variation of
107
α.
2. Geodesics and geodesic-like curves A geodesic on a Riemannian manifold ( M , g ) is a smooth curve
D ds
α = 0,
α (s) : [a, b] → M which satisfies the differential equation (1)
D d where ds is the covariant derivative associated to g and α = ds α (s). In a local coordinate system (x1 , . . . , xn ), geodesics can be represented by the solutions of the system of differential equations
d 2 xi dt 2
+
dxk dx j Γ jki x(t ) = 0 for each i ∈ {1, 2, . . . , n}. dt dt
j ,k
(2)
The system (2) is known as the differential equations of the geodesics of M. Now, we shall discuss geodesics on surfaces from the variation of calculus. Let S be a regular surface and let b be a parametrization on S, b : U ⊂ R2 → S. Hence, b is a diffeomorphism between U and b(U ). If α is a smooth curve on S with a parameter α : [a, b] → U , then the arc length of α is defined by
b d L (α ) = b α (s) ds ds
(3)
a
and the energy of
α is
b 2 1 d b α (s) ds. E (α ) = 2 ds
(4)
a
α , see Fig. 1, is a differentiable map h : [a, b] × [−ε, ε] → U such that h(s, 0) = α (s), s ∈ [a, b], (5) h(a, t ) = α (a), t ∈ [−ε , ε ], h(b, t ) = α (b), t ∈ [−ε , ε ]. Intuitively, a proper variation of α , ht = h(·, t ), is a family of curves which are differentiable on the parameter t ∈ (−ε , ε ) ∂(b◦h) ∂(b◦h) with the same initial point α (a) and the same end point α (b) and such that h0 agrees with α . Thus ∂ s and ∂ t are differentiable vector fields along α . The energy of the curve ht is represented as 2 b 1 ∂(b ◦ h) E (t ) = (s, t ) (6) ds, t ∈ [−ε , ε ], 2 ∂s A proper variation of
a
and the length of ht is represented as
L (t ) =
b ∂(b ◦ h) ds, ( s , t ) ∂s a
t ∈ [−ε , ε ].
(7)
108
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
From the first variation formula of E (t ),
E (0) =
b
D
V (t ),
dt
d
b◦α
dt
dt ,
(8)
a
where V (t ) is the variational field of b ◦ h, Theorem 1 is obvious. Theorem 1. Suppose that b : U → S is an parametrization on S and α is a curve parametrized by arc length on surface S. E (t ) and L (t ) are defined by Eqs. (6) and (7), respectively. The following statements are equivalents. (1) b(α (s)) is a geodesic on surface S. (2) α (s) has extreme energy on S, that is E (0) = 0. (3) α (s) has extreme length on S, i.e. L (0) = 0. Hence, the shortest path from p to q is a geodesic. We can find these properties in Do Carmo (1976, 1992). For simplicity,
α : [a, b] → U is called a geodesic on S if b(α ) is a geodesic on S in this study. According to the Weierstrass theorem (Marsden and Hoffman, 1993), any continuous real function can be approximated by the Bernstein polynomial. We summarize the Weierstrass theorem as follows. Theorem 2 (Weierstrass theorem). Let f : [0, 1] → R be a continuous function and let > 0. Then the sequence of Bernstein polynomials
p n (s) =
n
i
B ni (s) f
(9)
n
i =0
converges uniformly to f as n → ∞, where B ni (s) = C in (1 − s)n−i si is the Bernstein polynomial. Thus geodesics on surfaces can be approximated by the Bézier curves. Theorem 3. Let S be a parametric surface with a parametrization b : U → S and let K be a compact subset in U . Suppose that γ : [a, b] → K is a geodesic on S with no any conjugate point. If the curve c˜ n (s) = ni=0 B ni (s)˜c i in K is a solution of the optimal problem:
2
n b ∂ n min E (u 1 , . . . , un−1 , v 1 , . . . , v n−1 ) = B i (s)(u i , v i ) ds b 2 ∂s 1
s.t.
i =0
a
(u 0 , v 0 ) = γ (a), (un , v n ) = γ (b), (u i , v i ) ∈ K for each i ,
(10)
then limn→∞ c˜ n = γ . Proof. It is clear that {˜cn } is a uniform convergent sequence. Without the loss of generality, we assume that [a, b] = [0, 1]. By Theorem 2, a sequence of Bézier curves exists n
c (s) =
n
i
B ni (s)
γ
(11)
n
i =0
γ as n → ∞. As c˜ n is a solution of Eq. (10), we have n E (γ ) E c˜ E cn
converges uniformly to
(12)
for each n. In addition, {˜c } and {c } are uniform convergent sequences. Thus n
n
E (γ ) lim E c˜ n = E lim c˜ n lim E cn = E lim cn = E (γ ). n→∞
n→∞
n→∞
n→∞
(13)
This implies that E (limn→∞ c˜ n ) = E (γ ), that is limn→∞ c˜ n has minimal energy. Therefore, limn→∞ c˜ n is a geodesic between γ (a) and γ (b) on S. Because γ is the unique minimum geodesic between γ (a) and γ (b), limn→∞ c˜ n agrees with γ . The curve c˜ n in the above theorem is called a geodesic-like curve. The definition of the geodesic-like curve is stated as follows.
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
109
Definition 4. Let b(u , v ) be a parametrization of a parametric surface S, b : U ⊂ R2 → S and let { f i (s)| f i : [a, b] → R}ni=0 be a n sequence of real functions. A curve c˜ (s) on U is called a geodesic-like curve of order (n + 1) on S if c˜ (s) = i =0 f i (s)(u˜ i , v˜ i ) and satisfies
(∇ E )(u˜ i , v˜ j ) = 0,
(14)
where
2
n b ∂ E= b f i (s)(u i , v i ) ds 2 ∂s 1
(15)
i =0
a
is the energy of the curve c (s) =
n
i =0
f i (s)(u i , v i ).
The function f i (s) is called the base function of geodesic-like curve and c (s) is the base curve of geodesic-like curve. We know that any continuous function can be approximated by the Bernstein polynomials (see Theorem 2) or the B-spline functions. If f i (s) is a Bernstein polynomial or B-spline function, the geodesic-like curve approaches a geodesic on S when its order n is large enough. The base function f i is assumed a Bernstein polynomial or B-spline function in this paper. From Theorem 3, the following theorem is obvious. Theorem 5. The solutions of the following optimization problem are geodesic-like curves:
1 d 2 min E (u i , v j ) = b c (s) ds 2 ds 1
0
⎧ ⎨ c (s) = ni=0 f i (s)(u i , v i ) ∈ U , s.t. c (a) = x0 , ⎩ c (b) = x1 .
(16)
Moreover, the geodesic-like curve approaches the shortest path between x0 and x1 . Its length approaches the distance between x0 and x1 when n is large enough. 2 3 Let us present the gradient and the hessian of the nenergy function E in Eq. (15). Denote b : U ⊂ R → R is a parametrization of a parametric surface S and c (s) = i =0 f i (s)(u i , v i ) is a base curve on U . For simplicity, b denotes b(c (s)), c denotes c (s) and the notation means the derivative of function with respect to the parameter s. The energy of c is the function
E (u i , v j ) =
1
b
2
b , b ds.
a
Set c (s) = (u (s), v (s)) where u (s) =
(17)
n
i =0
∂ u (s) = f i (s), ∂ ui ∂ u ( s ) = 0, ∂ vi ∂ v ( s ) = 0, ∂ ui ∂ v (s) = f i (s), ∂ vi
f i (s)u i and v (s) =
n
i =0
f i (s) v i . Since
(18)
we obtain
∂ ∂ ∂ b = b u (s), v (s) ∂ ui ∂ ui ∂ s ∂ (bu u + b v v ) = ∂ ui = buu u f i + bu f i + buv v f i = (buu u + buv v ) f i + bu f i
(19)
110
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
and
∂ b = (buv u + b v v v ) f i + b v f i . ∂ vi Hence,
E ui = E vi =
b
(buu u a b (buv u a
+ buv v ) f i + bu f i , bu u + b v v ds, + b v v v ) f i + b v f i , bu u + b v v ds.
(20)
(21)
The gradient of E is
∇ E = ( E u 0 , E u 1 , . . . , E un , E v 0 , E v 1 , . . . , E v n ).
(22)
Therefore, c (s) is a geodesic-like curve if it is a solution of
b
(buu u + buv v ) f i + bu f i , bu u + b v v ds = 0,
a
b
(buv u + b v v v ) f i + b v f i , bu u + b v v ds = 0,
(23)
a
for each i ∈ {0, 1, . . . , n}. The system of equations in (15) or (23) is called the system of geodesic-like equations. Since
b E ui = a
b Ev j = a
∂ b (t ), b (t ) dt , ∂ ui ∂ ∂v j
b (t ), b (t ) dt ,
(24)
we have
⎧ b ∂ 2 ∂ ∂ ⎪ ⎪ E u i u j = a ∂ u i ∂ u j b (t ), b (t ) + ∂ u i b (t ), ∂ u j b (t ) dt , ⎪ ⎪ ⎪ b ∂ 2 ∂ ⎪ ∂ ⎨E ui v j = a ∂ u i ∂ v j b (t ), b (t ) + ∂ u i b (t ), ∂ v j b (t ) dt , ⎪ E vi u j = Eu j vi , ⎪ ⎪ ⎪ ⎪ b ∂ 2 ∂ ⎪ ⎩E ∂ vi v j = a ∂ v i ∂ v j b (t ), b (t ) + ∂ v i b (t ), ∂ v j b (t ) dt .
(25)
By Eqs. (19) and (20), we have
⎧ ∂ 2 b ⎪ ∂ u i ∂ u j = (buuu u + buuv v ) f i f j + buu ( f i f j + f i f j ) ⎪ ⎪ ⎪ ⎪ ⎪ = (buuu u + buuv v ) A i j + buu ( A i j ) , ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ∂ 2 b = (buuv u + buv v v ) f i f j + buv ( f i f + f f j ) j i ∂ ui ∂ v j ⎪ = (buuv u + buv v v ) A i j + buv ( A i j ) , ⎪ ⎪ ⎪ ⎪ ⎪ ∂ 2 b ⎪ ⎪ ∂ v i ∂ v j = (buv v u + b v v v v ) f i f j + b v v ( f i f j + f i f j ) ⎪ ⎪ ⎪ ⎩ = (buv v u + b v v v v ) A i j + b v v ( A i j ) ,
(26)
where A i j = f i f j . Thus the hessian of E is the matrix
E uu E vu
E uv Evv
,
(27)
where E uu = ( E u i u j ), E uv = E vu = ( E u i v j ) and E v v = ( E v i v j ). The constant curve is a trivial solution of the system of geodesic-like equations in Eq. (23), because the nonzero vector ω = (0, u 1 − u 0 , u 2 − u 0 , . . . , un − u 0 , v 1 − v 0 , . . . , v n − v 0 ) is always in a decreasing direction. In next section, different kinds of initial conditions of the system of geodesic-like equations are discussed and the proposed method is compared with others.
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
111
3. Some elegant methods and geodesic-like curves In this section, some elegant and accurate methods of solving geodesic problems are presented and compared with the proposed method. 3.1. Kasap’s numerical method Finding the shortest path between two points on surfaces is a classical problem and it has many important applications (Paluszny, 2008; Sprynski et al., 2008; Sánchez-Reyes and Dorado, 2008). A numerical research presented by Kasap et al. (2005) for the computation of geodesic between two points on surfaces is elegant and accurate. They solve the system (2) via the finite-difference method and the iterative method. By the centered-difference formulas, the differential equations of the geodesics
⎧ 2 ⎨ d 2u + Γ 1 du 2 + Γ 1 du dv + Γ 1 dv 2 = 0, 11 ds 12 ds ds 22 ds ds ⎩ d2 v + Γ 2 du 2 + Γ 2 du dv + Γ 2 dv 2 = 0 11 ds 12 ds ds 22 ds ds2
(28)
can be regard as the system
U i +1 − 2U i + U i −1 = h2 f (·), V i +1 − 2V i + V i −1 = h2 g (·),
(29)
where
1 f (·) = −Γ11 (U i , V i )
U i +1 − U i −1
2
2h
U i +1 − U i −1 V i +1 − V i −1 1 − Γ12 (U i , V i )
2 V i +1 − V i −1 1 − Γ22 (U i , V i )
2h
2h
(30)
2h
and
2 g (·) = −Γ11 (U i , V i )
U i +1 − U i −1 2h
2
U i +1 − U i −1 V i +1 − V i −1 2 − Γ12 (U i , V i )
2 V i +1 − V i −1 2 − Γ22 (U i , V i ) .
2h
2h
2h
(31)
The system (29) is a system of nonlinear equations with parameters {U i , V j }. It can be approximated by the Newton’s method, the iterative method or other numerical methods. In Figs. 2 and 3, the yellow curves are the Kasap’s geodesics between two points on the Witeny’s umbrella and a spherical surface. Let b be a parametrization on surface S and let p , q be two distinct points on S. There exists two points (u˜ 0 , v˜ 0 ), (u˜ n , v˜ n ) in U such that
b(u˜ 0 , v˜ 0 ) = p , b(u˜ n , v˜ n ) = q.
(32)
We consider the base curve c : [a, b] → U in the Definition 4 satisfies c (a) = (u˜ 0 , v˜ 0 ) and c (b) = (u˜ n , v˜ n ), i.e. the starting point and the terminal point of c are fixed. As E u 0 , E v 0 , E un and E v n are vanished, the system of geodesic-like equations in Eq. (15) becomes
∇ E ( u i , v j ) = ( E u 1 , E u 2 , . . . , E u n −1 , E v 1 , E v 2 , . . . , E v n −1 ) = 0.
(33)
The solutions of the system (33) are the geodesics between p and q. Furthermore, if the hessian of E is a positive-definite matrix, then this geodesic-like curve approaches the shortest path between p and q as n → ∞. The system (33) is also a system of nonlinear equations. This system can certainly improved by the Newton’s method, the fixed point theorem or other numerical methods. The Kasap’s geodesics and the geodesic-like curves on the Witeny umbrella and the spherical surface obtained are show in Figs. 2 and 3. The Kasap’s geodesics are yellow curves and the geodesic-like curves are red curves. In Fig. 3, the black curve is an exact geodesic on the spherical surface. In these simulations, the Kasap’s geodesic has 100 vertices and the geodesic-like curve comes from a cubic B-spline curve with 10 control points in U . In Fig. 2, the length of geodesic-like curve is less than the length of Kasap’s geodesic. The geodesic-like curve is closer to the exact geodesic than the Kasap’s geodesic on spherical surface and as shown in Fig. 3. Both of Kasap and our methods study the geodesic problem via a system of nonlinear equations. The number of control points of geodesic-like curve is less than the number of vertices in Kasap’s geodesic and the geodesic-like curve is more accurate than the Kasap’s geodesic.
112
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
Length Time (second)
Kasap’s geodesic (yellow curve)
Geodesic-like curve (red curve)
1.553 2.012
1.550 0.354
Fig. 2. The Kasap’s geodesic and geodesic-like curve on Witeny umbrella. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
Time (second)
Kasap’s geodesic
Geodesic-like curve
4.21
0.407
Fig. 3. The Kasap’s geodesic and geodesic-like curve on spherical surface. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
3.2. Ravi Kumar’s discrete method In this paragraph, we shall study the problem of starting point and its direction. It is well known that for each p ∈ S and y ∈ T p S, there exists an unique geodesic such that α (0) = p and α (0) = y locally (Do Carmo, 1976). Many researches improve this problem by a polygon and a triangular mesh rather than a regular surface. Martinez et al. (2005), Surazhsky et al. (2005) and Ravi Kumar et al. (2003b) present some different discrete methods for solving this problem. All of these methods are accurate. The Ravi Kumar’s method is described and compared with the proposed method. The main idea of Ravi Kumar’s method is the orthogonal projection. The method 2b in Ravi Kumar’s paper (Ravi Kumar et al., 2003b) is summarized as below. (1) Given a point p 0 and a tangent direction v 0 . (2) Add p 0 at the first point of discrete geodesic P. For each i 0 (3) Determine the point p i whether a vertex, in a edge or in a triangle and compute the normal vector ni at p i . (4) Let t i be the unit normalized projection of the vector v i onto some facet on the triangular mesh or tessellation and shown in Fig. 4. (5) Determine the next point p i +1 by the intersect of the line from p i with direction t i and the edges of triangle.
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
113
Fig. 4. The next geodesic direction t i .
Although Ravi Kumar’s method works on a tessellation or a triangular mesh, it can be extended to estimate geodesics on regular surfaces easily. It only needs to project each point of discrete geodesic to the regular surface via the normal vector field on surface S. We can also improve this problem by geodesic-like curve. Given a point p ∈ U and a vector y ∈ R2 such that p + y ∈ U . Consider the conditions
c (0) = p ∈ U , c (0) = y ∈ R2 .
(34) and c (0)
Suppose that c (s) is a base curve defined on [0, 1], we have c (0) = (u 0 , v 0 ) = n((u 1 , v 1 ) − (u 0 , v 0 )). Thus (u 0 , v 0 ) = p and (u 1 , v 1 ) = ny + p are fixed in U . The system of geodesic-like equations becomes
∇ E ( u i , v j ) = ( E u 2 , E u 3 , . . . , E u n , E v 2 , E v 3 , . . . , E v n ) = 0.
(35)
However, these equations have a trivial solution (u i , v i ) = (u 1 , v 1 ) for each i ∈ {2, 3, . . . , n}. To avoid this problem, we set the terminal point of c (s) fixed, the system of geodesic-like equations becomes
∇ E ( u i , v j ) = E ( u 2 , u 3 , . . . , u n −1 , v 2 , v 3 , . . . , v n −1 ) = 0.
(36)
If the set {u i }ni=2 or the set { v i }ni=2 is fixed, then the system of geodesic-like equations correspond to
∇ E (u i , v j ) = E ( v 2 , v 3 , . . . , v n ) = 0
(37)
∇ E ( u i , v j ) = E ( u 2 , u 3 , . . . , u n ) = 0.
(38)
or
Ravi Kumar’s geodesics and the geodesic-like curves on the Witeny umbrella and a torus are estimated to have 140 vertices and quadric B-spline curves with 20 control points respectively in our simulations. The Ravi Kumar’s geodesics are yellow curves and the geodesic-like curves are red curves in Figs. 5 and 6. In Figs. 5 and 6, the number of control points and the length of geodesic-like curve are fewer and shorter than the Ravi Kumar’s geodesic. Thus geodesic-like curve is more accurate that Ravi Kumar’s geodesic. 3.3. Geodesic between two curves The problem of the shortest path between two curves on surfaces is also important but most geodesic estimating methods can not improve this problem easily. We shall present the geodesic-like model for solving this problem. Let S be a parametric surface and let b : U → S be a parametrization on S. Assume that Γ1 and Γ2 are two curves on S with parametrization c 1 : [0, 1] → U and c 2 : [0, 1] → U , respectively. The mathematical model of the distance between Γ1 and Γ2 on S is
min E (c ) such that c : [a, b] → U , c ([a, b]) ∩ c 1 ([0, 1]) = c (a), c ([a, b]) ∩ c 2 ([0, 1]) = c (b).
(39)
114
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
Length Time (second)
Ravi Kumar’s geodesic (yellow curve)
Geodesic-like curve (red curve)
10.710 2.023
10.706 2.012
Fig. 5. The discrete geodesic and geodesic-like curve on Witeny umbrella surface. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
The length of exact geodesic (black curve) = 22.8828516
Ravi Kumar’s geodesic Geodesic-like curve
Length
Time (second)
22.9259979 22.8828616
2.047 2.051
Fig. 6. The discrete geodesic and geodesic-like curve on torus. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
Denote c 1 = c 1 (μ) = (u 0 (μ), v 0 (μ)) and c 2 = c 2 (ν ) = (un (ν ), v n (ν )). If the curve c is the base curve in Definition 4, then E is a real function with parameters μ, ν , u 1 , u 2 , . . . , un−1 and v 1 , v 2 , . . . , v n−1 . Hence, the solutions of (39) satisfy the system of geodesic-like equations
∇ E = ( E t 1 , E t 2 , E u 1 , E u 2 , . . . , E u n −1 , E v 1 , E v 2 , . . . , E v n −1 ) = 0.
(40)
where E μ and E ν can be computed by
E μ = E u0 E ν = E un
du 0 dμ dun dν
+ E v0 + E vn
dv 0 dμ dv n dν
, .
(41)
In other words, Eq. (40) is a system of geodesic-like equations with the conditions
c (0) ∈ C 1 , c (1) ∈ C 2 .
(42)
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
115
Fig. 7. The geodesic-like curve on closed surfaces (periodic surfaces). (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
Fig. 8. The minimal geodesic-like curve on a surface of revolution. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
Fig. 9. The closed geodesic on a surface of revolution. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
4. Conclusions and discussions Let S be a closed surface (periodic surface) with a parametrization b : U → S. As b is not an onto map, the geodesic-like curves should not defined on U . In this case, we consider the periodic map
b˜ : R2 → S
(43)
such that b = b˜ ◦ π where π is the projection map from R2 to U . And construct the geodesic-like curve by the energy function
E=
b ∂ 2 ˜ ◦ c (s) ds. b 2 ∂s 1
a
In this case, geodesic-likes curve should be defined on R2 . Geodesics on periodic surfaces are shown in Figs. 7–9.
(44)
116
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
Fig. 10. The geodesic-like curves with different number of control points. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
Number of control point Length of geodesic-like curve Time (second) Error Number of control point Length of geodesic-like curve Time (second) Error
4
5
6
7
8
31.1151331 0.04125 0.2318035 9 30.900999 0.323 0.0176694
30.9761609 0.0813 0.0928313 10 30.8892435 0.4404 0.0059139
30.9755519 0.13 0.0922223 20 30.8853687 2.64 0.0020391
30.9040219 0.18 0.0206923 40 30.8842772 6.39 0.0009476
30.9035892 0.234 0.0202596 60 30.8833296 20.12 <1.0e−7
Fig. 11. The length of geodesic-like curves by quadric B-spline curves.
The degree of geodesic-like curve Length of geodesic-like curve Time (second) Error
2
3
4
5
6
7
10
31.1151 0.015 0.2318
31.0047 0.062 0.12134
30.9756 0.219 0.09222
0.964 0.344 29.9193
30.9156 0.469 0.03231
30.9048 0.672 0.02151
30.8952 2.45 0.01184
Fig. 12. The length of geodesic-like curves by Bézier curves.
The geodesic-like curves in this section always come from uniform quadric B-spline curves and Bézier curves. The geodesic-like curves on periodic surfaces (closed surface) obtained and shown in Figs. 7, 8 and 9. The yellow curves on these figures are the boundary of surfaces and the red curves are the geodesic-like curves which come from the uniform quadric B-spline curves with 8 control points. In Fig. 7, the geodesic-like curves are not the shortest path on a torus and a surface of revolution. In Fig. 8, the geodesic-like curve agrees with the shortest path on surface. However, all of geodesiclike curves are close to geodesics on that surface. We also show a closed geodesic on a surface of revolution in Fig. 9. In Fig. 10, we show some geodesic-like curves with different orders. The order of red geodesic-like curve is 4 and that of yellow geodesic-like curve is 10. We list the results of the length of geodesic-like curves with different orders, the number of control points, in Fig. 11. Obviously, the geodesic-like curve approaches the minimal geodesic when the order is large enough. In Fig. 12, the geodesic-like curves are the Bézier curves. Similarly, the geodesic-like curve approaches the geodesic when the degree increasing. Because it is expensive to estimate the Bézier curve of large degree, the geodesic-like curve by the B-spline curve is more effective the Bézier curve. In these tests, the length of exact geodesic curve is 30.8833296. It is clear that the geodesic-likes curve also converges to the exact geodesic. From these figures, geodesic-like curves come from the uniform quadric B-spline curves are more accurate than the Bézier curves. The method of the geodesic-like curve is different from other traditional approaches. The advantages of geodesic-like curve are nature, accuracy and elegance. Furthermore, the geodesic-like curve is a mathematical curve. However, the accuracy of geodesic-like curve depends on different parametrizations on surface. It is well known that the critical points of length function in Eq. (7) agree with the critical points of energy function in Eq. (6) when the curves always have arc-length parameters (Do Carmo, 1976, 1992). Therefore, we can also define the geodesic-like curve via length function. Another definition is stated as follows. Definition 6 (Another definition of geodesic-like curve). Let b(u , v ) be a parametrization of a parametric surface S, b : U ⊂ R2 → S and let { f i (s) | f i : [a, b] → R}ni=0 be a sequence of real functions. A curve c˜ (s) on U is called a geodesic-like curve n of order (n + 1) on S if c˜ (s) = i =0 f i (s)(u˜ i , v˜ i ) and satisfies
(∇ L )(u˜ i , v˜ j ) = 0,
(45)
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
117
where
n b ∂ L= b f i (s)(u i , v i ) ds ∂s a
i =0
is the length of the curve c (s) =
n
i =0
(46)
f i (s)(u i , v i ).
Although all properties of geodesic-like curve in Definition 4 are still kept, the derivatives of L with respect to the parameters u i and v j are more complicated than that of E. Hence, the geodesic-like curve in Definition 4 is more effective than the geodesic-like curve in Definition 6. Finding the order of geodesic-like curve arriving that makes errors is less than ε for some ε > 0 is still a difficult task. We only know the optimum order depends on the parametrization. For example, when the parametrization is the exponential map at p on surfaces, every geodesic-like curve with the starting point p is a geodesic on the surface and the optimum order is 2. Further investigation on the mathematical analysis of geodesic-like curves and its applications is required in the near future. References Do Carmo, M.P., 1976. Differential Geometry of Curves and Surfaces. Prentice-Hall, Englewood Cliffs, NJ. Do Carmo, M.P., 1992. Riemannian Geometry. Hotz, I., Hagen, H., 2000. Visualizing geodesics. In: Proceedings IEEE Visualization. Salt Lake City, UT, pp. 311–318. Kasap, E., Yapici, M., Talay Akyildiz, F., 2005. A numerical study for computation of geodesic curves. Applied Mathematics and Computation 171 (2), 1206– 1213. Marsden, J.E., Hoffman, M.J., 1993. Elementary Classical Analysis. W.H. Freeman and Company. Martinez, D., Velho, L., Carvalho, P.C., 2005. Computing geodesics on triangular meshes. Computer & Graphics 29, 667–675. Paluszny, M., 2008. Cubic polynomial patches through geodesics. Computer-Aided Design 40, 56–61. Polthier, K., Schmies, M., 1998. In: Hege, H.C., Polthier, H.K. (Eds.), Straightest Geodesics on Polyhedral Surfaces in Mathematical Visualization. SpringerVerlag, Berlin. Ravi Kumar, G.V.V., Shastry, K.G., Prakash, B.G., 2003a. Computing offsets of trimmed NURBS surfaces. Computer-Aided Design 35, 411–420. Ravi Kumar, G.V.V., Srinivasan, P., Devaraja Holla, V., Shastry, K.G., Prakash, B.G., 2003b. Geodesic curve computations on surfaces. Computer Aided Geometric Design 20 (2), 119–133. Sánchez-Reyes, J., Dorado, R., 2008. Constrained design of polynomial surfaces from geodesic curves. CAD 40, 49–55. Sprynski, N., Szafran, N., Localle, B., Biard, L., 2008. Surface reconstruction via geodesic interpolation. CAD 40, 480–492. 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. Tucker, C.L., 1997. Forming of advanced composites. In: Gutowski, T.G. (Ed.), Advanced Composites Manufacturing. Wiley, New York.
Contents lists available at ScienceDirect
Computer Aided Geometric Design www.elsevier.com/locate/cagd
Geodesic-like curves on parametric surfaces Sheng-Gwo Chen 1 Department of Applied Mathematics, National Chiayi University, Chia-Yi 600, Taiwan
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 5 March 2009 Received in revised form 31 August 2009 Accepted 1 October 2009 Available online 3 October 2009
Finding the geodesic on a surface is an important task in many fields. As the traditional discrete algorithms or numerical approach can only find a list of discrete points, this study aims to propose a new, elegant and accurate method for approximating geodesics on a surface. The proposed geodesic-like curve is a mathematical curve that approaches to the geodesic on surface when the order of the curve reaches to infinity. The geodesic-like curve is also compared to some other well-known methods in this study. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Geodesics play an important key of differential geometry and a powerful tool to improve many problems in CAGD, CAD/CAM, computational geometry, computer vision, etc. (Paluszny, 2008; Ravi Kumar et al., 2003a; Sánchez-Reyes and Dorado, 2008; Sprynski et al., 2008). A curve α on a regular surface is a geodesic if its osculating plane at each point on the curve contains the normal vector of the surface at that point. The osculating plane at α (t 0 ) is a plane generated by the tangent vector and the normal vector of α at α (t 0 ). A geodesic on surface can be also defined by a curve with zero geodesic curvature at each point or a solution of the system of geodesic equations (2). In the calculus of variation, the geodesic curve is a critical point of the length function or the energy function in Eqs. (6) and (7). Hence, the geodesics joining two points on surfaces have extreme (maximal or minimal) lengths and extreme energies. More details of geodesics can be found in Do Carmo (1976). Analytical approaches are quite complex and difficult to find on surfaces (Do Carmo, 1976). In 2000, Hotz and Hagen (2000) improved the geodesic from the first, second fundamental forms and numerical method. The numerical approach presented by Emin Kasap et al. (2005) for the computation of geodesics between two fixed points on a surface is elegant through the use of the system of geodesic equations in Eq. (2) and finite-difference. Many accurate discrete methods approach geodesics or the shortest paths on tessellated surfaces (Tucker, 1997), polygonal surfaces (Polthier and Schmies, 1998) and triangular meshes (Martinez et al., 2005; Surazhsky et al., 2005). These methods are listed in the references. A new method for estimating the geodesic is introduced in this paper. The geodesic-like curve which is a linear combination of base functions is a critical point of our energy function in Eq. (15). The base functions are considered as the Bernstein polynomials or B-spline functions in this paper. Obviously, the geodesic-like curve approaches the geodesic on surface when its order is large enough. The definition of geodesic-like curves on surfaces and some properties of geodesiclike curves are provided in Section 2. In Section 3, conditions to improve different kinds of geodesic problems are discussed and the method of geodesic-like curves with other recent researches.
1
E-mail address: [email protected]. Partially supported by NSC, Taiwan.
0167-8396/$ – see front matter doi:10.1016/j.cagd.2009.10.001
© 2009 Elsevier B.V.
All rights reserved.
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
Fig. 1. A proper variation of
107
α.
2. Geodesics and geodesic-like curves A geodesic on a Riemannian manifold ( M , g ) is a smooth curve
D ds
α = 0,
α (s) : [a, b] → M which satisfies the differential equation (1)
D d where ds is the covariant derivative associated to g and α = ds α (s). In a local coordinate system (x1 , . . . , xn ), geodesics can be represented by the solutions of the system of differential equations
d 2 xi dt 2
+
dxk dx j Γ jki x(t ) = 0 for each i ∈ {1, 2, . . . , n}. dt dt
j ,k
(2)
The system (2) is known as the differential equations of the geodesics of M. Now, we shall discuss geodesics on surfaces from the variation of calculus. Let S be a regular surface and let b be a parametrization on S, b : U ⊂ R2 → S. Hence, b is a diffeomorphism between U and b(U ). If α is a smooth curve on S with a parameter α : [a, b] → U , then the arc length of α is defined by
b d L (α ) = b α (s) ds ds
(3)
a
and the energy of
α is
b 2 1 d b α (s) ds. E (α ) = 2 ds
(4)
a
α , see Fig. 1, is a differentiable map h : [a, b] × [−ε, ε] → U such that h(s, 0) = α (s), s ∈ [a, b], (5) h(a, t ) = α (a), t ∈ [−ε , ε ], h(b, t ) = α (b), t ∈ [−ε , ε ]. Intuitively, a proper variation of α , ht = h(·, t ), is a family of curves which are differentiable on the parameter t ∈ (−ε , ε ) ∂(b◦h) ∂(b◦h) with the same initial point α (a) and the same end point α (b) and such that h0 agrees with α . Thus ∂ s and ∂ t are differentiable vector fields along α . The energy of the curve ht is represented as 2 b 1 ∂(b ◦ h) E (t ) = (s, t ) (6) ds, t ∈ [−ε , ε ], 2 ∂s A proper variation of
a
and the length of ht is represented as
L (t ) =
b ∂(b ◦ h) ds, ( s , t ) ∂s a
t ∈ [−ε , ε ].
(7)
108
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
From the first variation formula of E (t ),
E (0) =
b
D
V (t ),
dt
d
b◦α
dt
dt ,
(8)
a
where V (t ) is the variational field of b ◦ h, Theorem 1 is obvious. Theorem 1. Suppose that b : U → S is an parametrization on S and α is a curve parametrized by arc length on surface S. E (t ) and L (t ) are defined by Eqs. (6) and (7), respectively. The following statements are equivalents. (1) b(α (s)) is a geodesic on surface S. (2) α (s) has extreme energy on S, that is E (0) = 0. (3) α (s) has extreme length on S, i.e. L (0) = 0. Hence, the shortest path from p to q is a geodesic. We can find these properties in Do Carmo (1976, 1992). For simplicity,
α : [a, b] → U is called a geodesic on S if b(α ) is a geodesic on S in this study. According to the Weierstrass theorem (Marsden and Hoffman, 1993), any continuous real function can be approximated by the Bernstein polynomial. We summarize the Weierstrass theorem as follows. Theorem 2 (Weierstrass theorem). Let f : [0, 1] → R be a continuous function and let > 0. Then the sequence of Bernstein polynomials
p n (s) =
n
i
B ni (s) f
(9)
n
i =0
converges uniformly to f as n → ∞, where B ni (s) = C in (1 − s)n−i si is the Bernstein polynomial. Thus geodesics on surfaces can be approximated by the Bézier curves. Theorem 3. Let S be a parametric surface with a parametrization b : U → S and let K be a compact subset in U . Suppose that γ : [a, b] → K is a geodesic on S with no any conjugate point. If the curve c˜ n (s) = ni=0 B ni (s)˜c i in K is a solution of the optimal problem:
2
n b ∂ n min E (u 1 , . . . , un−1 , v 1 , . . . , v n−1 ) = B i (s)(u i , v i ) ds b 2 ∂s 1
s.t.
i =0
a
(u 0 , v 0 ) = γ (a), (un , v n ) = γ (b), (u i , v i ) ∈ K for each i ,
(10)
then limn→∞ c˜ n = γ . Proof. It is clear that {˜cn } is a uniform convergent sequence. Without the loss of generality, we assume that [a, b] = [0, 1]. By Theorem 2, a sequence of Bézier curves exists n
c (s) =
n
i
B ni (s)
γ
(11)
n
i =0
γ as n → ∞. As c˜ n is a solution of Eq. (10), we have n E (γ ) E c˜ E cn
converges uniformly to
(12)
for each n. In addition, {˜c } and {c } are uniform convergent sequences. Thus n
n
E (γ ) lim E c˜ n = E lim c˜ n lim E cn = E lim cn = E (γ ). n→∞
n→∞
n→∞
n→∞
(13)
This implies that E (limn→∞ c˜ n ) = E (γ ), that is limn→∞ c˜ n has minimal energy. Therefore, limn→∞ c˜ n is a geodesic between γ (a) and γ (b) on S. Because γ is the unique minimum geodesic between γ (a) and γ (b), limn→∞ c˜ n agrees with γ . The curve c˜ n in the above theorem is called a geodesic-like curve. The definition of the geodesic-like curve is stated as follows.
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
109
Definition 4. Let b(u , v ) be a parametrization of a parametric surface S, b : U ⊂ R2 → S and let { f i (s)| f i : [a, b] → R}ni=0 be a n sequence of real functions. A curve c˜ (s) on U is called a geodesic-like curve of order (n + 1) on S if c˜ (s) = i =0 f i (s)(u˜ i , v˜ i ) and satisfies
(∇ E )(u˜ i , v˜ j ) = 0,
(14)
where
2
n b ∂ E= b f i (s)(u i , v i ) ds 2 ∂s 1
(15)
i =0
a
is the energy of the curve c (s) =
n
i =0
f i (s)(u i , v i ).
The function f i (s) is called the base function of geodesic-like curve and c (s) is the base curve of geodesic-like curve. We know that any continuous function can be approximated by the Bernstein polynomials (see Theorem 2) or the B-spline functions. If f i (s) is a Bernstein polynomial or B-spline function, the geodesic-like curve approaches a geodesic on S when its order n is large enough. The base function f i is assumed a Bernstein polynomial or B-spline function in this paper. From Theorem 3, the following theorem is obvious. Theorem 5. The solutions of the following optimization problem are geodesic-like curves:
1 d 2 min E (u i , v j ) = b c (s) ds 2 ds 1
0
⎧ ⎨ c (s) = ni=0 f i (s)(u i , v i ) ∈ U , s.t. c (a) = x0 , ⎩ c (b) = x1 .
(16)
Moreover, the geodesic-like curve approaches the shortest path between x0 and x1 . Its length approaches the distance between x0 and x1 when n is large enough. 2 3 Let us present the gradient and the hessian of the nenergy function E in Eq. (15). Denote b : U ⊂ R → R is a parametrization of a parametric surface S and c (s) = i =0 f i (s)(u i , v i ) is a base curve on U . For simplicity, b denotes b(c (s)), c denotes c (s) and the notation means the derivative of function with respect to the parameter s. The energy of c is the function
E (u i , v j ) =
1
b
2
b , b ds.
a
Set c (s) = (u (s), v (s)) where u (s) =
(17)
n
i =0
∂ u (s) = f i (s), ∂ ui ∂ u ( s ) = 0, ∂ vi ∂ v ( s ) = 0, ∂ ui ∂ v (s) = f i (s), ∂ vi
f i (s)u i and v (s) =
n
i =0
f i (s) v i . Since
(18)
we obtain
∂ ∂ ∂ b = b u (s), v (s) ∂ ui ∂ ui ∂ s ∂ (bu u + b v v ) = ∂ ui = buu u f i + bu f i + buv v f i = (buu u + buv v ) f i + bu f i
(19)
110
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
and
∂ b = (buv u + b v v v ) f i + b v f i . ∂ vi Hence,
E ui = E vi =
b
(buu u a b (buv u a
+ buv v ) f i + bu f i , bu u + b v v ds, + b v v v ) f i + b v f i , bu u + b v v ds.
(20)
(21)
The gradient of E is
∇ E = ( E u 0 , E u 1 , . . . , E un , E v 0 , E v 1 , . . . , E v n ).
(22)
Therefore, c (s) is a geodesic-like curve if it is a solution of
b
(buu u + buv v ) f i + bu f i , bu u + b v v ds = 0,
a
b
(buv u + b v v v ) f i + b v f i , bu u + b v v ds = 0,
(23)
a
for each i ∈ {0, 1, . . . , n}. The system of equations in (15) or (23) is called the system of geodesic-like equations. Since
b E ui = a
b Ev j = a
∂ b (t ), b (t ) dt , ∂ ui ∂ ∂v j
b (t ), b (t ) dt ,
(24)
we have
⎧ b ∂ 2 ∂ ∂ ⎪ ⎪ E u i u j = a ∂ u i ∂ u j b (t ), b (t ) + ∂ u i b (t ), ∂ u j b (t ) dt , ⎪ ⎪ ⎪ b ∂ 2 ∂ ⎪ ∂ ⎨E ui v j = a ∂ u i ∂ v j b (t ), b (t ) + ∂ u i b (t ), ∂ v j b (t ) dt , ⎪ E vi u j = Eu j vi , ⎪ ⎪ ⎪ ⎪ b ∂ 2 ∂ ⎪ ⎩E ∂ vi v j = a ∂ v i ∂ v j b (t ), b (t ) + ∂ v i b (t ), ∂ v j b (t ) dt .
(25)
By Eqs. (19) and (20), we have
⎧ ∂ 2 b ⎪ ∂ u i ∂ u j = (buuu u + buuv v ) f i f j + buu ( f i f j + f i f j ) ⎪ ⎪ ⎪ ⎪ ⎪ = (buuu u + buuv v ) A i j + buu ( A i j ) , ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ∂ 2 b = (buuv u + buv v v ) f i f j + buv ( f i f + f f j ) j i ∂ ui ∂ v j ⎪ = (buuv u + buv v v ) A i j + buv ( A i j ) , ⎪ ⎪ ⎪ ⎪ ⎪ ∂ 2 b ⎪ ⎪ ∂ v i ∂ v j = (buv v u + b v v v v ) f i f j + b v v ( f i f j + f i f j ) ⎪ ⎪ ⎪ ⎩ = (buv v u + b v v v v ) A i j + b v v ( A i j ) ,
(26)
where A i j = f i f j . Thus the hessian of E is the matrix
E uu E vu
E uv Evv
,
(27)
where E uu = ( E u i u j ), E uv = E vu = ( E u i v j ) and E v v = ( E v i v j ). The constant curve is a trivial solution of the system of geodesic-like equations in Eq. (23), because the nonzero vector ω = (0, u 1 − u 0 , u 2 − u 0 , . . . , un − u 0 , v 1 − v 0 , . . . , v n − v 0 ) is always in a decreasing direction. In next section, different kinds of initial conditions of the system of geodesic-like equations are discussed and the proposed method is compared with others.
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
111
3. Some elegant methods and geodesic-like curves In this section, some elegant and accurate methods of solving geodesic problems are presented and compared with the proposed method. 3.1. Kasap’s numerical method Finding the shortest path between two points on surfaces is a classical problem and it has many important applications (Paluszny, 2008; Sprynski et al., 2008; Sánchez-Reyes and Dorado, 2008). A numerical research presented by Kasap et al. (2005) for the computation of geodesic between two points on surfaces is elegant and accurate. They solve the system (2) via the finite-difference method and the iterative method. By the centered-difference formulas, the differential equations of the geodesics
⎧ 2 ⎨ d 2u + Γ 1 du 2 + Γ 1 du dv + Γ 1 dv 2 = 0, 11 ds 12 ds ds 22 ds ds ⎩ d2 v + Γ 2 du 2 + Γ 2 du dv + Γ 2 dv 2 = 0 11 ds 12 ds ds 22 ds ds2
(28)
can be regard as the system
U i +1 − 2U i + U i −1 = h2 f (·), V i +1 − 2V i + V i −1 = h2 g (·),
(29)
where
1 f (·) = −Γ11 (U i , V i )
U i +1 − U i −1
2
2h
U i +1 − U i −1 V i +1 − V i −1 1 − Γ12 (U i , V i )
2 V i +1 − V i −1 1 − Γ22 (U i , V i )
2h
2h
(30)
2h
and
2 g (·) = −Γ11 (U i , V i )
U i +1 − U i −1 2h
2
U i +1 − U i −1 V i +1 − V i −1 2 − Γ12 (U i , V i )
2 V i +1 − V i −1 2 − Γ22 (U i , V i ) .
2h
2h
2h
(31)
The system (29) is a system of nonlinear equations with parameters {U i , V j }. It can be approximated by the Newton’s method, the iterative method or other numerical methods. In Figs. 2 and 3, the yellow curves are the Kasap’s geodesics between two points on the Witeny’s umbrella and a spherical surface. Let b be a parametrization on surface S and let p , q be two distinct points on S. There exists two points (u˜ 0 , v˜ 0 ), (u˜ n , v˜ n ) in U such that
b(u˜ 0 , v˜ 0 ) = p , b(u˜ n , v˜ n ) = q.
(32)
We consider the base curve c : [a, b] → U in the Definition 4 satisfies c (a) = (u˜ 0 , v˜ 0 ) and c (b) = (u˜ n , v˜ n ), i.e. the starting point and the terminal point of c are fixed. As E u 0 , E v 0 , E un and E v n are vanished, the system of geodesic-like equations in Eq. (15) becomes
∇ E ( u i , v j ) = ( E u 1 , E u 2 , . . . , E u n −1 , E v 1 , E v 2 , . . . , E v n −1 ) = 0.
(33)
The solutions of the system (33) are the geodesics between p and q. Furthermore, if the hessian of E is a positive-definite matrix, then this geodesic-like curve approaches the shortest path between p and q as n → ∞. The system (33) is also a system of nonlinear equations. This system can certainly improved by the Newton’s method, the fixed point theorem or other numerical methods. The Kasap’s geodesics and the geodesic-like curves on the Witeny umbrella and the spherical surface obtained are show in Figs. 2 and 3. The Kasap’s geodesics are yellow curves and the geodesic-like curves are red curves. In Fig. 3, the black curve is an exact geodesic on the spherical surface. In these simulations, the Kasap’s geodesic has 100 vertices and the geodesic-like curve comes from a cubic B-spline curve with 10 control points in U . In Fig. 2, the length of geodesic-like curve is less than the length of Kasap’s geodesic. The geodesic-like curve is closer to the exact geodesic than the Kasap’s geodesic on spherical surface and as shown in Fig. 3. Both of Kasap and our methods study the geodesic problem via a system of nonlinear equations. The number of control points of geodesic-like curve is less than the number of vertices in Kasap’s geodesic and the geodesic-like curve is more accurate than the Kasap’s geodesic.
112
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
Length Time (second)
Kasap’s geodesic (yellow curve)
Geodesic-like curve (red curve)
1.553 2.012
1.550 0.354
Fig. 2. The Kasap’s geodesic and geodesic-like curve on Witeny umbrella. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
Time (second)
Kasap’s geodesic
Geodesic-like curve
4.21
0.407
Fig. 3. The Kasap’s geodesic and geodesic-like curve on spherical surface. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
3.2. Ravi Kumar’s discrete method In this paragraph, we shall study the problem of starting point and its direction. It is well known that for each p ∈ S and y ∈ T p S, there exists an unique geodesic such that α (0) = p and α (0) = y locally (Do Carmo, 1976). Many researches improve this problem by a polygon and a triangular mesh rather than a regular surface. Martinez et al. (2005), Surazhsky et al. (2005) and Ravi Kumar et al. (2003b) present some different discrete methods for solving this problem. All of these methods are accurate. The Ravi Kumar’s method is described and compared with the proposed method. The main idea of Ravi Kumar’s method is the orthogonal projection. The method 2b in Ravi Kumar’s paper (Ravi Kumar et al., 2003b) is summarized as below. (1) Given a point p 0 and a tangent direction v 0 . (2) Add p 0 at the first point of discrete geodesic P. For each i 0 (3) Determine the point p i whether a vertex, in a edge or in a triangle and compute the normal vector ni at p i . (4) Let t i be the unit normalized projection of the vector v i onto some facet on the triangular mesh or tessellation and shown in Fig. 4. (5) Determine the next point p i +1 by the intersect of the line from p i with direction t i and the edges of triangle.
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
113
Fig. 4. The next geodesic direction t i .
Although Ravi Kumar’s method works on a tessellation or a triangular mesh, it can be extended to estimate geodesics on regular surfaces easily. It only needs to project each point of discrete geodesic to the regular surface via the normal vector field on surface S. We can also improve this problem by geodesic-like curve. Given a point p ∈ U and a vector y ∈ R2 such that p + y ∈ U . Consider the conditions
c (0) = p ∈ U , c (0) = y ∈ R2 .
(34) and c (0)
Suppose that c (s) is a base curve defined on [0, 1], we have c (0) = (u 0 , v 0 ) = n((u 1 , v 1 ) − (u 0 , v 0 )). Thus (u 0 , v 0 ) = p and (u 1 , v 1 ) = ny + p are fixed in U . The system of geodesic-like equations becomes
∇ E ( u i , v j ) = ( E u 2 , E u 3 , . . . , E u n , E v 2 , E v 3 , . . . , E v n ) = 0.
(35)
However, these equations have a trivial solution (u i , v i ) = (u 1 , v 1 ) for each i ∈ {2, 3, . . . , n}. To avoid this problem, we set the terminal point of c (s) fixed, the system of geodesic-like equations becomes
∇ E ( u i , v j ) = E ( u 2 , u 3 , . . . , u n −1 , v 2 , v 3 , . . . , v n −1 ) = 0.
(36)
If the set {u i }ni=2 or the set { v i }ni=2 is fixed, then the system of geodesic-like equations correspond to
∇ E (u i , v j ) = E ( v 2 , v 3 , . . . , v n ) = 0
(37)
∇ E ( u i , v j ) = E ( u 2 , u 3 , . . . , u n ) = 0.
(38)
or
Ravi Kumar’s geodesics and the geodesic-like curves on the Witeny umbrella and a torus are estimated to have 140 vertices and quadric B-spline curves with 20 control points respectively in our simulations. The Ravi Kumar’s geodesics are yellow curves and the geodesic-like curves are red curves in Figs. 5 and 6. In Figs. 5 and 6, the number of control points and the length of geodesic-like curve are fewer and shorter than the Ravi Kumar’s geodesic. Thus geodesic-like curve is more accurate that Ravi Kumar’s geodesic. 3.3. Geodesic between two curves The problem of the shortest path between two curves on surfaces is also important but most geodesic estimating methods can not improve this problem easily. We shall present the geodesic-like model for solving this problem. Let S be a parametric surface and let b : U → S be a parametrization on S. Assume that Γ1 and Γ2 are two curves on S with parametrization c 1 : [0, 1] → U and c 2 : [0, 1] → U , respectively. The mathematical model of the distance between Γ1 and Γ2 on S is
min E (c ) such that c : [a, b] → U , c ([a, b]) ∩ c 1 ([0, 1]) = c (a), c ([a, b]) ∩ c 2 ([0, 1]) = c (b).
(39)
114
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
Length Time (second)
Ravi Kumar’s geodesic (yellow curve)
Geodesic-like curve (red curve)
10.710 2.023
10.706 2.012
Fig. 5. The discrete geodesic and geodesic-like curve on Witeny umbrella surface. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
The length of exact geodesic (black curve) = 22.8828516
Ravi Kumar’s geodesic Geodesic-like curve
Length
Time (second)
22.9259979 22.8828616
2.047 2.051
Fig. 6. The discrete geodesic and geodesic-like curve on torus. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
Denote c 1 = c 1 (μ) = (u 0 (μ), v 0 (μ)) and c 2 = c 2 (ν ) = (un (ν ), v n (ν )). If the curve c is the base curve in Definition 4, then E is a real function with parameters μ, ν , u 1 , u 2 , . . . , un−1 and v 1 , v 2 , . . . , v n−1 . Hence, the solutions of (39) satisfy the system of geodesic-like equations
∇ E = ( E t 1 , E t 2 , E u 1 , E u 2 , . . . , E u n −1 , E v 1 , E v 2 , . . . , E v n −1 ) = 0.
(40)
where E μ and E ν can be computed by
E μ = E u0 E ν = E un
du 0 dμ dun dν
+ E v0 + E vn
dv 0 dμ dv n dν
, .
(41)
In other words, Eq. (40) is a system of geodesic-like equations with the conditions
c (0) ∈ C 1 , c (1) ∈ C 2 .
(42)
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
115
Fig. 7. The geodesic-like curve on closed surfaces (periodic surfaces). (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
Fig. 8. The minimal geodesic-like curve on a surface of revolution. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
Fig. 9. The closed geodesic on a surface of revolution. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
4. Conclusions and discussions Let S be a closed surface (periodic surface) with a parametrization b : U → S. As b is not an onto map, the geodesic-like curves should not defined on U . In this case, we consider the periodic map
b˜ : R2 → S
(43)
such that b = b˜ ◦ π where π is the projection map from R2 to U . And construct the geodesic-like curve by the energy function
E=
b ∂ 2 ˜ ◦ c (s) ds. b 2 ∂s 1
a
In this case, geodesic-likes curve should be defined on R2 . Geodesics on periodic surfaces are shown in Figs. 7–9.
(44)
116
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
Fig. 10. The geodesic-like curves with different number of control points. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
Number of control point Length of geodesic-like curve Time (second) Error Number of control point Length of geodesic-like curve Time (second) Error
4
5
6
7
8
31.1151331 0.04125 0.2318035 9 30.900999 0.323 0.0176694
30.9761609 0.0813 0.0928313 10 30.8892435 0.4404 0.0059139
30.9755519 0.13 0.0922223 20 30.8853687 2.64 0.0020391
30.9040219 0.18 0.0206923 40 30.8842772 6.39 0.0009476
30.9035892 0.234 0.0202596 60 30.8833296 20.12 <1.0e−7
Fig. 11. The length of geodesic-like curves by quadric B-spline curves.
The degree of geodesic-like curve Length of geodesic-like curve Time (second) Error
2
3
4
5
6
7
10
31.1151 0.015 0.2318
31.0047 0.062 0.12134
30.9756 0.219 0.09222
0.964 0.344 29.9193
30.9156 0.469 0.03231
30.9048 0.672 0.02151
30.8952 2.45 0.01184
Fig. 12. The length of geodesic-like curves by Bézier curves.
The geodesic-like curves in this section always come from uniform quadric B-spline curves and Bézier curves. The geodesic-like curves on periodic surfaces (closed surface) obtained and shown in Figs. 7, 8 and 9. The yellow curves on these figures are the boundary of surfaces and the red curves are the geodesic-like curves which come from the uniform quadric B-spline curves with 8 control points. In Fig. 7, the geodesic-like curves are not the shortest path on a torus and a surface of revolution. In Fig. 8, the geodesic-like curve agrees with the shortest path on surface. However, all of geodesiclike curves are close to geodesics on that surface. We also show a closed geodesic on a surface of revolution in Fig. 9. In Fig. 10, we show some geodesic-like curves with different orders. The order of red geodesic-like curve is 4 and that of yellow geodesic-like curve is 10. We list the results of the length of geodesic-like curves with different orders, the number of control points, in Fig. 11. Obviously, the geodesic-like curve approaches the minimal geodesic when the order is large enough. In Fig. 12, the geodesic-like curves are the Bézier curves. Similarly, the geodesic-like curve approaches the geodesic when the degree increasing. Because it is expensive to estimate the Bézier curve of large degree, the geodesic-like curve by the B-spline curve is more effective the Bézier curve. In these tests, the length of exact geodesic curve is 30.8833296. It is clear that the geodesic-likes curve also converges to the exact geodesic. From these figures, geodesic-like curves come from the uniform quadric B-spline curves are more accurate than the Bézier curves. The method of the geodesic-like curve is different from other traditional approaches. The advantages of geodesic-like curve are nature, accuracy and elegance. Furthermore, the geodesic-like curve is a mathematical curve. However, the accuracy of geodesic-like curve depends on different parametrizations on surface. It is well known that the critical points of length function in Eq. (7) agree with the critical points of energy function in Eq. (6) when the curves always have arc-length parameters (Do Carmo, 1976, 1992). Therefore, we can also define the geodesic-like curve via length function. Another definition is stated as follows. Definition 6 (Another definition of geodesic-like curve). Let b(u , v ) be a parametrization of a parametric surface S, b : U ⊂ R2 → S and let { f i (s) | f i : [a, b] → R}ni=0 be a sequence of real functions. A curve c˜ (s) on U is called a geodesic-like curve n of order (n + 1) on S if c˜ (s) = i =0 f i (s)(u˜ i , v˜ i ) and satisfies
(∇ L )(u˜ i , v˜ j ) = 0,
(45)
S.-G. Chen / Computer Aided Geometric Design 27 (2010) 106–117
117
where
n b ∂ L= b f i (s)(u i , v i ) ds ∂s a
i =0
is the length of the curve c (s) =
n
i =0
(46)
f i (s)(u i , v i ).
Although all properties of geodesic-like curve in Definition 4 are still kept, the derivatives of L with respect to the parameters u i and v j are more complicated than that of E. Hence, the geodesic-like curve in Definition 4 is more effective than the geodesic-like curve in Definition 6. Finding the order of geodesic-like curve arriving that makes errors is less than ε for some ε > 0 is still a difficult task. We only know the optimum order depends on the parametrization. For example, when the parametrization is the exponential map at p on surfaces, every geodesic-like curve with the starting point p is a geodesic on the surface and the optimum order is 2. Further investigation on the mathematical analysis of geodesic-like curves and its applications is required in the near future. References Do Carmo, M.P., 1976. Differential Geometry of Curves and Surfaces. Prentice-Hall, Englewood Cliffs, NJ. Do Carmo, M.P., 1992. Riemannian Geometry. Hotz, I., Hagen, H., 2000. Visualizing geodesics. In: Proceedings IEEE Visualization. Salt Lake City, UT, pp. 311–318. Kasap, E., Yapici, M., Talay Akyildiz, F., 2005. A numerical study for computation of geodesic curves. Applied Mathematics and Computation 171 (2), 1206– 1213. Marsden, J.E., Hoffman, M.J., 1993. Elementary Classical Analysis. W.H. Freeman and Company. Martinez, D., Velho, L., Carvalho, P.C., 2005. Computing geodesics on triangular meshes. Computer & Graphics 29, 667–675. Paluszny, M., 2008. Cubic polynomial patches through geodesics. Computer-Aided Design 40, 56–61. Polthier, K., Schmies, M., 1998. In: Hege, H.C., Polthier, H.K. (Eds.), Straightest Geodesics on Polyhedral Surfaces in Mathematical Visualization. SpringerVerlag, Berlin. Ravi Kumar, G.V.V., Shastry, K.G., Prakash, B.G., 2003a. Computing offsets of trimmed NURBS surfaces. Computer-Aided Design 35, 411–420. Ravi Kumar, G.V.V., Srinivasan, P., Devaraja Holla, V., Shastry, K.G., Prakash, B.G., 2003b. Geodesic curve computations on surfaces. Computer Aided Geometric Design 20 (2), 119–133. Sánchez-Reyes, J., Dorado, R., 2008. Constrained design of polynomial surfaces from geodesic curves. CAD 40, 49–55. Sprynski, N., Szafran, N., Localle, B., Biard, L., 2008. Surface reconstruction via geodesic interpolation. CAD 40, 480–492. 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. Tucker, C.L., 1997. Forming of advanced composites. In: Gutowski, T.G. (Ed.), Advanced Composites Manufacturing. Wiley, New York.