Existence conditions for Coons patches ... - Semantic Scholar
Computer Aided Geometric Design 26 (2009) 599–614
Contents lists available at ScienceDirect
Computer Aided Geometric Design www.elsevier.com/locate/cagd
Existence conditions for Coons patches interpolating geodesic boundary curves R.T. Farouki b , N. Szafran a , L. Biard a,∗ a b
Laboratoire Jean Kuntzmann, Université Joseph Fourier – Grenoble, France Department of Mechanical and Aeronautical Engineering, University of California, Davis, CA 95616, USA
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 19 March 2008 Received in revised form 10 October 2008 Accepted 1 January 2009 Available online 7 January 2009 Keywords: Geodesic curves Surface reconstruction Hermite/Coons interpolation Surface smoothing
Given two pairs of regular space curves r1 (u ), r3 (u ) and r2 ( v ), r4 ( v ) that define a curvilinear rectangle, we consider the problem of constructing a C 2 surface patch R(u , v ) for which these four boundary curves correspond to geodesics of the surface. The possibility of constructing such a surface patch is shown to depend on the given boundary curves satisfying two types of consistency constraints. The first constraint is global in nature, and is concerned with compatibility of the variation of the principal normals along the four curves with the normal to an oriented surface. The second constraint is a local differential condition, relating the curvatures and torsions of the curves meeting at each of the four patch corners to the angle between those curves. For curves satisfying these constraints, the surface patch is constructed using a bicubically-blended Coons interpolation process. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Given four parametric space curves r1 (u ), r3 (u ) and r2 ( v ), r4 ( v ) specifying a curvilinear rectangle, the Coons interpolation procedure (Coons, 1964, 1974; Farin, 2002; Gordon, 1983) defines a surface patch R(u , v ) bounded by these four curves, such that R(u , 0) = r1 (u ), R(u , 1) = r3 (u ) and R(0, v ) = r2 ( v ), R(1, v ) = r4 ( v ). The Coons scheme is motivated by the fact that, in numerous surface design or reconstruction contexts, only the boundary curves of a surface patch are specified, and a means of smoothly “filling in” the interior is needed. In addition to the four boundary curves, the bicubically-blended Coons patch requires transverse derivative data along them — i.e., the four vector functions R v (u , 0), R v (u , 1) and Ru (0, v ), Ru (1, v ) must also be specified. In the Coons interpolation scheme, the boundary curves and derivative data (which together specify the tangent plane variation along the patch boundary) are considered to impart sufficient information to ensure the desired surface shape. In this paper, the Coons patch is modified to admit a more fundamental geometrical significance to the specified boundary curves. Namely, these curves are stipulated to be geodesics on the constructed surface. The motivation for this study comes from the problem of constructing analytic computer representations of free-form surfaces from positional/orientational measurements obtained with the Morphosense — a flexible ribbon-like device with embedded microsensors that, when placed on a physical surface, assumes the shape of a geodesic (Sprynski et al., 2008). By placing the Morphosense on the surface at regular intervals along two different directions, it is subdivided into rectangular patches with geodesic boundary curves.
*
Corresponding author. E-mail address: [email protected] (L. Biard).
0167-8396/$ – see front matter doi:10.1016/j.cagd.2009.01.003
© 2009
Elsevier B.V. All rights reserved.
600
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
The construction of surfaces that incorporate one or two given space curves as geodesics has been considered by several authors (Bennis et al., 1991; Sánchez-Reyes and Dorado, 2008; Tucker, 1997; Wang et al., 2004), in the context of applications such as distortion-free mapping of textures onto free-form surfaces; specifying fabric shapes for garment and shoe design; and in the layout of fibers in composite material structures. However, the problem of constructing rectangular surface patches, when all four prescribed boundary curves are required to be geodesics of the resulting surface, does not appear to have been previously studied. This paper is concerned with extending prior studies that address the problem of interpolating two boundary geodesics by a surface “strip” (Paluszny, 2008; Sprynski et al., 2008) to the case where all four boundary curves of a rectangular patch are specified as geodesics of the constructed surface. In particular, we are concerned with identifying constraints on the boundary curves, whose satisfaction constitutes a sufficient-and-necessary condition for the existence of analytic surfaces that interpolate those curves as geodesics. A companion paper (Farouki et al., 2008) will address the construction of (polynomial or rational) Bézier surface patches that interpolate four given Bézier curves, satisfying the existence conditions, as geodesic boundaries. The plan for this paper is as follows. After reviewing basic concepts concerning the geometry of curves on surfaces in Section 2, the ideas of Sprynski et al. (2008) are extended in Section 3 to identify a local condition that must be satisfied by two intersecting space curves if those curves are to be geodesics of an analytic surface. A global compatibility condition, that must be satisfied by the geodesic boundary curves of a rectangular patch to ensure a continuous oriented surface normal, is then identified in Section 4. Section 5 describes a modified Coons algorithm, for boundary data satisfying the conditions of Sections 3 and 4, to construct surface patches with the given boundary curves as geodesics. Finally, Section 6 presents some representative surfaces computed using this algorithm, and Section 7 summarizes our main results and identifies key issues that warrant further investigation. 2. Background on curves and surfaces In the following discussion, all curves and surfaces are considered to be regular and “sufficiently smooth.” A curve is regular if it admits a tangent line at each point, while a surface is regular if it admits a tangent plane at each point. All surfaces are considered to be oriented. A surface is said to be oriented if its unit normal vector (4) is continuous on each closed regular curve on the surface. The inner product of two vectors u, v in R3 is denoted by u, v. Similarly, the plane through a point p in R3 spanned by two linearly-independent vectors u, v is denoted by [p, u, v]. For linearly-independent unit vectors u, v and a unit vector n such that n ⊥ u and n ⊥ v, we denote by (u, v)n the oriented angle between u and v in the sense of n. Precisely, the angle A = (u, v)n is defined (see Fig. 1) by sin A = det(u, v, n),
cos A = u, v.
(1)
The variable s is employed to denote arc length along a space curve. Note that the arc-length parameterization r : s → r(s) of a curve satisfies r (s) = 1 and r (s) ⊥ r (s) for all s. However, in this paper, a general parameterization r : t → r(t ) is often used in the surface construction problem. The parameters of functions may sometimes be omitted when no confusion can arise.
• With each point r(s) of a curve satisfying r (s) = 0, we associate the Serret–Frenet frame (e(s), n(s), b(s)) where e(s) = r (s), n(s) = r (s)/r (s), and b(s) = e(s)× n(s) are, respectively, the unit tangent, principal normal, and binormal vectors of the curve at the point r(s). The arc-length derivative of the Serret–Frenet frame is governed by the relations e( s ) 0 k(s) 0 e( s ) d (2) n(s) = −k(s) 0 τ (s) n(s) , ds b( s ) 0 −τ (s) 0 b( s ) where the curvature k(s) and torsion
k(s) = r (s)
and
τ (s) =
τ (s) of the curve r(s) are defined by
det(r (s), r (s), r (s))
r (s)2
.
Fig. 1. Angle measurement.
(3)
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
601
The osculating plane at each curve point r(s) is spanned by the two vectors e(s), n(s) and does not depend on the curve parameterization. If k(s) = 0 for some s, then r (s) = 0 and the normal vector n(s) and osculating plane are undefined at that point. This condition identifies an inflection of the curve. • On a regular oriented surface (u , v ) → R(u , v ), the unit normal is defined at each point in terms of the partial derivatives Ru = ∂ R/∂ u, R v = ∂ R/∂ v by N( u , v ) =
Ru (u , v ) × R v (u , v ) . Ru (u , v ) × R v (u , v )
(4)
• Consider a curve r(s) = R(u (s), v (s)) on a surface R(u , v ), where s denotes arc length for the space curve r(s), but not necessarily for the plane curve defined by s → (u (s), v (s)). With each point r(s), we associate the Darboux frame (e(s), h(s), N(s)) — where e(s) is the unit tangent vector of the curve, N(s) is the unit normal vector of the surface at the point R(u (s), v (s)) = r(s), and h(s) = N(s) × e(s). The arc-length derivative of the Darboux frame is given by the relations d
ds
e( s ) h(s) N( s )
=
k g (s) kn (s) 0 0 −τ g (s) −k g (s) 0 −kn (s) τ g (s)
e( s ) h(s) , N( s )
(5)
which define the normal curvature kn (s), the geodesic curvature k g (s), and the geodesic torsion curve r(s) as
kn =
de ds
,N ,
kg =
de ds
,h ,
τg =
dN ds
,h .
τ g (s) at each point of the (6)
• A regular curve r(t ) on a surface R(u , v ) is called a geodesic of the surface if its geodesic curvature is identically zero. From (2) and (5), this is equivalent to requiring that N(s) = ±n(s),
h(s) = ∓b(s)
(7)
— i.e., the Frenet and Darboux frames agree modulo signs. Hence, we have the following useful characterizations of geodesic curves. A regular curve t → r(t ) is a geodesic on the surface R(u , v ) if and only if (D1) the geodesic curvature of r(t ) is identically zero; (D2) the principal normal at each non-inflection point of r(t ) is orthogonal to the surface tangent plane at the point R(u (t ), v (t )) = r(t ); (D3) the osculating plane at each non-inflection point of r(t ) is orthogonal to the surface tangent plane at the point R(u (t ), v (t )) = r(t ). In the case of an isolated inflection point of the curve, we may encounter cases (Do Carmo, 1976) where the osculating plane cannot be defined by continuity, since the limit planes approaching it from the left and the right disagree. For nonisolated inflection points, the vanishing of the curvature over an interval in the parameter implies that the corresponding curve segment degenerates to a straight line (and hence a geodesic) on the surface. For such a segment, the Darboux frame can be defined at each point, but the Frenet frame is undefined. 3. Geodesic curves crossing on a regular surface The normal curvature kn and geodesic torsion τ g at any point of a curve r(t ) on a regular C 2 oriented surface R(u , v ) depend only on the tangent direction of r(t ) at that point. Namely, if Πp is the tangent plane at a given point p of the surface R(u , v ), and L is any line in Πp passing through p, all curves on R(u , v ) that go through p, and have L as their tangent line there, possess the same normal curvature kn and geodesic torsion τ g . This property allows us to define the normal curvature K n ( L ) and geodesic torsion T g ( L ) at each point p of a surface R(u , v ), in the direction of any line L in the tangent plane Πp . At a fixed surface point p, the normal curvature K n ( L ) exhibits two extrema K 1 , K 2 (the principal curvatures), corresponding to two orthogonal directions L 1 , L 2 (the principal directions). Consider an orthonormal basis (e1 , e2 ) for Πp aligned with the principal directions L 1 , L 2 such that e1 × e2 defines the unit surface normal N at p. Then each line L in the tangent plane Πp corresponds to an angle α (with 0 α < π ), such that the normal curvature of the surface at point p can be regarded as a π -periodic function of the angle α . It can be shown (Do Carmo, 1976) that K n (α ) = K 1 cos2 α + K 2 sin2 α . Similarly, the geodesic torsion of the surface at point p can be regarded as a given (Do Carmo, 1976) by
T g (α ) = ( K 1 − K 2 ) sin α cos α .
(8)
π -periodic function T g (α ) of the angle α , (9)
602
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
Remark. A point p at which K 1 = K 2 is called an umbilic of the surface. At such a point, the normal curvature K n (α ) defined by (8) is the same in all directions α (so the surface shape is locally like a sphere), while the geodesic torsion (9) vanishes identically. For subsequent use, we now derive Proposition 3, which identifies a necessary condition for two intersecting space curves to be geodesics on a regular surface. First, notice that we have only three degrees of freedom in the curvature distribution, namely the two principal curvatures K 1 , K 2 and the orientation of the principal directions. Thus, the normal curvature and geodesic torsion are locally determined in any direction from relations (8) and (9). So, the normal curvatures and geodesic torsions along two directions in a point of a regular surface must satisfy one compatibility equation, which is established in the following lemma. Lemma 1. For two directions in the surface tangent plane Πp specified by the angles torsion satisfy the relation
α1 , α2 the normal curvature and the geodesic
sin(α2 − α1 ) T g (α1 ) + T g (α2 ) = cos(α2 − α1 ) K n (α1 ) − K n (α2 ) .
(10)
Proof. By direct computation, we have
sin(α2 − α1 ) T g (α1 ) + T g (α2 )
= ( K 1 − K 2 )(sin α2 cos α1 − sin α1 cos α2 )(sin α1 cos α1 + sin α2 cos α2 )
= ( K 1 − K 2 ) sin α1 sin α2 cos2 α1 − cos2 α2 + cos α1 cos α2 sin2 α2 − sin2 α1 1 = ( K 1 − K 2 ) sin α1 sin α2 (cos 2α1 − cos 2α2 ) + cos α1 cos α2 (cos 2α1 − cos 2α2 ) = =
2 1 2 1 2
( K 1 − K 2 ) cos(α2 − α1 )(cos 2α1 − cos 2α2 )
cos(α2 − α1 ) K 1 2 cos2 α1 − 1 − 2 cos2 α2 + 1 − K 2 1 − 2 sin2 α1 − 1 + 2 sin2 α2
= cos(α2 − α1 ) K n (α1 ) − K n (α2 ) .
2
We now consider in greater detail the relations (7) that characterize a geodesic curve on a surface. Lemma 2. Consider a regular curve r(s) on a surface R(u , v ) with notations as introduced in Section 2. Then, setting N(s) = σ n(s) with σ = ±1, we have h(s) = −σ b(s),
k g (s) = 0,
τ g (s) = −τ (s),
kn (s) = σ k(s).
(11)
Proof. Setting N(s) = σ n(s) so that h(s) = −σ b(s), a direct calculation using relations (2) gives d ds
e( s ) h(s) N( s )
e( s ) 0 k(s) 0 = σ τ (s) 0 0 −σ b(s) = ds σ n(s) −σ k(s) 0 σ τ (s) 0 0 σ k ( s ) e( s ) = 0 0 τ (s) h(s) , 0 N( s ) −σ k(s) −τ (s) d
e( s ) n(s) b( s )
which corresponds to the stated results by comparison with relations (5).
(12)
2
Proposition 3. Consider two geodesics r1 (s) and r2 (s) parameterized by arc length on the surface R(u , v ), with principal normal n1 (s) and n2 (s), crossing at the point p = r1 (s1 ) = r2 (s2 ). Assuming that p is not an inflection on either of the curves r1 (s) and r2 (s), let α = (r 1 (s1 ), r 2 (s2 ))Np be the oriented angle between them at p, in the sense of the surface normal Np at that point. Also, let ki (s) and τi (s) be the curvature and torsion of ri (s) for i = 1, 2. Then, for the values σ1 , σ2 ∈ {−1, +1} such that Np = σ1 n1 (s1 ) = σ2 n2 (s2 ), we have
τ1 (s1 ) + τ2 (s2 ) sin α = σ2k2 (s2 ) − σ1k1 (s1 ) cos α .
(13)
Proof. Consider first any geodesic r(s), free of inflections, parameterized by arc length on the surface R(u , v ). Its geodesic curvature satisfies k g (s) ≡ 0, and its principal normal n(s) agrees (modulo sign) with the surface normal N(s) (here N(s) is the restriction of the unit surface normal to the curve r(s)). Setting N(s) = σ n(s) with σ = ±1, Lemma 2 gives
τ g (s) = −τ (s) and kn (s) = σ k(s). Hence, if the tangents to the two geodesics r1 (s) and r2 (s) at the point p are identified by angles the result follows directly from Lemma 1, with α = α2 − α1 . 2
(14)
α1 and α2 in the plane Πp ,
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
603
Fig. 2. Geodesics on a cylinder.
Proposition 4 (Consistency with surface orientation). The geodesic crossing relation (13) in Proposition 3 does not depend on the surface orientation. Proof. Consider a new parameterization of the surface in Proposition 3, that induces at each point a surface normal opposite to the original normal. Then the angle α and the values σ1 , σ2 will change sign, but the curvatures and the torsions, which are intrinsic to the curves r1 (s) and r2 (s), remain unchanged. Hence, we obtain
τ1 (s1 ) + τ2 (s2 ) sin(−α ) = (−σ2 )k2 (s2 ) − (−σ1 )k1 (s1 ) cos(−α ),
which is identical to (13).
(15)
2
Corollary 5. If two space curves intersecting at a point p do not satisfy the relation (13), no regular C 2 oriented surface can interpolate these curves in such a manner that they are geodesics of the surface. Example. Consider the cylinder x2 + y 2 = R 2 , and for a1 , a2 ∈ R (with a1 = a2 ) the geodesic curves (circular helices) defined by ri (t ) = ( R cos t , R sin t , ai t ) on it, crossing at p = r1 (0) = r2 (0) = ( R , 0, 0) with the derivatives r i (0) = (0, R , ai ) at that point. Let αi be the angle that r i (0) makes with the horizontal plane. Note that, for each circular helix, the principal normal n(t ) = (− cos t , − sin t , 0) points to the inside of the cylinder. From Proposition 4, we can choose an orientation of the cylinder such that σ1 = σ2 = 1, so that positive angles αi must be counted clockwise in Fig. 2. Thus, the positive angle α between the two vectors r i (0) is α = α1 − α2 , and we have tan α = tan(α1 − α2 ) =
a1 / R − a2 / R 1 + (a1 / R )(a2 / R )
=
R (a1 − a2 ) R 2 + a1 a2
For each helical curve ri (t ), the curvature ki (t ) and torsion direct computation gives k2 (0) − k1 (0)
τ1 (0) + τ2 (0)
=
R (a1 − a2 ) R 2 + a1 a2
= tan α =
sin α cos α
.
τi (t ) are constant, ki = R /(a2i + R 2 ) and τi = ai /(a2i + R 2 ), and a
,
which is relation (13). As another simple example, one can easily check that this relation is satisfied by great-circle geodesics crossing on a sphere. Proposition 6. (See Struik, 1976.) Through each point of a smooth surface, there exists a geodesic in every direction, and a geodesic is uniquely determined by an initial point and tangent direction at that point. Henceforth we shall consider only distinct crossing geodesics, with a non-zero angle of intersection. For such curves, the tangent plane Πp of the surface R(u , v ) at the point p is spanned by the two derivatives r 1 (s1 ), r 2 (s2 ). 4. Geodesic interpolation problem We now consider the problem of defining a smooth rectangular surface patch, bounded by four given parametric curves, in such a manner that these border curves are geodesics of the surface. This problem, introduced in Section 4.1, is called the “geodesic interpolation problem.” After introducing notations in Section 4.2 and identifying a local corner compatibility condition on the curve principal normals in Section 4.3, we introduce in Section 4.4 the global normal compatibility constraint, which ensures that the surface normal vector varies continuously around the patch boundary. Finally, Section 4.5 specifies the local differential constraints that must be satisfied by the border curves meeting at each patch corner, in order to ensure the existence of a surface for which these curves are geodesics. The details of the surface construction procedure will be presented in Section 5.
604
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
Fig. 3. Surface patch boundary curves.
4.1. Specification of the problem Consider, as illustrated in Fig. 3, four regular space curves r1 (u ), r3 (u ) and r2 ( v ), r4 ( v ) with u , v ∈ [0, 1] such that r1 (0) = r2 (0) = p00 , r1 (1) = r4 (0) = p10 , r2 (1) = r3 (0) = p01 , r3 (1) = r4 (1) = p11 .
(16)
As noted earlier, these curves are assumed to be “sufficiently smooth.” Also, consistent with the last remark in Section 3, their derivatives at each corner pi j , i , j = 0, 1, are assumed to be linearly independent, in order to define the tangent plane Πi j of the interpolating surface R(u , v ) (see below) at pi j . Our goal is to construct a surface that interpolates these four curves in such a manner that they are geodesics on the constructed surface. More precisely, if R(u , v ) for (u , v ) ∈ [0, 1]2 denotes the interpolating surface, it must exhibit the following three properties.
• South-North interpolating property: R(u , 0) = r1 (u ) and R(u , 1) = r3 (u ) for u ∈ [0, 1]; • West-East interpolating property: R(0, v ) = r2 ( v ) and R(1, v ) = r4 ( v ) for v ∈ [0, 1]; • Geodesic property: The principal normal at each non-inflectional point of the curves r1 (u ), r3 (u ) and r2 ( v ), r4 ( v ) is normal to the surface R(u , v ). Note that each of the geodesic curves r1 (u ), r3 (u ) and r2 ( v ), r4 ( v ) is also an isoparametric curve of the surface R(u , v ). 4.2. Notations Let t be either of the parametric variables u or v. For i = 1, 2, 3, 4 we denote by (ei (t ), ni (t ), bi (t )) and ki (t ), Serret–Frenet frame and the curvature and torsion of the curves ri (t ) at each non-inflection point. Namely, ei (t ) = ki (t ) =
r i (t )
||r i (t )||
,
bi (t ) =
r i (t ) × r i (t ) , r i (t )3
r i (t ) × r i (t )
||r i (t ) × r i (t )||
τi (t ) =
,
ni (t ) = bi (t ) × ei (t ),
det(r i (t ), r i (t ), r (t )) i
r i (t ) × r i (t )2
.
τi (t ) the
(17) (18)
At inflection points, the curvature is zero and the torsion is undefined. Note that, whereas the principal normal vector is intrinsic to the curve, the tangent and binormal vectors depend on the sense of the parameterization. 4.3. Condition (C1): osculating constraints at corners From the geodesic definitions in Section 2 we see that, at each corner pi j , the principal normals of the boundary curves that meet there must agree modulo sign. Hence, the boundary curves must satisfy the following constraints: corner p00 : n1 (0) = ±n2 (0),
corner p10 : n1 (1) = ±n4 (0),
corner p01 : n2 (1) = ±n3 (0),
corner p11 : n3 (1) = ±n4 (1).
(19)
These conditions imply that, at each corner pi j , the boundary curves meeting there have osculating planes orthogonal to the surface tangent plane Πi j .
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
605
4.4. Condition (C2): global normal orientation constraint Let r(t ) for t ∈ [0, 4] denote the “concatenation” of the four boundary curves r2 (t ), r3 (t ), r4 (t ), r1 (t ), defined as follows: r(t ) = r2 (t ),
t ∈ [0, 1],
r(t ) = r3 (t − 1), t ∈ [1, 2],
r(t ) = r4 (3 − t ), t ∈ [2, 3],
r(t ) = r1 (4 − t ), t ∈ [3, 4].
(20)
Consider the principal normal n(t ) of the concatenated curve r(t ), defined on the interval t ∈ [0, 4]: n(t ) is simply the concatenation of the principal normals ni (t ) of the four boundary curves ri (t ), according to the parameterization (20) used for the definition of r(t ). Let N(u , v ) be the unit normal (4) to the interpolating surface R(u , v ) and let N(t ) for t ∈ [0, 4] be its restriction to the boundary r(t ) of the surface patch. The curve r(t ) : R → R3 and the vectors n(t ) : R → S 2 and N(t ) : R → S 2 (where S 2 is the unit sphere) are regarded as periodic functions, of period 4. Since we desire a regular oriented surface, N(t ) must be a continuous vector function, with N(4) = N(0). Furthermore, since the curves ri (t ) are specified to be geodesics on R(u , v ) we must have N(t ) = ±n(t ) for t ∈ [0, 4]. Now the vector function n(t ) depends only on the prescribed boundary curves ri (t ). If these curves are regular, n(t ) can be discontinuous — i.e., exhibit a sudden reversal — only at their inflection points, or at the parameter values t = 0, 1, 2, 3, 4 identifying the patch corners, where the curves ri (t ) meet. Hence, the existence of a regular oriented surface R(u , v ) that interpolates the concatenated boundary r(t ), with the four individual curves ri (t ) as geodesics, is contingent on the existence of a continuous unit vector function N(t ) such that N(t ) = ±n(t ) for all t ∈ R. Specifically, a regular oriented interpolating surface can exist only if the unit vector function n(t ) exhibits an even number of reversals on the interval t ∈ [0, 4). Remark. We specifically exclude the case of “pathological” inflections — i.e., points where k = 0, and the left and right limits n− and n+ of the principal normal are not parallel or anti-parallel (see Section 2), since no solution can be found if such points are present. Fig. 4 shows an example in which the variation of n(t ) is consistent with the existence of an oriented surface patch with a normal N(t ) that is continuous around the patch boundary, while Fig. 5 shows a case for which the variation of n(t ) is incompatible with such a surface. In the Fig. 4 example, the normal n(t ) exhibits two reversals — one at corner p00 , and the other at an inflection I in the middle of the curve r2 (t ). In this case, a globally continuous solution for N(t ) can be specified by taking N(t ) = −n(t ) for t ∈ (0, 0.5) and N(t ) = n(t ) otherwise. In the example of Fig. 5, there is only one reversal of n(t ), at the corner p00 . Consequently, there is no continuous solution for N(t ). 4.5. Condition (C3): geodesic crossing constraints at corners By Corollary 5, the boundary curves ri (t ) must satisfy the crossing constraint (13) at each corner pi j , with crossing angle A i j defined (see Fig. 6) by:
A 00 = e1 (0), e2 (0)
N(0)
:
sin A 00 = det e1 (0), e2 (0), N(0) ,
cos A 00 = e1 (0), e2 (0) ,
A 01 = e3 (0), e2 (1)
N(1)
:
sin A 01 = det e3 (0), e2 (1), N(1) ,
cos A 01 = e3 (0), e2 (1) ,
A 11 = e3 (1), e4 (1)
N(2)
:
sin A 11 = det e3 (1), e4 (1), N(2) ,
cos A 11 = e3 (1), e4 (1) ,
N(3)
:
sin A 10 = det e1 (1), e4 (0), N(3) ,
cos A 10 = e1 (1), e4 (0) .
A 10 = e1 (1), e4 (0)
(21)
Fig. 4. Left: patch boundaries that satisfy the global normal orientation constraint. Right: continuous surface normal N(t ) along the patch boundary.
606
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
Fig. 5. A set of surface patch boundary curves that are inconsistent with the global normal orientation constraint.
Fig. 6. Geodesic crossing constraints at the four corners of the surface patch. At each corner pi j the surface normal N is directed toward the observer.
Hence, considering scalars
σiL , σi R ∈ {−1, +1} defined by
N(1) = σ2R n2 (1) = σ3L n3 (0),
N(2) = σ3R n3 (1) = σ4R n4 (1),
N(0) = σ1L n1 (0) = σ2L n2 (0),
N(3) = σ1R n1 (1) = σ4L n4 (0),
(22)
the geodesic crossing constraints at the four patch corners become p00 : p01 : p11 : p10 :
σ1L k1 (0) − σ2L k2 (0) cos A 00 + τ1 (0) + τ2 (0) sin A 00 = 0,
σ3L k3 (0) − σ2R k2 (1) cos A 01 + τ3 (0) + τ2 (1) sin A 01 = 0,
σ3R k3 (1) − σ4R k4 (1) cos A 11 + τ3 (1) + τ4 (1) sin A 11 = 0,
σ1R k1 (1) − σ4L k4 (0) cos A 10 + τ1 (1) + τ4 (0) sin A 10 = 0.
(23)
We consider the construction of surface patches with given geodesic boundary curves, knowing a priori that these curves satisfy the above constraints. This stipulation can be met by, for example, selecting the boundary curves as known geodesics on simple analytic surfaces. For general free-form boundary curves, their construction so as to satisfy the system of constraints (23) — or the modification of initial boundary curves so as to satisfy them — is a substantive problem in its own right, which shall be addressed separately in another paper.
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
607
Fig. 7. Coons interpolation of four geodesic boundary curves.
5. Main steps of the construction The algorithm comprises the following steps (see Fig. 7) and is based on Coons interpolation. The approach is similar to that used in Sprynski et al. (2008). • Step 1. Tangent plane. Along each of the four boundary curves ri (t ), the tangent plane Πi (t ) of the surface R(u , v ) is determined by the surface normal vector N(t ) along the patch boundary. Specifically,
Πi (t ) = ri (t ), ei (t ), hi (t ) .
• Step 2. Departure/arrival vectors for Hermite interpolation. A unit vector Ti (t ) in the tangent plane Πi (t ) at each point of ri (t ), and another unit vector Ti +2 (t ) in the tangent plane Πi +2 (t ) at each point of ri +2 (t ), must be specified for i = 1, 2. These vectors specify the departure/arrival directions between the corresponding points on opposite sides ri (t ), ri +2 (t ) of the patch boundary. Consider the situation shown in Fig. 8. A natural choice is to assume Ti (t ) = bi (t ) for i = 1, 3, 4 and Ti (t ) = −bi (t ) for i = 2. However, we shall see in step 4 that the vectors Ti must satisfy one more constraint at the four patch corners. So
608
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
Fig. 8. Vectors along patch boundary.
we adopt a more general form for the departure/arrival vectors, based on the Darboux frame, yielding a larger family of interpolating surfaces. Namely, we set Ti (t ) = cos αi (t )ei (t ) + sin αi (t )hi (t ),
i = 1, 2, 3, 4.
(24)
Here, the angle functions αi (t ) specify the inclination of the vectors Ti relative to the boundary curve tangents ei . Note that we specify unit departure/arrival vectors Ti , since we shall subsequently introduce scalar magnitude functions to control the “parametric speed” of the isoparametric curves on the Coons interpolating surface patch. These functions can be used to manipulate the surface shape for optimum smoothness. • Step 3. Surface construction by bicubic Coons interpolation. Let d1 (u ), d3 (u ) and d2 ( v ), d4 ( v ) be scalar functions, used to modulate the magnitude of the vectors T1 (u ), T3 (u ) and T2 ( v ), T4 ( v ). Consider the three surfaces R13 (u , v ), R24 (u , v ), and R0 (u , v ) defined by
R13 (u , v ) = H 0 ( v )
H 1(v )
R24 (u , v ) = H 0 (u )
H 1 (u )
H 2(v )
H 2 (u )
⎡
⎤
⎡
⎤
r1 (u ) ⎢ d1 (u )T1 (u ) ⎥ ⎥, ⎢ H 3(v ) ⎣ d3 (u )T3 (u ) ⎦ r3 (u ) r2 ( v ) ⎢ d2 ( v )T2 ( v ) ⎥ ⎥, ⎢ H 3 (u ) ⎣ d4 ( v )T4 ( v ) ⎦ r4 ( v )
⎤⎡ ⎡ p ⎤ r 2 (0) r 2 (1) p01 00 H 0(v ) ⎢ r1 (0) Ruv (0, 0) Ruv (0, 1) r3 (0) ⎥ ⎢ H 1 ( v ) ⎥ ⎥ R0 (u , v ) = H 0 (u ) H 1 (u ) H 2 (u ) H 3 (u ) ⎢ ⎣ r (1) Ruv (1, 0) Ruv (1, 1) r (1) ⎦ ⎣ H 2 ( v ) ⎦ . 1 3 H 3(v ) p10 r 4 (0) r 4 (1) p11 In the above expressions for R13 (u , v ), R24 (u , v ), R0 (u , v ), the functions H 0 (t ) = 2t 3 − 3t 2 + 1, H 2 (t ) = t 3 − t 2 ,
H 1 (t ) = t 3 − 2t 2 + t ,
H 3 (t ) = −2t 3 + 3t 2 ,
are the cubic Hermite polynomials (Farin, 2002). The twist vectors Ruv (i , j ) at the four patch corners i , j = 0, 1 will be determined in step 5 below. The desired interpolating surface is then defined (Farin, 2002) by R(u , v ) = R13 (u , v ) + R24 (u , v ) − R0 (u , v ).
• Step 4. Interpolating properties of tangent vectors at corners. Along the boundary v = 0, the equation of the surface R(u , v ) reduces to ⎡ ⎤ r2 (0) − p00 ⎢ d2 (0)T2 (0) − r 1 (0) ⎥ R(u , 0) = r1 (u ) + H 0 (u ) H 1 (u ) H 2 (u ) H 3 (u ) ⎣ ⎦ d4 (0)T4 (0) − r 1 (1) r4 (0) − p10
(25)
and the partial derivatives Ru = ∂ R/∂ u, R v = ∂ R/∂ v along this curve are
⎡
Ru (u , 0) = r 1 (u ) + H 0 (u )
H 1 (u )
H 2 (u )
⎤ r2 (0) − p00 ⎢ d2 (0)T2 (0) − r1 (0) ⎥
H 3 (u ) ⎣
d4 (0)T4 (0) − r 1 (1) r4 (0) − p10
⎦,
(26)
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
⎡
R v (u , 0) = d1 (u )T1 (u ) + H 0 (u ),
H 1 (u ),
H 2 (u ),
609
⎤
r 2 (0) − r 2 (0) ⎢ d ⎢ d2 ( v )T2 ( v )| v =0 − Ruv (0, 0) ⎥ ⎥ H 3 (u ) ⎢ dv ⎥. ⎣ d d4 ( v )T4 ( v )| v =0 − Ruv (1, 0) ⎦ dv
(27)
r 4 (0) − r 4 (0)
From expressions (25) and (26), we deduce that the interpolation constraints R(u , 0) = r1 (u ) and Ru (u , 0) = r 1 (u ) are satisfied if and only if d2 (0)T2 (0) = r 1 (0) and d4 (0)T4 (0) = r 1 (1). Hence, for the four boundary curves, we obtain the following natural constraints: d1 (0)T1 (0) = r 2 (0)
and
d1 (1)T1 (1) = r 4 (0),
(28)
d2 (0)T2 (0) = r1 (0)
and
d2 (1)T2 (1) = r3 (0),
(29)
d3 (0)T3 (0) = r 2 (1)
and
d3 (1)T3 (1) = r 4 (1),
(30)
d4 (0)T4 (0) = r1 (1)
and
d4 (1)T4 (1) = r3 (1).
(31)
In other words, the products di (t )Ti (t ) of the magnitude functions with the corresponding departure/arrival vectors must coincide with the derivatives of the two adjacent curves at the extremities of each boundary curve ri (t ). To study these conditions in greater detail we consider, for example, the first condition in (28). Invoking the definition (24) of vectors Ti (t ), and noting that the magnitude functions di (t ) are required to be positive, we have d1 (0)T1 (0) = r 2 (0)
⇐⇒ ⇐⇒
d1 (0) = r 2 (0)
cos α1 (0)e1 (0) + sin α1 (0)h1 (0) = e2 (0)
d1 (0) = r 2 (0) and
α1 (0) = A 00 .
The relations (28)–(31) thus lead us to introduce the following scalar coefficients diL , di R , Fig. 9):
d1L := d1 (0) = r 2 (0),
d1R := d1 (1) = r 4 (0),
d2L := d2 (0) = r 1 (0),
d2R := d2 (1) = r 3 (0),
d3L := d3 (0) = r 2 (1),
d3R := d3 (1) = r 4 (1),
d4L
d4R
:= d4 (0) = r 1 (1),
α1L := α1 (0) = + A 00 , α2L := α2 (0) = − A 00 , α3L := α3 (0) = + A 01 , α4L := α4 (0) = − A 10 ,
:= d4 (1) = r 3 (1),
α1R := α1 (1) = + A 10 , α2R := α2 (1) = − A 01 , α3R := α3 (1) = + A 11 , α4R := α4 (1) = − A 11 .
αiL , αi R for i = 1, . . . , 4 (see
(32)
(33)
Thus, the angle functions αi (t ) used in (24) to define the unit departure/arrival vectors Ti (t ) and the corresponding magnitude functions di (t ) must satisfy
αi (0) = αiL , αi (1) = αi R and di (0) = diL , di (1) = di R .
Fig. 9. Binormal vectors (blue) and arrival vectors (red) on the curve r4 (t ). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
610
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
Fig. 10. Twist vector constraint at the patch corner p00 .
We shall use the functions
αi (t ) = αiL H 0 (t ) + αi1 H 1 (t ) + αi2 H 2 (t ) + αi R H 3 (t ),
(34)
di (t ) = diL H 0 (t ) + di1 H 1 (t ) + di2 H 2 (t ) + di R H 3 (t ),
(35)
where the parameters αi1 , αi2 , di1 , di2 will be specified in steps 5 and 6. • Step 5. Geodesic crossing property and twist vectors. From Eq. (27), we observe that the cross-derivative property R v (u , 0) = d1 (u )T1 (u ) along the curve r1 (t ) is satisfied if Ruv (0, 0) =
d2 ( v )T2 ( v )
d dv
Ruv (1, 0) =
and v =0
d dv
d4 ( v )T4 ( v )
. v =0
Considering the cross-derivative property along each boundary curve ri (t ), we obtain the following constraints and definitions: d du d dv
d1 (u )T1 (u )
u =0
d2 ( v )T2 ( v )
v =1
= =
d dv d du
d2 ( v )T2 ( v )
d3 (u )T3 (u )
v =0
u =0
=: Ruv (0, 0),
(36)
=: Ruv (0, 1),
(37)
d d d3 (u )T3 (u ) d4 ( v )T4 ( v ) = =: Ruv (1, 1), du dv u =1 v =1 d d d4 ( v )T4 ( v ) d1 (u )T1 (u ) = =: Ruv (1, 0). dv
v =0
du
(38) (39)
u =1
Consider for example the “twist constraint” (36) at the corner p00 . Since the four unit vectors e1 (0) = T2 (0), e2 (0) = T1 (0), h1 (0), h2 (0) are coplanar (see Fig. 10), we can write e1 (0) = cos A 00 e2 (0) − sin A 00 h2 (0), e2 (0) = cos A 00 e1 (0) + sin A 00 h1 (0).
(40)
Hence (noting from Section 4.1 that A 00 = 0, π ) we deduce that h1 (0) =
− cos A 00 e1 (0) + e2 (0) sin A 00
,
h2 (0) =
−e1 (0) + cos A 00 e2 (0) sin A 00
(41)
.
Using definitions (24), (34), (35), (22), relations (32), (33), (41), and relations (12) from Lemma 2, we express each of the derivatives occurring in (36) with respect to the frame (e1 (0), e2 (0), N(0, 0)), with N(0, 0) being the unit normal to the surface R(u , v ) at the corner p00 . d du
d1 (u )T1 (u )
u =0
= d 1 (0)T1 (0) + d1 (0)
d du
cos α1 (u )e1 (u ) + sin α1 (u )h1 (u )
u =0
= d11 e2 (0) + d1L −α1 (0) sin α1 (0)e1 (0) + α1 (0) cos α1 (0)h1 (0) + cos α1 (0)e 1 (0) + sin α1 (0)h 1 (0) = d11 e2 (0) + d1L α11 − sin α1L e1 (0) + cos α1L h1 (0) + d1L r 1 (0) (cos α1L )σ1L k1 (0) + (sin α1L )τ1 (0) N(0, 0) 1 cos A 00 = −d1L α11 e1 (0) + d11 + d1L α11 e2 (0) sin A 00
sin A 00
+ d1L d2L (cos A 00 )σ1L k1 (0) + (sin A 00 )τ1 (0) N(0, 0).
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
In the same way, we obtain d dv
d2 ( v )T2 ( v )
v =0
= d 2 (0)T2 (0) + d2 (0)
d dv
cos A 00
611
cos α2 ( v )e2 ( v ) + sin α2 ( v )h2 ( v )
v =0
1
= d21 − d2L α21 e1 (0) + d2L α21 e2 (0) sin A 00 sin A 00 + d2L d1L (cos A 00 )σ2L k2 (0) − (sin A 00 )τ2 (0) N(0, 0). Thus, in order to satisfy the constraint (36) at corner p00 , and be able to define the twist vector Ruv (0, 0), we must find parameters α11 , d11 , α21 , d21 such that the three following relations are satisfied:
⎧ 1 cos A 00 ⎪ ⎪ −d1L α11 = d21 − d2L α21 , ⎪ ⎪ sin A 00 sin A 00 ⎨ cos A 00
1
⎪ d11 + d1L α11 = d2L α21 , ⎪ ⎪ sin A 00 sin A 00 ⎪ ⎩ (cos A 00 )σ1L k1 (0) + (sin A 00 )τ1 (0) = (cos A 00 )σ2L k2 (0) − (sin A 00 )τ2 (0).
(42)
The third equation in (42) is the geodesic crossing relation (13), and is satisfied through the assumption that the boundary curves obey the relations (23) set up in Section 4.5. Multiplying by sin A 00 , the system (42) is equivalent to the following linear equations in the four unknowns α11 , α21 , d11 , d21 :
α11 d1L − α21 d2L cos A 00 + d21 sin A 00 = 0, α11 d1L cos A 00 − α21 d2L + d11 sin A 00 = 0.
(43)
This system admits solutions. For example, the parameters d11 , d21 can be freely chosen, and the parameters α11 , α21 are then computed from these equations (see examples in Fig. 14 which illustrate the influence of parameters di j ). Hence, we are able to define the twist vector Ruv (0, 0) at the patch corner p00 . In the same way, the twist vectors Ruv (i , j ) at each corner pi j for i , j ∈ {0, 1} can be defined. We can now state the following result. Proposition 7. Given four regular space curves, as specified in Section 4.1, satisfying the following conditions: (C1) the osculating constraints (19); (C2) the global normal orientation constraint set up in Section 4.4; (C3) the crossing constraints (23) at corners pi j ; there exists a regular oriented surface R(u , v ) interpolating these four curves in such a way that these curves are geodesics of the surface. Conversely, if any of the conditions (C1)–(C3) is not satisfied, such an interpolating surface can not be constructed.
• Step 6. Surface smoothing. In this last step we propose to find optimal parameters di j for i = 1, 2, 3, 4 and j = 1, 2 [parameters αi j being deduced from Eqs. (43)] which provide smooth interpolating surfaces. The criteria involving the curvature and the torsion of the isoparametric curves, applied in Sprynski et al. (2008), are not efficient here, due the corner constraints which mainly impose the departure and ending direction of the interpolating curves. So, we propose to minimize the following functionals. Criterion 1. This criterion minimizes the parametric speed variation along each isoparametric curve.
1 1 min
d
dv
di j
0
0
R v (u¯ , v )
2
1 1 dv du¯ + 0
d Ru (u , v¯ ) du
2
du d v¯ .
0
Criterion 2. Minimization of the Dirichlet energy.
1 1 min di j
0
2 2
Ru (u , v ) + R v (u , v ) du dv .
0
Criterion 3. Minimization of the thin plate spline energy.
1 1 min di j
0
Ruu (u , v )2 + 2Ruv (u , v )2 + R v v (u , v )2 du dv .
0
Computed examples are presented and discussed in Section 6.
612
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
Fig. 11. Geodesic on a cone.
Fig. 12. Four geodesic curves are evaluated on a circular cone in (1). The Coons patch interpolating these curves as boundary geodesics is shown superposed on the cone in (2). Finally, the Coons patch is shown on its own — with and without the departure/arrival vectors — in (3) and (4).
6. Examples The following examples illustrate the application of the existence conditions for surfaces interpolating geodesic boundary curves. A companion paper (Farouki et al., 2008) gives a detailed analysis in the case of Bézier boundary curves, in terms of constraints on the control points that ensure satisfaction of these conditions. 6.1. Analytical geodesics from elementary surfaces We consider here elementary surfaces, which admit analytical expressions for their geodesic curves. These geodesic curves naturally satisfy the constraints (C1), (C2), (C3). The interpolation of the four boundary curves as geodesics is accomplished using the procedure described in Section 5, and we compare the constructed surface patch with the original analytic surface. For simplicity, the parameters αi1 , αi2 and di1 , di2 are all set to zero in these examples.
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
613
Fig. 13. Left and center: two views of the Coons patch interpolant to four great-circle boundary segments on the sphere. Right: an approximation of the whole sphere using 6 patches. Each patch interpolates its boundaries in such a way that they are geodesics of that patch. Note that the construction is not symmetric, since each of the six patches is different from the others.
Fig. 14. Six different surfaces interpolating given geodesic boundary curves, corresponding to different values of the parameters di1 , di2 for i = 1, 2, 3, 4. Surface (n) is associated with the parameters di1 = −8(n − 1) and di2 = 8(n − 1) — these are related to αi1 , αi2 through Eqs. (43).
• Circular cone. Consider a circular cone generated by rotating a straight line about the z-axis. If this line has inclination a with the z-axis, we can parameterize the cone by R(u , v ) = (u cos v , u sin v , Au ) for u 0 and 0 v < 2π , where A = cot a. Geodesic curves on this cone may be parameterized in terms of the angle v as r( v ) =
u0
cos[( v − v 0 ) sin a]
cos v sin v A
where (u 0 , v 0 ) are parameter values such that the geodesic is tangent at r( v 0 ) to the circle of latitude of radius u 0 on the cone (see Fig. 11). Fig. 12 illustrates the construction of a surface patch bounded by four geodesic curves on a circular cone. Although these boundary curves are geodesics of the constructed surface patch, the patch is not (in general) an exact subset of the cone. In Fig. 12, portions of the patch beneath the cone are shown hatched, while those above the cone are shown smoothly shaded. • Sphere. Geodesics on the sphere are simply great-circle arcs. Fig. 13 shows a four-sided patch constructed so as to incorporate four great-circle segments as geodesic boundary curves — again, this patch is not (in general) an exact subset of the sphere. Fig. 13 also illustrates a covering of the entire sphere by six patches. 6.2. Influence of parameters The examples in Fig. 14 illustrate the influence on the surface shape of varying the parameters di1 , di2 for i = 1, 2, 3, 4. These examples show the importance of the smoothing procedures in step 6 of the algorithm, since the shape quality is evidently quite sensitive to these parameters, and ad hoc choices for them can yield surfaces of poor quality. Finally, Fig. 15 shows the outcome of the three different smoothing procedures on a surface constructed using the geodesic boundary data in Fig. 4.
614
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
Fig. 15. Coons patches interpolating the geodesic boundary curves shown in Fig. 4, as smoothed according to the criteria 1, 2, and 3 in step 6 of the algorithm (two different views of each smoothed surface are presented).
7. Conclusion The problem of constructing a rectangular Coons surface patch from given boundary curves was investigated, under the stipulation that these boundary curves correspond to geodesic curves of the constructed surface. The existence of solutions to this problem was shown to be contingent on the satisfaction of two types of constraints by the given boundary curves. A global (topological) constraint concerns the compatibility of the variation of the boundary curve principal normals with the normal to an oriented surface. Local constraints, applicable to each of the four patch corners, relate the angles at which these curves meet to their curvatures and torsions. For boundary curve data satisfying the above constraints, the construction of geodesic-bounded tensor-product surface patches can be accomplished using the bicubically-blended Coons interpolation scheme. By way of illustration, examples of the construction of such patches on elementary analytic surfaces were presented. When the boundaries are specified as compatible free-form polynomial or rational Bézier curves, the resulting surface is not (in general) a polynomial or rational patch, because of the need to unitize certain vectors in its construction. In a companion paper (Farouki et al., 2008) we relax the requirement of unit arrival/departure vectors, and thereby identify additional constraints on the control points and/or weights of given Bézier boundary curves, so as to obtain polynomial or rational surface patches for which these boundaries are geodesic curves. Acknowledgements This work has been accomplished during the visit of the third author to the Department of Mechanical and Aeronautical Engineering, University of California, Davis. References Bennis, C., Vézien, J.-M., Iglésias, G., 1991. Piecewise surface flattening for non-distorted texture mapping. Computer Graphics 25 (4), 237–246. Coons, S., 1964. Surfaces for Computer Aided Design. Technical Report, MIT (available as AD 663 504 from National Technical Information Service, Springfield, VA 22161). Coons, S., 1974. Surface patches and B-spline curves. In: Barnhill, R., Riesenfeld, R. (Eds.), Computer Aided Geometric Design. Academic Press. Do Carmo, M.P., 1976. Differential Geometry of Curves and Surfaces. Prentice-Hall. Farin, G., 2002. Curves and Surfaces for CAGD, 5th edition. Academic Press. Farouki, R.T., Szafran, N., Biard, L., 2008. Construction of Bézier surface patches with Bézier curves as geodesic boundaries. Computer-Aided Design, submitted for publication. Gordon, W., 1983. An operator calculus for surface and volume modeling. IEEE Computer Graphics and Applications 3, 18–22. Paluszny, M., 2008. Cubic polynomial patches through geodesics. Computer-Aided Design 40 (1), 56–61. Sánchez-Reyes, J., Dorado, R., 2008. Constrained design of polynomial surfaces from geodesic curves. Computer-Aided Design 40 (1), 49–55. Sprynski, N., Szafran, N., Lacolle, B., Biard, L., 2008. Surface reconstruction via geodesic interpolation. Computer-Aided Design 40 (4), 480–492. Struik, D.J., 1976. Lectures on Classical Differential Geometry. Dover (reprint). Tucker III, C.L., 1997. Forming of advanced composites. In: Gutowski, T.G. (Ed.), Advanced Composites Manufacturing. Wiley, New York, pp. 297–372. Wang, G., Tang, K., Tai, C.H., 2004. Parametric representation of a surface pencil with a common spatial geodesic. Computer-Aided Design 36 (5), 447–459.
Contents lists available at ScienceDirect
Computer Aided Geometric Design www.elsevier.com/locate/cagd
Existence conditions for Coons patches interpolating geodesic boundary curves R.T. Farouki b , N. Szafran a , L. Biard a,∗ a b
Laboratoire Jean Kuntzmann, Université Joseph Fourier – Grenoble, France Department of Mechanical and Aeronautical Engineering, University of California, Davis, CA 95616, USA
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 19 March 2008 Received in revised form 10 October 2008 Accepted 1 January 2009 Available online 7 January 2009 Keywords: Geodesic curves Surface reconstruction Hermite/Coons interpolation Surface smoothing
Given two pairs of regular space curves r1 (u ), r3 (u ) and r2 ( v ), r4 ( v ) that define a curvilinear rectangle, we consider the problem of constructing a C 2 surface patch R(u , v ) for which these four boundary curves correspond to geodesics of the surface. The possibility of constructing such a surface patch is shown to depend on the given boundary curves satisfying two types of consistency constraints. The first constraint is global in nature, and is concerned with compatibility of the variation of the principal normals along the four curves with the normal to an oriented surface. The second constraint is a local differential condition, relating the curvatures and torsions of the curves meeting at each of the four patch corners to the angle between those curves. For curves satisfying these constraints, the surface patch is constructed using a bicubically-blended Coons interpolation process. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Given four parametric space curves r1 (u ), r3 (u ) and r2 ( v ), r4 ( v ) specifying a curvilinear rectangle, the Coons interpolation procedure (Coons, 1964, 1974; Farin, 2002; Gordon, 1983) defines a surface patch R(u , v ) bounded by these four curves, such that R(u , 0) = r1 (u ), R(u , 1) = r3 (u ) and R(0, v ) = r2 ( v ), R(1, v ) = r4 ( v ). The Coons scheme is motivated by the fact that, in numerous surface design or reconstruction contexts, only the boundary curves of a surface patch are specified, and a means of smoothly “filling in” the interior is needed. In addition to the four boundary curves, the bicubically-blended Coons patch requires transverse derivative data along them — i.e., the four vector functions R v (u , 0), R v (u , 1) and Ru (0, v ), Ru (1, v ) must also be specified. In the Coons interpolation scheme, the boundary curves and derivative data (which together specify the tangent plane variation along the patch boundary) are considered to impart sufficient information to ensure the desired surface shape. In this paper, the Coons patch is modified to admit a more fundamental geometrical significance to the specified boundary curves. Namely, these curves are stipulated to be geodesics on the constructed surface. The motivation for this study comes from the problem of constructing analytic computer representations of free-form surfaces from positional/orientational measurements obtained with the Morphosense — a flexible ribbon-like device with embedded microsensors that, when placed on a physical surface, assumes the shape of a geodesic (Sprynski et al., 2008). By placing the Morphosense on the surface at regular intervals along two different directions, it is subdivided into rectangular patches with geodesic boundary curves.
*
Corresponding author. E-mail address: [email protected] (L. Biard).
0167-8396/$ – see front matter doi:10.1016/j.cagd.2009.01.003
© 2009
Elsevier B.V. All rights reserved.
600
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
The construction of surfaces that incorporate one or two given space curves as geodesics has been considered by several authors (Bennis et al., 1991; Sánchez-Reyes and Dorado, 2008; Tucker, 1997; Wang et al., 2004), in the context of applications such as distortion-free mapping of textures onto free-form surfaces; specifying fabric shapes for garment and shoe design; and in the layout of fibers in composite material structures. However, the problem of constructing rectangular surface patches, when all four prescribed boundary curves are required to be geodesics of the resulting surface, does not appear to have been previously studied. This paper is concerned with extending prior studies that address the problem of interpolating two boundary geodesics by a surface “strip” (Paluszny, 2008; Sprynski et al., 2008) to the case where all four boundary curves of a rectangular patch are specified as geodesics of the constructed surface. In particular, we are concerned with identifying constraints on the boundary curves, whose satisfaction constitutes a sufficient-and-necessary condition for the existence of analytic surfaces that interpolate those curves as geodesics. A companion paper (Farouki et al., 2008) will address the construction of (polynomial or rational) Bézier surface patches that interpolate four given Bézier curves, satisfying the existence conditions, as geodesic boundaries. The plan for this paper is as follows. After reviewing basic concepts concerning the geometry of curves on surfaces in Section 2, the ideas of Sprynski et al. (2008) are extended in Section 3 to identify a local condition that must be satisfied by two intersecting space curves if those curves are to be geodesics of an analytic surface. A global compatibility condition, that must be satisfied by the geodesic boundary curves of a rectangular patch to ensure a continuous oriented surface normal, is then identified in Section 4. Section 5 describes a modified Coons algorithm, for boundary data satisfying the conditions of Sections 3 and 4, to construct surface patches with the given boundary curves as geodesics. Finally, Section 6 presents some representative surfaces computed using this algorithm, and Section 7 summarizes our main results and identifies key issues that warrant further investigation. 2. Background on curves and surfaces In the following discussion, all curves and surfaces are considered to be regular and “sufficiently smooth.” A curve is regular if it admits a tangent line at each point, while a surface is regular if it admits a tangent plane at each point. All surfaces are considered to be oriented. A surface is said to be oriented if its unit normal vector (4) is continuous on each closed regular curve on the surface. The inner product of two vectors u, v in R3 is denoted by u, v. Similarly, the plane through a point p in R3 spanned by two linearly-independent vectors u, v is denoted by [p, u, v]. For linearly-independent unit vectors u, v and a unit vector n such that n ⊥ u and n ⊥ v, we denote by (u, v)n the oriented angle between u and v in the sense of n. Precisely, the angle A = (u, v)n is defined (see Fig. 1) by sin A = det(u, v, n),
cos A = u, v.
(1)
The variable s is employed to denote arc length along a space curve. Note that the arc-length parameterization r : s → r(s) of a curve satisfies r (s) = 1 and r (s) ⊥ r (s) for all s. However, in this paper, a general parameterization r : t → r(t ) is often used in the surface construction problem. The parameters of functions may sometimes be omitted when no confusion can arise.
• With each point r(s) of a curve satisfying r (s) = 0, we associate the Serret–Frenet frame (e(s), n(s), b(s)) where e(s) = r (s), n(s) = r (s)/r (s), and b(s) = e(s)× n(s) are, respectively, the unit tangent, principal normal, and binormal vectors of the curve at the point r(s). The arc-length derivative of the Serret–Frenet frame is governed by the relations e( s ) 0 k(s) 0 e( s ) d (2) n(s) = −k(s) 0 τ (s) n(s) , ds b( s ) 0 −τ (s) 0 b( s ) where the curvature k(s) and torsion
k(s) = r (s)
and
τ (s) =
τ (s) of the curve r(s) are defined by
det(r (s), r (s), r (s))
r (s)2
.
Fig. 1. Angle measurement.
(3)
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
601
The osculating plane at each curve point r(s) is spanned by the two vectors e(s), n(s) and does not depend on the curve parameterization. If k(s) = 0 for some s, then r (s) = 0 and the normal vector n(s) and osculating plane are undefined at that point. This condition identifies an inflection of the curve. • On a regular oriented surface (u , v ) → R(u , v ), the unit normal is defined at each point in terms of the partial derivatives Ru = ∂ R/∂ u, R v = ∂ R/∂ v by N( u , v ) =
Ru (u , v ) × R v (u , v ) . Ru (u , v ) × R v (u , v )
(4)
• Consider a curve r(s) = R(u (s), v (s)) on a surface R(u , v ), where s denotes arc length for the space curve r(s), but not necessarily for the plane curve defined by s → (u (s), v (s)). With each point r(s), we associate the Darboux frame (e(s), h(s), N(s)) — where e(s) is the unit tangent vector of the curve, N(s) is the unit normal vector of the surface at the point R(u (s), v (s)) = r(s), and h(s) = N(s) × e(s). The arc-length derivative of the Darboux frame is given by the relations d
ds
e( s ) h(s) N( s )
=
k g (s) kn (s) 0 0 −τ g (s) −k g (s) 0 −kn (s) τ g (s)
e( s ) h(s) , N( s )
(5)
which define the normal curvature kn (s), the geodesic curvature k g (s), and the geodesic torsion curve r(s) as
kn =
de ds
,N ,
kg =
de ds
,h ,
τg =
dN ds
,h .
τ g (s) at each point of the (6)
• A regular curve r(t ) on a surface R(u , v ) is called a geodesic of the surface if its geodesic curvature is identically zero. From (2) and (5), this is equivalent to requiring that N(s) = ±n(s),
h(s) = ∓b(s)
(7)
— i.e., the Frenet and Darboux frames agree modulo signs. Hence, we have the following useful characterizations of geodesic curves. A regular curve t → r(t ) is a geodesic on the surface R(u , v ) if and only if (D1) the geodesic curvature of r(t ) is identically zero; (D2) the principal normal at each non-inflection point of r(t ) is orthogonal to the surface tangent plane at the point R(u (t ), v (t )) = r(t ); (D3) the osculating plane at each non-inflection point of r(t ) is orthogonal to the surface tangent plane at the point R(u (t ), v (t )) = r(t ). In the case of an isolated inflection point of the curve, we may encounter cases (Do Carmo, 1976) where the osculating plane cannot be defined by continuity, since the limit planes approaching it from the left and the right disagree. For nonisolated inflection points, the vanishing of the curvature over an interval in the parameter implies that the corresponding curve segment degenerates to a straight line (and hence a geodesic) on the surface. For such a segment, the Darboux frame can be defined at each point, but the Frenet frame is undefined. 3. Geodesic curves crossing on a regular surface The normal curvature kn and geodesic torsion τ g at any point of a curve r(t ) on a regular C 2 oriented surface R(u , v ) depend only on the tangent direction of r(t ) at that point. Namely, if Πp is the tangent plane at a given point p of the surface R(u , v ), and L is any line in Πp passing through p, all curves on R(u , v ) that go through p, and have L as their tangent line there, possess the same normal curvature kn and geodesic torsion τ g . This property allows us to define the normal curvature K n ( L ) and geodesic torsion T g ( L ) at each point p of a surface R(u , v ), in the direction of any line L in the tangent plane Πp . At a fixed surface point p, the normal curvature K n ( L ) exhibits two extrema K 1 , K 2 (the principal curvatures), corresponding to two orthogonal directions L 1 , L 2 (the principal directions). Consider an orthonormal basis (e1 , e2 ) for Πp aligned with the principal directions L 1 , L 2 such that e1 × e2 defines the unit surface normal N at p. Then each line L in the tangent plane Πp corresponds to an angle α (with 0 α < π ), such that the normal curvature of the surface at point p can be regarded as a π -periodic function of the angle α . It can be shown (Do Carmo, 1976) that K n (α ) = K 1 cos2 α + K 2 sin2 α . Similarly, the geodesic torsion of the surface at point p can be regarded as a given (Do Carmo, 1976) by
T g (α ) = ( K 1 − K 2 ) sin α cos α .
(8)
π -periodic function T g (α ) of the angle α , (9)
602
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
Remark. A point p at which K 1 = K 2 is called an umbilic of the surface. At such a point, the normal curvature K n (α ) defined by (8) is the same in all directions α (so the surface shape is locally like a sphere), while the geodesic torsion (9) vanishes identically. For subsequent use, we now derive Proposition 3, which identifies a necessary condition for two intersecting space curves to be geodesics on a regular surface. First, notice that we have only three degrees of freedom in the curvature distribution, namely the two principal curvatures K 1 , K 2 and the orientation of the principal directions. Thus, the normal curvature and geodesic torsion are locally determined in any direction from relations (8) and (9). So, the normal curvatures and geodesic torsions along two directions in a point of a regular surface must satisfy one compatibility equation, which is established in the following lemma. Lemma 1. For two directions in the surface tangent plane Πp specified by the angles torsion satisfy the relation
α1 , α2 the normal curvature and the geodesic
sin(α2 − α1 ) T g (α1 ) + T g (α2 ) = cos(α2 − α1 ) K n (α1 ) − K n (α2 ) .
(10)
Proof. By direct computation, we have
sin(α2 − α1 ) T g (α1 ) + T g (α2 )
= ( K 1 − K 2 )(sin α2 cos α1 − sin α1 cos α2 )(sin α1 cos α1 + sin α2 cos α2 )
= ( K 1 − K 2 ) sin α1 sin α2 cos2 α1 − cos2 α2 + cos α1 cos α2 sin2 α2 − sin2 α1 1 = ( K 1 − K 2 ) sin α1 sin α2 (cos 2α1 − cos 2α2 ) + cos α1 cos α2 (cos 2α1 − cos 2α2 ) = =
2 1 2 1 2
( K 1 − K 2 ) cos(α2 − α1 )(cos 2α1 − cos 2α2 )
cos(α2 − α1 ) K 1 2 cos2 α1 − 1 − 2 cos2 α2 + 1 − K 2 1 − 2 sin2 α1 − 1 + 2 sin2 α2
= cos(α2 − α1 ) K n (α1 ) − K n (α2 ) .
2
We now consider in greater detail the relations (7) that characterize a geodesic curve on a surface. Lemma 2. Consider a regular curve r(s) on a surface R(u , v ) with notations as introduced in Section 2. Then, setting N(s) = σ n(s) with σ = ±1, we have h(s) = −σ b(s),
k g (s) = 0,
τ g (s) = −τ (s),
kn (s) = σ k(s).
(11)
Proof. Setting N(s) = σ n(s) so that h(s) = −σ b(s), a direct calculation using relations (2) gives d ds
e( s ) h(s) N( s )
e( s ) 0 k(s) 0 = σ τ (s) 0 0 −σ b(s) = ds σ n(s) −σ k(s) 0 σ τ (s) 0 0 σ k ( s ) e( s ) = 0 0 τ (s) h(s) , 0 N( s ) −σ k(s) −τ (s) d
e( s ) n(s) b( s )
which corresponds to the stated results by comparison with relations (5).
(12)
2
Proposition 3. Consider two geodesics r1 (s) and r2 (s) parameterized by arc length on the surface R(u , v ), with principal normal n1 (s) and n2 (s), crossing at the point p = r1 (s1 ) = r2 (s2 ). Assuming that p is not an inflection on either of the curves r1 (s) and r2 (s), let α = (r 1 (s1 ), r 2 (s2 ))Np be the oriented angle between them at p, in the sense of the surface normal Np at that point. Also, let ki (s) and τi (s) be the curvature and torsion of ri (s) for i = 1, 2. Then, for the values σ1 , σ2 ∈ {−1, +1} such that Np = σ1 n1 (s1 ) = σ2 n2 (s2 ), we have
τ1 (s1 ) + τ2 (s2 ) sin α = σ2k2 (s2 ) − σ1k1 (s1 ) cos α .
(13)
Proof. Consider first any geodesic r(s), free of inflections, parameterized by arc length on the surface R(u , v ). Its geodesic curvature satisfies k g (s) ≡ 0, and its principal normal n(s) agrees (modulo sign) with the surface normal N(s) (here N(s) is the restriction of the unit surface normal to the curve r(s)). Setting N(s) = σ n(s) with σ = ±1, Lemma 2 gives
τ g (s) = −τ (s) and kn (s) = σ k(s). Hence, if the tangents to the two geodesics r1 (s) and r2 (s) at the point p are identified by angles the result follows directly from Lemma 1, with α = α2 − α1 . 2
(14)
α1 and α2 in the plane Πp ,
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
603
Fig. 2. Geodesics on a cylinder.
Proposition 4 (Consistency with surface orientation). The geodesic crossing relation (13) in Proposition 3 does not depend on the surface orientation. Proof. Consider a new parameterization of the surface in Proposition 3, that induces at each point a surface normal opposite to the original normal. Then the angle α and the values σ1 , σ2 will change sign, but the curvatures and the torsions, which are intrinsic to the curves r1 (s) and r2 (s), remain unchanged. Hence, we obtain
τ1 (s1 ) + τ2 (s2 ) sin(−α ) = (−σ2 )k2 (s2 ) − (−σ1 )k1 (s1 ) cos(−α ),
which is identical to (13).
(15)
2
Corollary 5. If two space curves intersecting at a point p do not satisfy the relation (13), no regular C 2 oriented surface can interpolate these curves in such a manner that they are geodesics of the surface. Example. Consider the cylinder x2 + y 2 = R 2 , and for a1 , a2 ∈ R (with a1 = a2 ) the geodesic curves (circular helices) defined by ri (t ) = ( R cos t , R sin t , ai t ) on it, crossing at p = r1 (0) = r2 (0) = ( R , 0, 0) with the derivatives r i (0) = (0, R , ai ) at that point. Let αi be the angle that r i (0) makes with the horizontal plane. Note that, for each circular helix, the principal normal n(t ) = (− cos t , − sin t , 0) points to the inside of the cylinder. From Proposition 4, we can choose an orientation of the cylinder such that σ1 = σ2 = 1, so that positive angles αi must be counted clockwise in Fig. 2. Thus, the positive angle α between the two vectors r i (0) is α = α1 − α2 , and we have tan α = tan(α1 − α2 ) =
a1 / R − a2 / R 1 + (a1 / R )(a2 / R )
=
R (a1 − a2 ) R 2 + a1 a2
For each helical curve ri (t ), the curvature ki (t ) and torsion direct computation gives k2 (0) − k1 (0)
τ1 (0) + τ2 (0)
=
R (a1 − a2 ) R 2 + a1 a2
= tan α =
sin α cos α
.
τi (t ) are constant, ki = R /(a2i + R 2 ) and τi = ai /(a2i + R 2 ), and a
,
which is relation (13). As another simple example, one can easily check that this relation is satisfied by great-circle geodesics crossing on a sphere. Proposition 6. (See Struik, 1976.) Through each point of a smooth surface, there exists a geodesic in every direction, and a geodesic is uniquely determined by an initial point and tangent direction at that point. Henceforth we shall consider only distinct crossing geodesics, with a non-zero angle of intersection. For such curves, the tangent plane Πp of the surface R(u , v ) at the point p is spanned by the two derivatives r 1 (s1 ), r 2 (s2 ). 4. Geodesic interpolation problem We now consider the problem of defining a smooth rectangular surface patch, bounded by four given parametric curves, in such a manner that these border curves are geodesics of the surface. This problem, introduced in Section 4.1, is called the “geodesic interpolation problem.” After introducing notations in Section 4.2 and identifying a local corner compatibility condition on the curve principal normals in Section 4.3, we introduce in Section 4.4 the global normal compatibility constraint, which ensures that the surface normal vector varies continuously around the patch boundary. Finally, Section 4.5 specifies the local differential constraints that must be satisfied by the border curves meeting at each patch corner, in order to ensure the existence of a surface for which these curves are geodesics. The details of the surface construction procedure will be presented in Section 5.
604
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
Fig. 3. Surface patch boundary curves.
4.1. Specification of the problem Consider, as illustrated in Fig. 3, four regular space curves r1 (u ), r3 (u ) and r2 ( v ), r4 ( v ) with u , v ∈ [0, 1] such that r1 (0) = r2 (0) = p00 , r1 (1) = r4 (0) = p10 , r2 (1) = r3 (0) = p01 , r3 (1) = r4 (1) = p11 .
(16)
As noted earlier, these curves are assumed to be “sufficiently smooth.” Also, consistent with the last remark in Section 3, their derivatives at each corner pi j , i , j = 0, 1, are assumed to be linearly independent, in order to define the tangent plane Πi j of the interpolating surface R(u , v ) (see below) at pi j . Our goal is to construct a surface that interpolates these four curves in such a manner that they are geodesics on the constructed surface. More precisely, if R(u , v ) for (u , v ) ∈ [0, 1]2 denotes the interpolating surface, it must exhibit the following three properties.
• South-North interpolating property: R(u , 0) = r1 (u ) and R(u , 1) = r3 (u ) for u ∈ [0, 1]; • West-East interpolating property: R(0, v ) = r2 ( v ) and R(1, v ) = r4 ( v ) for v ∈ [0, 1]; • Geodesic property: The principal normal at each non-inflectional point of the curves r1 (u ), r3 (u ) and r2 ( v ), r4 ( v ) is normal to the surface R(u , v ). Note that each of the geodesic curves r1 (u ), r3 (u ) and r2 ( v ), r4 ( v ) is also an isoparametric curve of the surface R(u , v ). 4.2. Notations Let t be either of the parametric variables u or v. For i = 1, 2, 3, 4 we denote by (ei (t ), ni (t ), bi (t )) and ki (t ), Serret–Frenet frame and the curvature and torsion of the curves ri (t ) at each non-inflection point. Namely, ei (t ) = ki (t ) =
r i (t )
||r i (t )||
,
bi (t ) =
r i (t ) × r i (t ) , r i (t )3
r i (t ) × r i (t )
||r i (t ) × r i (t )||
τi (t ) =
,
ni (t ) = bi (t ) × ei (t ),
det(r i (t ), r i (t ), r (t )) i
r i (t ) × r i (t )2
.
τi (t ) the
(17) (18)
At inflection points, the curvature is zero and the torsion is undefined. Note that, whereas the principal normal vector is intrinsic to the curve, the tangent and binormal vectors depend on the sense of the parameterization. 4.3. Condition (C1): osculating constraints at corners From the geodesic definitions in Section 2 we see that, at each corner pi j , the principal normals of the boundary curves that meet there must agree modulo sign. Hence, the boundary curves must satisfy the following constraints: corner p00 : n1 (0) = ±n2 (0),
corner p10 : n1 (1) = ±n4 (0),
corner p01 : n2 (1) = ±n3 (0),
corner p11 : n3 (1) = ±n4 (1).
(19)
These conditions imply that, at each corner pi j , the boundary curves meeting there have osculating planes orthogonal to the surface tangent plane Πi j .
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
605
4.4. Condition (C2): global normal orientation constraint Let r(t ) for t ∈ [0, 4] denote the “concatenation” of the four boundary curves r2 (t ), r3 (t ), r4 (t ), r1 (t ), defined as follows: r(t ) = r2 (t ),
t ∈ [0, 1],
r(t ) = r3 (t − 1), t ∈ [1, 2],
r(t ) = r4 (3 − t ), t ∈ [2, 3],
r(t ) = r1 (4 − t ), t ∈ [3, 4].
(20)
Consider the principal normal n(t ) of the concatenated curve r(t ), defined on the interval t ∈ [0, 4]: n(t ) is simply the concatenation of the principal normals ni (t ) of the four boundary curves ri (t ), according to the parameterization (20) used for the definition of r(t ). Let N(u , v ) be the unit normal (4) to the interpolating surface R(u , v ) and let N(t ) for t ∈ [0, 4] be its restriction to the boundary r(t ) of the surface patch. The curve r(t ) : R → R3 and the vectors n(t ) : R → S 2 and N(t ) : R → S 2 (where S 2 is the unit sphere) are regarded as periodic functions, of period 4. Since we desire a regular oriented surface, N(t ) must be a continuous vector function, with N(4) = N(0). Furthermore, since the curves ri (t ) are specified to be geodesics on R(u , v ) we must have N(t ) = ±n(t ) for t ∈ [0, 4]. Now the vector function n(t ) depends only on the prescribed boundary curves ri (t ). If these curves are regular, n(t ) can be discontinuous — i.e., exhibit a sudden reversal — only at their inflection points, or at the parameter values t = 0, 1, 2, 3, 4 identifying the patch corners, where the curves ri (t ) meet. Hence, the existence of a regular oriented surface R(u , v ) that interpolates the concatenated boundary r(t ), with the four individual curves ri (t ) as geodesics, is contingent on the existence of a continuous unit vector function N(t ) such that N(t ) = ±n(t ) for all t ∈ R. Specifically, a regular oriented interpolating surface can exist only if the unit vector function n(t ) exhibits an even number of reversals on the interval t ∈ [0, 4). Remark. We specifically exclude the case of “pathological” inflections — i.e., points where k = 0, and the left and right limits n− and n+ of the principal normal are not parallel or anti-parallel (see Section 2), since no solution can be found if such points are present. Fig. 4 shows an example in which the variation of n(t ) is consistent with the existence of an oriented surface patch with a normal N(t ) that is continuous around the patch boundary, while Fig. 5 shows a case for which the variation of n(t ) is incompatible with such a surface. In the Fig. 4 example, the normal n(t ) exhibits two reversals — one at corner p00 , and the other at an inflection I in the middle of the curve r2 (t ). In this case, a globally continuous solution for N(t ) can be specified by taking N(t ) = −n(t ) for t ∈ (0, 0.5) and N(t ) = n(t ) otherwise. In the example of Fig. 5, there is only one reversal of n(t ), at the corner p00 . Consequently, there is no continuous solution for N(t ). 4.5. Condition (C3): geodesic crossing constraints at corners By Corollary 5, the boundary curves ri (t ) must satisfy the crossing constraint (13) at each corner pi j , with crossing angle A i j defined (see Fig. 6) by:
A 00 = e1 (0), e2 (0)
N(0)
:
sin A 00 = det e1 (0), e2 (0), N(0) ,
cos A 00 = e1 (0), e2 (0) ,
A 01 = e3 (0), e2 (1)
N(1)
:
sin A 01 = det e3 (0), e2 (1), N(1) ,
cos A 01 = e3 (0), e2 (1) ,
A 11 = e3 (1), e4 (1)
N(2)
:
sin A 11 = det e3 (1), e4 (1), N(2) ,
cos A 11 = e3 (1), e4 (1) ,
N(3)
:
sin A 10 = det e1 (1), e4 (0), N(3) ,
cos A 10 = e1 (1), e4 (0) .
A 10 = e1 (1), e4 (0)
(21)
Fig. 4. Left: patch boundaries that satisfy the global normal orientation constraint. Right: continuous surface normal N(t ) along the patch boundary.
606
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
Fig. 5. A set of surface patch boundary curves that are inconsistent with the global normal orientation constraint.
Fig. 6. Geodesic crossing constraints at the four corners of the surface patch. At each corner pi j the surface normal N is directed toward the observer.
Hence, considering scalars
σiL , σi R ∈ {−1, +1} defined by
N(1) = σ2R n2 (1) = σ3L n3 (0),
N(2) = σ3R n3 (1) = σ4R n4 (1),
N(0) = σ1L n1 (0) = σ2L n2 (0),
N(3) = σ1R n1 (1) = σ4L n4 (0),
(22)
the geodesic crossing constraints at the four patch corners become p00 : p01 : p11 : p10 :
σ1L k1 (0) − σ2L k2 (0) cos A 00 + τ1 (0) + τ2 (0) sin A 00 = 0,
σ3L k3 (0) − σ2R k2 (1) cos A 01 + τ3 (0) + τ2 (1) sin A 01 = 0,
σ3R k3 (1) − σ4R k4 (1) cos A 11 + τ3 (1) + τ4 (1) sin A 11 = 0,
σ1R k1 (1) − σ4L k4 (0) cos A 10 + τ1 (1) + τ4 (0) sin A 10 = 0.
(23)
We consider the construction of surface patches with given geodesic boundary curves, knowing a priori that these curves satisfy the above constraints. This stipulation can be met by, for example, selecting the boundary curves as known geodesics on simple analytic surfaces. For general free-form boundary curves, their construction so as to satisfy the system of constraints (23) — or the modification of initial boundary curves so as to satisfy them — is a substantive problem in its own right, which shall be addressed separately in another paper.
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
607
Fig. 7. Coons interpolation of four geodesic boundary curves.
5. Main steps of the construction The algorithm comprises the following steps (see Fig. 7) and is based on Coons interpolation. The approach is similar to that used in Sprynski et al. (2008). • Step 1. Tangent plane. Along each of the four boundary curves ri (t ), the tangent plane Πi (t ) of the surface R(u , v ) is determined by the surface normal vector N(t ) along the patch boundary. Specifically,
Πi (t ) = ri (t ), ei (t ), hi (t ) .
• Step 2. Departure/arrival vectors for Hermite interpolation. A unit vector Ti (t ) in the tangent plane Πi (t ) at each point of ri (t ), and another unit vector Ti +2 (t ) in the tangent plane Πi +2 (t ) at each point of ri +2 (t ), must be specified for i = 1, 2. These vectors specify the departure/arrival directions between the corresponding points on opposite sides ri (t ), ri +2 (t ) of the patch boundary. Consider the situation shown in Fig. 8. A natural choice is to assume Ti (t ) = bi (t ) for i = 1, 3, 4 and Ti (t ) = −bi (t ) for i = 2. However, we shall see in step 4 that the vectors Ti must satisfy one more constraint at the four patch corners. So
608
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
Fig. 8. Vectors along patch boundary.
we adopt a more general form for the departure/arrival vectors, based on the Darboux frame, yielding a larger family of interpolating surfaces. Namely, we set Ti (t ) = cos αi (t )ei (t ) + sin αi (t )hi (t ),
i = 1, 2, 3, 4.
(24)
Here, the angle functions αi (t ) specify the inclination of the vectors Ti relative to the boundary curve tangents ei . Note that we specify unit departure/arrival vectors Ti , since we shall subsequently introduce scalar magnitude functions to control the “parametric speed” of the isoparametric curves on the Coons interpolating surface patch. These functions can be used to manipulate the surface shape for optimum smoothness. • Step 3. Surface construction by bicubic Coons interpolation. Let d1 (u ), d3 (u ) and d2 ( v ), d4 ( v ) be scalar functions, used to modulate the magnitude of the vectors T1 (u ), T3 (u ) and T2 ( v ), T4 ( v ). Consider the three surfaces R13 (u , v ), R24 (u , v ), and R0 (u , v ) defined by
R13 (u , v ) = H 0 ( v )
H 1(v )
R24 (u , v ) = H 0 (u )
H 1 (u )
H 2(v )
H 2 (u )
⎡
⎤
⎡
⎤
r1 (u ) ⎢ d1 (u )T1 (u ) ⎥ ⎥, ⎢ H 3(v ) ⎣ d3 (u )T3 (u ) ⎦ r3 (u ) r2 ( v ) ⎢ d2 ( v )T2 ( v ) ⎥ ⎥, ⎢ H 3 (u ) ⎣ d4 ( v )T4 ( v ) ⎦ r4 ( v )
⎤⎡ ⎡ p ⎤ r 2 (0) r 2 (1) p01 00 H 0(v ) ⎢ r1 (0) Ruv (0, 0) Ruv (0, 1) r3 (0) ⎥ ⎢ H 1 ( v ) ⎥ ⎥ R0 (u , v ) = H 0 (u ) H 1 (u ) H 2 (u ) H 3 (u ) ⎢ ⎣ r (1) Ruv (1, 0) Ruv (1, 1) r (1) ⎦ ⎣ H 2 ( v ) ⎦ . 1 3 H 3(v ) p10 r 4 (0) r 4 (1) p11 In the above expressions for R13 (u , v ), R24 (u , v ), R0 (u , v ), the functions H 0 (t ) = 2t 3 − 3t 2 + 1, H 2 (t ) = t 3 − t 2 ,
H 1 (t ) = t 3 − 2t 2 + t ,
H 3 (t ) = −2t 3 + 3t 2 ,
are the cubic Hermite polynomials (Farin, 2002). The twist vectors Ruv (i , j ) at the four patch corners i , j = 0, 1 will be determined in step 5 below. The desired interpolating surface is then defined (Farin, 2002) by R(u , v ) = R13 (u , v ) + R24 (u , v ) − R0 (u , v ).
• Step 4. Interpolating properties of tangent vectors at corners. Along the boundary v = 0, the equation of the surface R(u , v ) reduces to ⎡ ⎤ r2 (0) − p00 ⎢ d2 (0)T2 (0) − r 1 (0) ⎥ R(u , 0) = r1 (u ) + H 0 (u ) H 1 (u ) H 2 (u ) H 3 (u ) ⎣ ⎦ d4 (0)T4 (0) − r 1 (1) r4 (0) − p10
(25)
and the partial derivatives Ru = ∂ R/∂ u, R v = ∂ R/∂ v along this curve are
⎡
Ru (u , 0) = r 1 (u ) + H 0 (u )
H 1 (u )
H 2 (u )
⎤ r2 (0) − p00 ⎢ d2 (0)T2 (0) − r1 (0) ⎥
H 3 (u ) ⎣
d4 (0)T4 (0) − r 1 (1) r4 (0) − p10
⎦,
(26)
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
⎡
R v (u , 0) = d1 (u )T1 (u ) + H 0 (u ),
H 1 (u ),
H 2 (u ),
609
⎤
r 2 (0) − r 2 (0) ⎢ d ⎢ d2 ( v )T2 ( v )| v =0 − Ruv (0, 0) ⎥ ⎥ H 3 (u ) ⎢ dv ⎥. ⎣ d d4 ( v )T4 ( v )| v =0 − Ruv (1, 0) ⎦ dv
(27)
r 4 (0) − r 4 (0)
From expressions (25) and (26), we deduce that the interpolation constraints R(u , 0) = r1 (u ) and Ru (u , 0) = r 1 (u ) are satisfied if and only if d2 (0)T2 (0) = r 1 (0) and d4 (0)T4 (0) = r 1 (1). Hence, for the four boundary curves, we obtain the following natural constraints: d1 (0)T1 (0) = r 2 (0)
and
d1 (1)T1 (1) = r 4 (0),
(28)
d2 (0)T2 (0) = r1 (0)
and
d2 (1)T2 (1) = r3 (0),
(29)
d3 (0)T3 (0) = r 2 (1)
and
d3 (1)T3 (1) = r 4 (1),
(30)
d4 (0)T4 (0) = r1 (1)
and
d4 (1)T4 (1) = r3 (1).
(31)
In other words, the products di (t )Ti (t ) of the magnitude functions with the corresponding departure/arrival vectors must coincide with the derivatives of the two adjacent curves at the extremities of each boundary curve ri (t ). To study these conditions in greater detail we consider, for example, the first condition in (28). Invoking the definition (24) of vectors Ti (t ), and noting that the magnitude functions di (t ) are required to be positive, we have d1 (0)T1 (0) = r 2 (0)
⇐⇒ ⇐⇒
d1 (0) = r 2 (0)
cos α1 (0)e1 (0) + sin α1 (0)h1 (0) = e2 (0)
d1 (0) = r 2 (0) and
α1 (0) = A 00 .
The relations (28)–(31) thus lead us to introduce the following scalar coefficients diL , di R , Fig. 9):
d1L := d1 (0) = r 2 (0),
d1R := d1 (1) = r 4 (0),
d2L := d2 (0) = r 1 (0),
d2R := d2 (1) = r 3 (0),
d3L := d3 (0) = r 2 (1),
d3R := d3 (1) = r 4 (1),
d4L
d4R
:= d4 (0) = r 1 (1),
α1L := α1 (0) = + A 00 , α2L := α2 (0) = − A 00 , α3L := α3 (0) = + A 01 , α4L := α4 (0) = − A 10 ,
:= d4 (1) = r 3 (1),
α1R := α1 (1) = + A 10 , α2R := α2 (1) = − A 01 , α3R := α3 (1) = + A 11 , α4R := α4 (1) = − A 11 .
αiL , αi R for i = 1, . . . , 4 (see
(32)
(33)
Thus, the angle functions αi (t ) used in (24) to define the unit departure/arrival vectors Ti (t ) and the corresponding magnitude functions di (t ) must satisfy
αi (0) = αiL , αi (1) = αi R and di (0) = diL , di (1) = di R .
Fig. 9. Binormal vectors (blue) and arrival vectors (red) on the curve r4 (t ). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
610
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
Fig. 10. Twist vector constraint at the patch corner p00 .
We shall use the functions
αi (t ) = αiL H 0 (t ) + αi1 H 1 (t ) + αi2 H 2 (t ) + αi R H 3 (t ),
(34)
di (t ) = diL H 0 (t ) + di1 H 1 (t ) + di2 H 2 (t ) + di R H 3 (t ),
(35)
where the parameters αi1 , αi2 , di1 , di2 will be specified in steps 5 and 6. • Step 5. Geodesic crossing property and twist vectors. From Eq. (27), we observe that the cross-derivative property R v (u , 0) = d1 (u )T1 (u ) along the curve r1 (t ) is satisfied if Ruv (0, 0) =
d2 ( v )T2 ( v )
d dv
Ruv (1, 0) =
and v =0
d dv
d4 ( v )T4 ( v )
. v =0
Considering the cross-derivative property along each boundary curve ri (t ), we obtain the following constraints and definitions: d du d dv
d1 (u )T1 (u )
u =0
d2 ( v )T2 ( v )
v =1
= =
d dv d du
d2 ( v )T2 ( v )
d3 (u )T3 (u )
v =0
u =0
=: Ruv (0, 0),
(36)
=: Ruv (0, 1),
(37)
d d d3 (u )T3 (u ) d4 ( v )T4 ( v ) = =: Ruv (1, 1), du dv u =1 v =1 d d d4 ( v )T4 ( v ) d1 (u )T1 (u ) = =: Ruv (1, 0). dv
v =0
du
(38) (39)
u =1
Consider for example the “twist constraint” (36) at the corner p00 . Since the four unit vectors e1 (0) = T2 (0), e2 (0) = T1 (0), h1 (0), h2 (0) are coplanar (see Fig. 10), we can write e1 (0) = cos A 00 e2 (0) − sin A 00 h2 (0), e2 (0) = cos A 00 e1 (0) + sin A 00 h1 (0).
(40)
Hence (noting from Section 4.1 that A 00 = 0, π ) we deduce that h1 (0) =
− cos A 00 e1 (0) + e2 (0) sin A 00
,
h2 (0) =
−e1 (0) + cos A 00 e2 (0) sin A 00
(41)
.
Using definitions (24), (34), (35), (22), relations (32), (33), (41), and relations (12) from Lemma 2, we express each of the derivatives occurring in (36) with respect to the frame (e1 (0), e2 (0), N(0, 0)), with N(0, 0) being the unit normal to the surface R(u , v ) at the corner p00 . d du
d1 (u )T1 (u )
u =0
= d 1 (0)T1 (0) + d1 (0)
d du
cos α1 (u )e1 (u ) + sin α1 (u )h1 (u )
u =0
= d11 e2 (0) + d1L −α1 (0) sin α1 (0)e1 (0) + α1 (0) cos α1 (0)h1 (0) + cos α1 (0)e 1 (0) + sin α1 (0)h 1 (0) = d11 e2 (0) + d1L α11 − sin α1L e1 (0) + cos α1L h1 (0) + d1L r 1 (0) (cos α1L )σ1L k1 (0) + (sin α1L )τ1 (0) N(0, 0) 1 cos A 00 = −d1L α11 e1 (0) + d11 + d1L α11 e2 (0) sin A 00
sin A 00
+ d1L d2L (cos A 00 )σ1L k1 (0) + (sin A 00 )τ1 (0) N(0, 0).
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
In the same way, we obtain d dv
d2 ( v )T2 ( v )
v =0
= d 2 (0)T2 (0) + d2 (0)
d dv
cos A 00
611
cos α2 ( v )e2 ( v ) + sin α2 ( v )h2 ( v )
v =0
1
= d21 − d2L α21 e1 (0) + d2L α21 e2 (0) sin A 00 sin A 00 + d2L d1L (cos A 00 )σ2L k2 (0) − (sin A 00 )τ2 (0) N(0, 0). Thus, in order to satisfy the constraint (36) at corner p00 , and be able to define the twist vector Ruv (0, 0), we must find parameters α11 , d11 , α21 , d21 such that the three following relations are satisfied:
⎧ 1 cos A 00 ⎪ ⎪ −d1L α11 = d21 − d2L α21 , ⎪ ⎪ sin A 00 sin A 00 ⎨ cos A 00
1
⎪ d11 + d1L α11 = d2L α21 , ⎪ ⎪ sin A 00 sin A 00 ⎪ ⎩ (cos A 00 )σ1L k1 (0) + (sin A 00 )τ1 (0) = (cos A 00 )σ2L k2 (0) − (sin A 00 )τ2 (0).
(42)
The third equation in (42) is the geodesic crossing relation (13), and is satisfied through the assumption that the boundary curves obey the relations (23) set up in Section 4.5. Multiplying by sin A 00 , the system (42) is equivalent to the following linear equations in the four unknowns α11 , α21 , d11 , d21 :
α11 d1L − α21 d2L cos A 00 + d21 sin A 00 = 0, α11 d1L cos A 00 − α21 d2L + d11 sin A 00 = 0.
(43)
This system admits solutions. For example, the parameters d11 , d21 can be freely chosen, and the parameters α11 , α21 are then computed from these equations (see examples in Fig. 14 which illustrate the influence of parameters di j ). Hence, we are able to define the twist vector Ruv (0, 0) at the patch corner p00 . In the same way, the twist vectors Ruv (i , j ) at each corner pi j for i , j ∈ {0, 1} can be defined. We can now state the following result. Proposition 7. Given four regular space curves, as specified in Section 4.1, satisfying the following conditions: (C1) the osculating constraints (19); (C2) the global normal orientation constraint set up in Section 4.4; (C3) the crossing constraints (23) at corners pi j ; there exists a regular oriented surface R(u , v ) interpolating these four curves in such a way that these curves are geodesics of the surface. Conversely, if any of the conditions (C1)–(C3) is not satisfied, such an interpolating surface can not be constructed.
• Step 6. Surface smoothing. In this last step we propose to find optimal parameters di j for i = 1, 2, 3, 4 and j = 1, 2 [parameters αi j being deduced from Eqs. (43)] which provide smooth interpolating surfaces. The criteria involving the curvature and the torsion of the isoparametric curves, applied in Sprynski et al. (2008), are not efficient here, due the corner constraints which mainly impose the departure and ending direction of the interpolating curves. So, we propose to minimize the following functionals. Criterion 1. This criterion minimizes the parametric speed variation along each isoparametric curve.
1 1 min
d
dv
di j
0
0
R v (u¯ , v )
2
1 1 dv du¯ + 0
d Ru (u , v¯ ) du
2
du d v¯ .
0
Criterion 2. Minimization of the Dirichlet energy.
1 1 min di j
0
2 2
Ru (u , v ) + R v (u , v ) du dv .
0
Criterion 3. Minimization of the thin plate spline energy.
1 1 min di j
0
Ruu (u , v )2 + 2Ruv (u , v )2 + R v v (u , v )2 du dv .
0
Computed examples are presented and discussed in Section 6.
612
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
Fig. 11. Geodesic on a cone.
Fig. 12. Four geodesic curves are evaluated on a circular cone in (1). The Coons patch interpolating these curves as boundary geodesics is shown superposed on the cone in (2). Finally, the Coons patch is shown on its own — with and without the departure/arrival vectors — in (3) and (4).
6. Examples The following examples illustrate the application of the existence conditions for surfaces interpolating geodesic boundary curves. A companion paper (Farouki et al., 2008) gives a detailed analysis in the case of Bézier boundary curves, in terms of constraints on the control points that ensure satisfaction of these conditions. 6.1. Analytical geodesics from elementary surfaces We consider here elementary surfaces, which admit analytical expressions for their geodesic curves. These geodesic curves naturally satisfy the constraints (C1), (C2), (C3). The interpolation of the four boundary curves as geodesics is accomplished using the procedure described in Section 5, and we compare the constructed surface patch with the original analytic surface. For simplicity, the parameters αi1 , αi2 and di1 , di2 are all set to zero in these examples.
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
613
Fig. 13. Left and center: two views of the Coons patch interpolant to four great-circle boundary segments on the sphere. Right: an approximation of the whole sphere using 6 patches. Each patch interpolates its boundaries in such a way that they are geodesics of that patch. Note that the construction is not symmetric, since each of the six patches is different from the others.
Fig. 14. Six different surfaces interpolating given geodesic boundary curves, corresponding to different values of the parameters di1 , di2 for i = 1, 2, 3, 4. Surface (n) is associated with the parameters di1 = −8(n − 1) and di2 = 8(n − 1) — these are related to αi1 , αi2 through Eqs. (43).
• Circular cone. Consider a circular cone generated by rotating a straight line about the z-axis. If this line has inclination a with the z-axis, we can parameterize the cone by R(u , v ) = (u cos v , u sin v , Au ) for u 0 and 0 v < 2π , where A = cot a. Geodesic curves on this cone may be parameterized in terms of the angle v as r( v ) =
u0
cos[( v − v 0 ) sin a]
cos v sin v A
where (u 0 , v 0 ) are parameter values such that the geodesic is tangent at r( v 0 ) to the circle of latitude of radius u 0 on the cone (see Fig. 11). Fig. 12 illustrates the construction of a surface patch bounded by four geodesic curves on a circular cone. Although these boundary curves are geodesics of the constructed surface patch, the patch is not (in general) an exact subset of the cone. In Fig. 12, portions of the patch beneath the cone are shown hatched, while those above the cone are shown smoothly shaded. • Sphere. Geodesics on the sphere are simply great-circle arcs. Fig. 13 shows a four-sided patch constructed so as to incorporate four great-circle segments as geodesic boundary curves — again, this patch is not (in general) an exact subset of the sphere. Fig. 13 also illustrates a covering of the entire sphere by six patches. 6.2. Influence of parameters The examples in Fig. 14 illustrate the influence on the surface shape of varying the parameters di1 , di2 for i = 1, 2, 3, 4. These examples show the importance of the smoothing procedures in step 6 of the algorithm, since the shape quality is evidently quite sensitive to these parameters, and ad hoc choices for them can yield surfaces of poor quality. Finally, Fig. 15 shows the outcome of the three different smoothing procedures on a surface constructed using the geodesic boundary data in Fig. 4.
614
R.T. Farouki et al. / Computer Aided Geometric Design 26 (2009) 599–614
Fig. 15. Coons patches interpolating the geodesic boundary curves shown in Fig. 4, as smoothed according to the criteria 1, 2, and 3 in step 6 of the algorithm (two different views of each smoothed surface are presented).
7. Conclusion The problem of constructing a rectangular Coons surface patch from given boundary curves was investigated, under the stipulation that these boundary curves correspond to geodesic curves of the constructed surface. The existence of solutions to this problem was shown to be contingent on the satisfaction of two types of constraints by the given boundary curves. A global (topological) constraint concerns the compatibility of the variation of the boundary curve principal normals with the normal to an oriented surface. Local constraints, applicable to each of the four patch corners, relate the angles at which these curves meet to their curvatures and torsions. For boundary curve data satisfying the above constraints, the construction of geodesic-bounded tensor-product surface patches can be accomplished using the bicubically-blended Coons interpolation scheme. By way of illustration, examples of the construction of such patches on elementary analytic surfaces were presented. When the boundaries are specified as compatible free-form polynomial or rational Bézier curves, the resulting surface is not (in general) a polynomial or rational patch, because of the need to unitize certain vectors in its construction. In a companion paper (Farouki et al., 2008) we relax the requirement of unit arrival/departure vectors, and thereby identify additional constraints on the control points and/or weights of given Bézier boundary curves, so as to obtain polynomial or rational surface patches for which these boundaries are geodesic curves. Acknowledgements This work has been accomplished during the visit of the third author to the Department of Mechanical and Aeronautical Engineering, University of California, Davis. References Bennis, C., Vézien, J.-M., Iglésias, G., 1991. Piecewise surface flattening for non-distorted texture mapping. Computer Graphics 25 (4), 237–246. Coons, S., 1964. Surfaces for Computer Aided Design. Technical Report, MIT (available as AD 663 504 from National Technical Information Service, Springfield, VA 22161). Coons, S., 1974. Surface patches and B-spline curves. In: Barnhill, R., Riesenfeld, R. (Eds.), Computer Aided Geometric Design. Academic Press. Do Carmo, M.P., 1976. Differential Geometry of Curves and Surfaces. Prentice-Hall. Farin, G., 2002. Curves and Surfaces for CAGD, 5th edition. Academic Press. Farouki, R.T., Szafran, N., Biard, L., 2008. Construction of Bézier surface patches with Bézier curves as geodesic boundaries. Computer-Aided Design, submitted for publication. Gordon, W., 1983. An operator calculus for surface and volume modeling. IEEE Computer Graphics and Applications 3, 18–22. Paluszny, M., 2008. Cubic polynomial patches through geodesics. Computer-Aided Design 40 (1), 56–61. Sánchez-Reyes, J., Dorado, R., 2008. Constrained design of polynomial surfaces from geodesic curves. Computer-Aided Design 40 (1), 49–55. Sprynski, N., Szafran, N., Lacolle, B., Biard, L., 2008. Surface reconstruction via geodesic interpolation. Computer-Aided Design 40 (4), 480–492. Struik, D.J., 1976. Lectures on Classical Differential Geometry. Dover (reprint). Tucker III, C.L., 1997. Forming of advanced composites. In: Gutowski, T.G. (Ed.), Advanced Composites Manufacturing. Wiley, New York, pp. 297–372. Wang, G., Tang, K., Tai, C.H., 2004. Parametric representation of a surface pencil with a common spatial geodesic. Computer-Aided Design 36 (5), 447–459.
