- Email: [email protected]
A variational approach to spline curves on surfaces
Computer Aided Geometric Design 22 (2005) 693–709 www.elsevier.com/locate/cagd
A variational approach to spline curves on surfaces Helmut Pottmann ∗ , Michael Hofer Research Group ‘Geometric Modeling and Industrial Geometry’, Vienna University of Technology, Wiedner Hauptstr. 8-10, A-1040 Wien, Austria Received 4 March 2004; accepted 13 January 2005 Available online 15 July 2005
Abstract Given an m-dimensional surface Φ in Rn , we characterize parametric curves in Φ, which interpolate or approximate a sequence of given points pi ∈ Φ and minimize a given energy functional. As energy functionals we study familiar functionals from spline theory, which are linear combinations of L2 norms of certain derivatives. The characterization of the solution curves is similar to the well-known unrestricted case. The counterparts to cubic splines on a given surface, defined as interpolating curves minimizing the L2 norm of the second derivative, are C 2 ; their segments possess fourth derivative vectors, which are orthogonal to Φ; at an end point, the second derivative is orthogonal to Φ. Analogously, we characterize counterparts to splines in tension, quintic C 4 splines and smoothing splines. On very special surfaces, some spline segments can be determined explicitly. In general, the computation has to be based on numerical optimization. Applications of splines on surfaces go beyond geometric modeling, and arise for example also in motion planning and image processing. 2005 Elsevier B.V. All rights reserved. Keywords: Spline interpolation; Spline approximation; Variational curve design; Curves in manifolds; Shape optimization
1. Introduction Curve design using splines is one of the most fundamental topics in CAGD. B-spline curves possess a beautiful shape preserving connection to their control polygon. They allow us the formulation of efficient algorithms for processing, especially subdivision algorithms. Moreover, at least the curves of odd degree * Corresponding author.
E-mail addresses: [email protected] (H. Pottmann), [email protected] (M. Hofer). 0167-8396/$ – see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cagd.2005.06.006
694
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
and maximal smoothness also arise as solutions of variational problems. There is a huge body of literature on these curves, and they received many generalizations (Hoschek and Lasser, 1993). In extension of the standard spline methods, variational curve design has been investigated in a large number of contributions (see (Brunnett et al., 1993; Moreton and Sequin, 1993) and the references therein). It is quite surprising that there seem to be relatively few contributions on the design of spline curves which are restricted to surfaces or more general surfaces (manifolds) in arbitrary dimensions. Of course, these curves can no longer be expressed as linear combinations of control points with suitable basis functions. Moreover, shape preserving schemes on manifolds, e.g., based on subdivision, are different from the solutions of variational problems. In the present paper, we will investigate the latter topic. We will study the functionals which are minimized by splines in affine spaces, but restrict the candidate curves to a given surface. The solution is not as simple and explicit as in the unrestricted case. However, we will find remarkable counterparts to known results. 1.1. Previous work Let us briefly sketch the literature dealing with splines on surfaces. We are focussing here only on curves which have some shape design handles or are generated as solutions of variational problems. Shoemake (1985) introduced spherical counterparts of Bézier curves. He extended de Casteljau’s algorithm to the sphere, replacing straight line segments by geodesic arcs and replacing ratios of Euclidean distances by ratios of geodesic distances. Using this on the 3-sphere in R4 , he generated the spherical component of rigid body motions and applied it to Computer Animation. Soon, it became apparent that the transition to the sphere has the following problem: One can use the spherical de Casteljau algorithm for evaluation of the curve, but one looses the powerful and important subdivision property. Thus, alternative constructions on the sphere have been proposed, see, e.g., (Kim et al., 1995; Jüttler and Wagner, 2002; Ramamoorthi and Barr, 1997). Again motivated by the application of motion design, Park and Ravani (1995) defined Bézier curves on Riemannian manifolds. Moreover, Sprott and Ravani (1997) extended B-spline algorithms to Lie groups and applied them to motion design. However, smoothness properties have not been proved. An algebraic approach to NURBS curves on quadrics, which is also rooted in line geometry and kinematics, has been proposed by Dietz et al. (1993). It has been shown that this method delivers a useful control structure and is suitable for interactive design (Pottmann and Wallner, 2001). The algebraic approach is extendable to low degree algebraic surfaces. In particular, it can be used for motion design (see e.g. (Jüttler and Wagner, 2002; Röschel, 1998)), but it is not suitable for curve design on general surfaces and manifolds. Among the few contributions dealing with subdivision schemes on manifolds we point to papers by L. Noakes (1998, 1999). His results concern C 1 smoothness of extensions of familiar subdivision schemes to manifolds. For example, he proved that the counterpart to the Lane–Riesenfeld algorithm for cubic B-spline curves on Riemannian manifolds yields a limit curve, which is differentiable and its derivative Lipschitz. Very recently, Wallner and Dyn (2005) found a technique to analyze nonlinear schemes based on certain ‘proximity’ conditions to known linear schemes. Their universally applicable results give remarkable extensions of known work; for example, it turns out that the Riemannian extension of cubic subdivision gives actually C 2 curves. We will not pursue subdivision in this paper, but concentrate on variational design. The ‘intersecting topic’, variational subdivision for curves on surfaces, has been computationally addressed, but not analyzed in (Hofer et al., 2003).
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
695
A number of papers dealing with variational design on surfaces is considering the sphere (Brunnett and Crouch, 1994; Brunnett et al., 1995), partially in view of the application to motion design (Park and Ravani, 1997; Ramamoorthi and Barr, 1997). An early contribution to variational curve design on more general surfaces is due to Noakes et al. (1989). The authors characterize the minimizers of an intrinsic geometric counterpart to the L2 norm of the second derivative. This is the integral of the squared covariant derivative of the first derivative with respect to arc length. For recent results on these intrinsic splines, see (Noakes, 2003; Camarinha et al., 2001). We will not pursue the intrinsic formulation in the present paper. The most closely related work to our contribution is the PhD thesis of H. Bohl (1999). It deals with the usual energy functionals used in CAGD, but restricts the curves to a given surface in R3 . Bohl proves the existence of a solution and gives computed examples on parametric surfaces, based on a quasi Newton optimization algorithm. A large part of his work deals with patch boundaries. If we consider on a surface the minimizers of the curve length, we obtain the well-studied geodesics. Their geometric properties have been investigated in classical differential geometry. Geodesics in a scaled arc length parameterization also arise as minimizers of the L2 norm of the first derivative. A variety of applications of geodesics has been described within Computer Vision and Image Processing (see e.g. (Osher and Paragios, 2003; Sapiro, 2001)). The present work has also been inspired by research on active contours (Blake and Isard, 1998), especially by work on active curves and geometric flows of curves on surfaces (see, e.g. (Cheng et al., 2002; Memoli and Sapiro, 2001; Osher and Paragios, 2003)). 1.2. Contributions and outline of the present paper We will study the minimizers of standard energy functionals of spline theory, which interpolate or approximate given data points, and are restricted to a given m-dimensional surface Φ ⊂ Rn . The restriction to Φ is the new aspect, and it is also the source of arising complications. However, we will find nice characterizations of the minimizers. In Section 2, we characterize the counterparts on surfaces to C 2 cubic splines, to splines in tension and to C 4 quintic splines. It will turn out that differentiability at the interpolation points is the familiar one, but the characterization of spline segments is slightly different. We give an example: Whereas the minimizers of the L2 norm of the second derivative have cubic segments (vanishing fourth derivative) in the unrestricted case, the corresponding splines on surfaces have segments with vanishing tangential component of the fourth derivative; we call such segments ‘tangentially cubic’. Hence, the differential equation c(4) = 0 changes to tpr c(4) = 0, with tpr denoting the tangential projection (orthogonal projection into the corresponding tangent space of Φ). Analogous changes are found in the characterizing differential equations for the other spline schemes. Section 3 is devoted to approximating curves, in particular to the counterparts of smoothing splines. In view of the importance of tangentially cubic curves, we give in Section 4 some explicit representations of such curves on special surfaces, namely certain cylinder surfaces. In particular, we address special tangentially cubic parametrizations of a curve; these are parametrizations whose fourth derivative is orthogonal to the curve. Computational algorithms are not the focus of the present paper. We address the computation just briefly in Section 5, where we also conclude the paper with examples and pointers to a number of interesting applications in Geometric Modeling, Computer Graphics, Robotics and Image Processing.
696
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
2. Interpolating spline curves on surfaces In view of the applications we have in mind, it is necessary to work on surfaces of arbitrary dimension and codimension. Thus, we consider an m-dimensional surface Φ in Euclidean Rn , m < n. Moreover, a sequence of points pi ∈ Φ, i = 1, . . . , N , and real numbers u1 < · · · < uN shall be given. We are seeking interpolating splines on the surface Φ. Sometimes we will also call Φ a manifold; if not stated otherwise, we work with a manifold without boundary. Thus, we have closed or unbounded surfaces. As we will point out in Section 5, the computation of energy-minimizing splines which respect patch boundaries or holes in a surface, can be reduced to the computation of a spline on (another) surface without boundary. 2.1. The counterparts to cubic splines Let us recall the situation, where we are not confined to a surface: Among all curves x(u) ⊂ Rn , whose first derivative is absolutely continuous, x˙ ∈ AC(I ), and whose second derivative satisfies x¨ ∈ L2 (I ) on I = [u1 , uN ], and which interpolate the given data, x(ui ) = pi , the unique minimizer of uN 2 (1) E2 (x) = x¨ (u) du, u1
is the interpolating C 2 cubic spline c(u). In the following we would like to extend this well-known result to the case where the admissible curves x(u) are restricted to the given surface Φ. We are not changing the functional E2 , whose interpretation requires an embedding space. We are considering the restriction to the surface as a constraint, rather than formulating the problem in terms of the intrinsic geometry of the manifold. As we will see later, this is a suitable formulation for the problems we would like to solve. Theorem 1. Consider real numbers u1 < · · · < uN and points p1 , . . . , pN on an m-dimensional C 4 surface Φ in Euclidean Rn . We let I = [u1 , uN ]. Then among all C 1 curves x : I → Φ ⊂ Rn , which interpolate the given data, i.e., x(ui ) = pi , i = 1, . . . , N , and whose restrictions to the intervals [ui , ui+1 ], i = 1, . . . , N − 1, are C 4 , a curve c which minimizes the functional E2 of Eq. (1) is C 2 and possesses segments c|[ui , ui+1 ], whose fourth derivative vectors are orthogonal to Φ. Moreover, at the end points p1 = c(u1 ) and pN = c(uN ) of the solution curve, the second derivative vector is orthogonal to Φ. Proof. If a solution curve c exists, the first variation of the energy functional must vanish there (Gelfand and Fomin, 1963). To express this condition, we consider neighboring curves x(u) ⊂ Φ, written as x(u, ε) = c(u) + h(u, ε).
(2)
For any fixed u˜ ∈ I , the curve h(u, ˜ ε) is a curve in Φ with h(u, ˜ 0) = 0. Its Taylor expansion at ε = 0 reads ε2 ˜ 0) + hεε (u, ˜ 0) + · · · , h(u, ˜ ε) = εhε (u, 2 ˜ 0) =: t(u) ˜ is a tangent vector of Φ at c(u). ˜ where the subscripts indicate differentiation. Note that hε (u, The displacement curves h to the given interpolation points vanish, h(ui , ε) = 0, for all ε. In particular,
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
697
Fig. 1. Displacement of the solution curve for deriving the first variation of the energy.
we have hε (ui , 0) = t(ui ) = 0. For the following, it is important to see that the mixed partial derivative vector hεu (ui , 0) = tu (ui ) is a tangent vector of Φ at c(ui ) = pi . Geometrically, this follows from the fact that the curve c(u) + t(u) is a curve on the ruled surface r(u, v) = c(u) + vt(u), which touches Φ along c. This curve passes through each interpolation point pi = c(ui ) + t(ui ), and thus its derivative vector cu (ui ) + tu (ui ) is a tangent vector of Φ. Together with the tangency of cu (ui ) this implies tangency of tu (ui ). Whatever field of displacement curves h(u, ε) we have chosen, the energy functional must assume a stationary value at ε = 0, d E2 x(u, ε) ε=0 = 0. (3) dε In view of 2 2 c¨ (u) + huu (u, ε) du = c¨ + εtuu + ε2 (· · ·) du, E2 x(u, ε) = I
I
this is equivalent to c¨ · tuu du = 0.
(4)
I
Integration by parts on each segment yields ui+1 ui+1 ui+1 (3) ui+1 c¨ · tuu du = c¨ · tu |ui − c · t|ui + c(4) · t du. ui
ui
We first note that the middle term vanishes, since t(ui ) = 0. Summing up over all intervals and denoting left and right derivatives at the knots ui by a subscript − and +, respectively, the condition for a solution c reads, N −1 (5) c¨ − (ui ) − c¨ + (ui ) · tu (ui ) + c¨ − (uN ) · tu (uN ) + c(4) · t du. 0 = −¨c+ (u1 ) · tu (u1 ) + i=2
I
698
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
Since t(u) is an arbitrary tangent vector field along c, and tu (ui ) are arbitrary tangent vectors at the given interpolation points, a solution curve c(u) must satisfy the following orthogonality conditions to the surface Φ. To formulate them, we denote by tpr the orthogonal projection of a vector at a point p ∈ Φ onto the tangent space τ of Φ at p, tpr c¨ + (u1 ) = tpr c¨ − (uN ) = 0, tpr c¨ − (ui ) − c¨ + (ui ) = 0, i = 2, . . . , N − 1,
(6)
tpr c(4) (u) = 0 on each interval [ui , ui+1 ].
(8)
(7)
Exactly these conditions are stated in Theorem 1. Eq. (6) means that the second derivative vector at the end points is orthogonal to the surface; this implies (but is not equivalent to) vanishing geodesic curvature at the end points. Eq. (7) shows that the tangential component of the second derivative is continuous along c; we call this behaviour tangentially C 2 or briefly T C 2 . It will be shown later that in view of the C 4 manifold Φ this implies C 2 . Finally, by Eq. (8), the fourth derivative vectors of the curve (which may be discontinuous at the interpolation points) are orthogonal to the manifold Φ. Vanishing fourth derivative characterizes a cubic curve; since here only the tangential component vanishes, we speak of a tangentially cubic curve in the following. To complete the proof, we show the following result: Lemma 2. On an m-dimensional C k manifold Φ ⊂ Rn , m < n, we consider a curve c(u), which is C k (k 2) everywhere except at a point c(u0 ). At c(u0 ), all derivatives up to order k have a continuous tangential component, i.e., (i) i = 1, . . . , k. (9) tpr c(i) − (u0 ) − c+ (u0 ) = 0, Then, the curve is also C k at c(u0 ). Proof. We may represent the manifold at least in a neighborhood of c(u0 ) as an intersection of n − m hypersurfaces, whose implicit representations shall be Fj (x) = 0,
j = 1, . . . , n − m.
Since the curve lies in all these hypersurfaces, we have the identity Fj (c(u)) = 0 in u. By differentiation, we obtain ∇Fj (c) · c˙ = 0, and for the ith derivative we find an expression of the form Gi−1 (c, . . . , c(i−1) ) + ∇Fj (c) · c(i) = 0.
(10)
Here, Gi−1 (c, . . . , c(i−1) ) abbreviates a function, which depends on Fj , up to differentiation order i, and on c, but only up to differentiation order i − 1. At c(u0 ) this holds for the left and right derivatives. The first derivative is tangential, and thus c is C 1 at u0 . This is the basis of an induction on the differentiation order i. Assuming continuity of the derivative c(i−1) at u0 , Eq. (10) expresses that the components ∇Fj (c) · c(i) of the ith derivative are continuous at u0 , since they equal Gi−1 (c, . . . , c(i−1) ), which is continuous by the induction hypothesis. Since the vectors ∇Fj (c(u0 )), j = 1, . . . , n − m, span the normal space of Φ at c(u0 ), we have shown that the normal component of the ith derivative is continuous. The tangential component is continuous by (9) and thus we have shown continuity of c(i) . This completes the proof of Lemma 2 and thus also the proof of Theorem 1. Note that the case k = 2 follows also immediately by Meusnier’s theorem. 2
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
699
Fig. 2. For a minimizer of (left) E1 and (right) E2 , the second respectively fourth derivative vectors (yellow) are orthogonal to the surface. (For interpretation of the references in color in this figure legend, the reader is referred to the web version of this article.)
The minimizers possess a characterization which is very similar to the familiar cubic C 2 splines. They are C 2 and tangentially cubic. However, the problem is nonlinear. Tangentially cubic curves can in general not be computed explicitly, but only by a numerical algorithm. Existence of the solution follows from the work of Bohl (1999) and Wallner (2004). Uniqueness cannot be expected, as will be clear from the following considerations. 2.2. Geodesics Let us replace the energy in (1) by the L2 norm of the first derivative, uN 2 E1 (x) = x˙ (u) du.
(11)
u1
Moreover, we just prescribe the two end points. Then we find with the same approach as in the proof of Theorem 1 that the curve’s second derivative vectors c¨ are orthogonal to Φ; we may say the curve is tangentially linear. In particular, c˙ and c¨ are orthogonal, and hence ˙c2 = const. This proves, that the curve is parameterized by a constant multiple of its arc length. Moreover, it shows that c¨ represents the principal normal, and orthogonality of the principal normal to Φ characterizes a geodesic. Thus we find the known result that the minimizers of E1 are geodesics in a scaled arc length parameterization. Note that the functionals we are considering are also optimizing the parameterization. Since even for geodesics we do not have simple results on uniqueness, the more involved case of splines will hardly allow us a characterization of situations with a unique solution.
700
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
Fig. 3. Spline curves c1 , ct , c2 on a surface Φ, interpolating the same points p1 , . . . , p4 and minimizing energies E1 , Et , and E2 , respectively.
2.3. Splines in tension on surfaces A linear combination of energies (11) and (1) leads in the unrestricted case to the well known splines in tension as minimizers (Schweikert, 1966). The counterpart on manifolds is characterized in the following theorem, whose proof can be omitted since it is completely analogous to that of Theorem 1. Theorem 3. Consider a sequence of data points, parameter values and admissible curves on a surface Φ as in Theorem 1. Then, a minimizer of the functional uN (12) Et (x) = (¨x2 + wx˙ 2 ) du, w = const > 0, u1
is a C 2 curve which satisfies tpr c(4) (u) − w¨c(u) = 0
(13)
on all segments. The end conditions are tpr c¨ + (u1 ) = tpr c¨ − (uN ) = 0. Increasing w brings the solution closer to the curve which minimizes E1 , i.e., the curve composed by geodesic segments between the interpolation points. Hence, the parameter w controls the tension of the curve, see Fig. 3. 2.4. Counterparts to higher order splines on surfaces It is known from spline theory, that the minimizers of the L2 norm of the third derivative, uN 2 E3 (x) = x(3) (u) du, u1
(14)
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
701
under interpolation conditions are quintic C 4 splines. It is not difficult to extend this result to manifolds. Theorem 4. Consider real numbers u1 < · · · < uN and points p1 , . . . , pN on an m-dimensional C 6 manifold Φ in Euclidean Rn . Then among all C 2 curves x : [u1 , uN ] → Φ ⊂ Rn , which interpolate the given data and whose restrictions to the intervals [ui , ui+1 ], i = 1, . . . , N − 1, are C 6 , a curve c which minimizes the L2 norm E3 of the third derivative is C 4 and possesses segments c|[ui , ui+1 ], whose sixth derivative vectors are orthogonal to Φ. Moreover, at the end points p1 = c(u1 ) and pN = c(uN ) of the solution curve, the third derivative vanishes and the fourth derivative vector is orthogonal to Φ. Proof. The proof resembles the one of Theorem 1. We use the same notation and just mention the few essential differences. A solution c(u) must satisfy c(3) · tuuu du = 0. (15) I
Integration by parts on each segment yields ui+1 ui+1 (3) (3) ui+1 (4) ui+1 (5) ui+1 c · tuuu du = c · tuu |ui − c · tu |ui + c · t|ui − c(6) · t du. ui
ui
From the four terms on the right hand side, the third one vanishes because of the interpolation conditions, expressed by t(ui ) = 0. Since t is an arbitrary tangent vector field along c, the last term vanishes iff tpr c(6) = 0. Now we can sum up over all intervals and collect the contributions from the first two terms c(3) · tuu and c(4) · tu . Since tu (ui ) is tangential, we find a T C 4 condition at the inner knots and orthogonality of the fourth derivative to the manifold at the end points. It is not hard to see that tuu (ui ) can be an arbitrary vector and this shows that a solution must have vanishing third derivative at the ends and it must be C 3 at the inner knots. Application of Lemma 2 shows that the curve is even C 4 at the inner knots, which completes the proof. 2 Remark 5. So far we always considered an open curve and natural end conditions. It is clear from the derivation of our results how to handle other cases. Let us discuss this at hand of the minimization of E3 . The segments are of course always tangentially quintic, i.e., tpr c(6) = 0. A closed spline is C 4 at all knots. Moreover, instead of working for an open curve with natural end conditions, we could prescribe first and second derivative vectors at the end points. The extension of the results to the minimization of the L2 norm of even higher derivatives, or combinations of L2 norms of different derivatives, is straightforward and thus it can be omitted.
3. Approximating spline curves on surfaces Reinsch (1967) relaxed the interpolation conditions x(ui ) − pi = 0 by adding the sum of squared interpolation errors as a penalty term to E2 , uN N 2 s x(ui ) − pi 2 , λ, µ > 0. ¨ x(u) du + µ (16) E2 (x) = λ u1
i=1
702
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
The minimizers, called smoothing splines, found a variety of applications in the analysis of observational data (Wahba, 1990). The choice of the smoothing parameter λ : µ in this method for data approximation is not a simple problem. Often one uses statistical methods for this critical task (Wahba, 1990). We should also mention that an analogous functional is used for freeform surface fitting in reverse engineering applications (Várady and Martin, 2002). We can use the same method as in the proof of Theorem 1 to derive the following result on the counterparts to cubic smoothing splines on surfaces. Theorem 6. Consider a sequence of data points, parameter values and differentiability conditions on an admissible curve on a surface Φ as in Theorem 1. Then, a minimizer of the smoothing spline functional E2s has all the properties of the minimizer of E2 in Theorem 1. However, the interpolation conditions c(ui ) = pi are replaced by the following transition conditions between segments, (3) = 0. (17) tpr λ c(3) + (ui ) − c− (ui ) + µ c(ui ) − pi In an analogous way, one can investigate the minimizers of smoothing spline counterparts Ets , E3s , . . . , to other energy functionals. For example, the case of E3s yields in modification of Theorem 4 curves, where the interpolation conditions are replaced by the smoothness conditions (5) = 0. (18) tpr λ c(5) − (ui ) − c+ (ui ) + µ c(ui ) − pi 4. Examples of tangentially cubic curves 4.1. Tangentially cubic curves on cylinder surfaces in R3 Let us consider a general cylinder surface Φ ⊂ R3 . We may choose a directrix curve l(u) = (x1 (u), x2 (u), 0) and rulings parallel to (0, 0, 1), i.e., a parameterization of Φ in the form Φ : x(u, v) = (x1 (u), x2 (u), v). For a tangentially cubic spline curve c(t) = (c1 (t), c2 (t), c3 (t)) ⊂ Φ, the fourth derivative c(4) must be orthogonal to Φ, and thus also orthogonal to the corresponding ruling of Φ. This requires a vanishing third component c3(4) (t) = 0 and shows that the function c3 is cubic, c3 (t) = a0 + a1 t + · · · + a3 t 3 . The other two coordinate functions (c1 (t), c2 (t)) have its fourth derivative orthogonal to a parameterization (c1 (t), c2 (t), 0) of the (planar) orthogonal cross section l of Φ. This is therefore a tangentially cubic parameterization of the cross section curve. Choosing the data points for a spline on a cross section, we find a result on energy minimizing parameterizations of a given curve. Corollary 7. Given a C 4 curve l in Rn , and a sequence of points pi with parameters ui , i = 1, . . . , N, on it. Then, among all C 1 , piecewise C 4 parameterizations of the curve which satisfy the interpolation conditions c(ui ) = pi , a minimizer of the L2 norm of the second derivative is tangentially cubic, i.e., a parameterization whose fourth derivative is orthogonal to l. Moreover, c is C 2 and the second derivative at the end points is orthogonal to l.
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
703
The derivation in Rn can be done with a cylinder surface in Rn+1 . We may also note that the case of a manifold of dimension 1 is actually included in the formulation and in the proof of Theorem 1. Since the study of splines on cylinders is reduced to tangentially cubic parameterizations of curves, we will study those in the following in more detail. 4.2. Tangentially cubic arc length parameterizations In view of the important role of arc length parameterizations for curves, we answer the question on which curves in Rn the arc length parameterization c(s) (or a scaled version of it), is tangentially cubic. Denoting the Frenet frame of c by e1 (s), . . . , en (s), and the curvatures by κ1 , . . . , κn−1 , we have c = e1 ,
e1 = κ1 e2 ,
e1 = κ1 e2 − κ 2 e1 + κ1 κ2 e3 .
Hence, the condition for a tangentially cubic parameterization becomes 0 = c(4) · c = e(3) 1 · e1 = −3κ1 κ1 .
Proposition 8. The arc length parameterization c(s) of a curve in Rn is tangentially cubic iff the curve possesses constant curvature κ1 . In the plane, we get only circles and straight lines. However, already in R3 this family is relatively rich, since the torsion τ = κ2 can be arbitrary. The most prominent examples are helices which even have constant torsion. 4.3. Tangentially cubic parameterizations of the circle; spline curves on a right circular cylinder Consider the unit circle x2 = 1 and its parameterization c(t) = cos φ(t), sin φ(t) . We want to determine φ(t) such that c(t) is a tangentially cubic parameterization of the circle. A simple computation shows that orthogonality of first and fourth derivative is equivalent to (19) φ (4) − 6φ˙ 2 φ¨ = 0. ˙ and thus have to solve the differential equation We set ψ(t) := φ(t) ˙ ψ (3) = 6ψ 2 ψ. We can immediately perform a first integration, ψ¨ = 2ψ 3 + A0 ,
(20)
with a constant A0 ∈ R. Since the particular solution ψ˙ = 0, i.e., the scaled arc length parameterization φ = at + b is already known from the previous considerations, we may assume ψ˙ = 0 from now on. Multiplying Eq. (20) by 2ψ˙ and integrating again, we arrive with A := 2A0 at ψ˙ 2 = ψ 4 + Aψ + B. Separation of the variables and integration leads to dψ . t +C =± ψ 4 + Aψ + B
(21)
704
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
Fig. 4. Examples of tangentially cubic curves with parameterization (23), (a) and (b), and parameterization (24), (c) and (d).
The integral on the right hand side requires Legendre normal integrals. We omit its discussion in the general case and just look at the special case A = B = 0. There, we obtain ψ = ±1/(t + C) and φ(t) = ln |C ± t| + D.
(22)
Thus, even in the very simple case of a right circular cylinder, just two very special families of tangentially cubic curves can be written down explicitly. On the unit cylinder x12 + x22 = 1, these are the curves derived from a scaled arc length parameterization c(t) = cos(at + b), sin(at + b), a0 + a1 t + · · · + a3 t 3 , (23) and the curves to the circle parameterization (22), c(t) = cos ln |C ± t| + D , sin ln |C ± t| + D , a0 + a1 t + · · · + a3 t 3 .
(24)
A few examples of curves (23) and (24) are depicted in Fig. 4. 4.4. Tangentially cubic parameterization by the direction angle of the curve tangent It turns out that nice explicit representations of tangentially cubic curves on cylinder surfaces arise as follows. We would like to determine those planar curves, whose parameterization c(φ) with respect to the direction angle φ of the curve tangent is tangentially cubic. Therefore, we write the family of tangent lines of such a curve in the form −x1 sin φ + x2 cos φ + h(φ) = 0.
(25)
Here, φ is the angle between the tangent, oriented by the vector (cos φ, sin φ), and the x1 -axis; h denotes the signed distance of the origin to the oriented tangent. The function h(φ) is known as support function. The curve parameterization c(φ) is computed as envelope of the family of tangents, i.e., from Eq. (25) and its derivative with respect to φ. We find c(φ) = (h sin φ + h˙ cos φ, −h cos φ + h˙ sin φ). (26) A tangentially cubic parameterization requires c˙ · c(4) = 0, which is equivalent to h(5) − 2h(3) − 3h˙ = 0.
(27)
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
The solution of this differential equation is, with real constants bi , √ h(φ) = b0 + b1 cos φ + b2 sin φ + b3 cosh( 3φ + b4 ).
705
(28)
To understand these curves, we note that an appropriate translation of the origin can eliminate the term b1 cos φ + b2 sin φ, and b4 just accounts for a rotation. Thus it is sufficient to study the case b1 = b2 = b4 = 0. Addition of the constant b0 to the support function means construction of the offset at signed distance b0 . The curve with support function √ (29) h(φ) = cosh 3φ, is a special instance of a so-called hypercycloid. Hypercycloids are a certain counterpart to the familiar epicycloids or hypocycloids, which are obtained as locus of a point on a circle, which rolls on another circle. A hypercycloid possesses such a kinematic generation, but the generating circles are not real: the moving circle has a complex radius and the fixed circle a purely imaginary radius. Hypercycloids arise for example as orthogonal projections of geodesics on a paraboloid of revolution, when projecting parallel to the rotational axis (see (Strubecker, 1969, pp. 227–232)). All other nontrivial solutions of our problem arise from the special hypercycloid with support function (29) by offsetting and a similarity transformation. The trivial solutions belong to b3 = 0 and are circles. On a cylinder surface with the hypercycloid to (29) as cross section, the curves √ √ √ cosh 3φ sin φ + 3 sinh 3φ cos φ √ √ √ c(φ) = − cosh 3φ cos φ + 3 sinh 3φ sin φ (30) a0 + a1 φ + a2 φ 2 + a3 φ 3 are tangentially cubic. Fig. 5 shows for φ ∈ [−2.5, 2.5] tangentially cubic curves (30) on a cylinder surface with the hypercycloid (29) as cross section, with (a) a0 = a1 = a3 = 0, a2 = 5 and (b) a0 = 0, a1 = 1, a2 = 2, a3 = 1.
Fig. 5. Tangentially cubic curves (30) on a cylinder surface with the hypercycloid (29) as cross section.
706
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
5. Computation based on an optimization algorithm The presented spline curves possess nice geometric characterizations. However, we found out in Section 4 that the problem is not even explicitly solvable on very simple surfaces. Hence, we have to use numerical algorithms based on a discretization. Geometrically this means that we replace the curve c(t) ⊂ Φ by a sufficiently dense polygon Pc with vertices ci ∈ Φ. Discretizing a quadratic functional, such as E1 , E2 , Et , results in a quadratic function in the coordinates of the vertices ci of the approximating polygon Pc . Since the polygon vertices are restricted to Φ, we end up with the constrained minimization of a quadratic function in a rather high dimensional space. From a geometric viewpoint, we have to construct the closest point p∗ on some manifold S to a given point p. Here, p represents the unrestricted (discretized) spline curve in Rn . Basically, one can use any algorithm for constrained optimization (Geiger and Kanzow, 2002). If the surface is given in parametric form, we can of course eliminate the constraints by the use of a polygon in the parameter domain. Then, an appropriate algorithm for unconstrained optimization over the parameter domain can be employed (Kelley, 1999). For example, H. Bohl (1999) uses the standard BFGS algorithm. However, several other efficient algorithms (quasi-Newton methods, Newton method with improvement of global behaviour, etc. (Kelley, 1999)) are available. We have developed a geometrically motivated optimization algorithm for the solution of the present optimization problem (Hofer and Pottmann, 2004). It extends and stabilizes the simple algorithm for variational curve design in (Hofer et al., 2003), and can also be interpreted within the Sobolev gradient method (Neuberber, 1997). This algorithm can be used for the computation of splines on parametric surfaces, level sets, triangle meshes and point set surfaces, see (Hofer and Pottmann, 2004). Fig. 3 depicts energy minimizing curves on a parametric surface. Fig. 6 shows energy minimizing curves on a point set surface and on a level set. All of them have been computed with the optimization algorithm presented in (Hofer and Pottmann, 2004).
Fig. 6. Curve minimizing E2 on a point set surface (left), and level set (right).
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
707
Another computational approach could directly use the characterization of the splines in a discretized way. For example, a spline minimizing E2 has 4th central difference vectors, which are orthogonal to Φ. Formulation of this property and the other characterizations gives a system of nonlinear equations, which can be solved with a Newton iteration. One will not apply the pure Newton algorithm, but use one of the known methods for improving the global convergence behaviour (Kelley, 1999). Remark 9. For applications in CAD, the surface Φ and the preimage of the desired spline curve in the parameter domain, have to be in NURBS form. To achieve a NURBS representation of the energy minimizing curve, we can use numerical optimization of the energy over a finite-dimensional linear space of B-spline curves in the parameter domain. Another approach separates the computation of the shape from the computation of a B-spline representation. One first computes the energy minimizing curve, represented by a polygon, with the numerical methods indicated above; then the resulting polygon vertices with their parameter values are approximated by one of the known schemes from CAGD. For the state of the art on the latter topic, we refer to (Renner and Weiss, 2004). 5.1. Respecting boundaries and holes Throughout the paper we assumed surfaces without boundary. This might appear as a severe limitation, especially in view of the frequent presence of boundaries in CAD related applications. Holes and boundaries can be treated as obstacles in a curve design problem. Fortunately, variational curve design in the presence of obstacles can be reformulated as curve design on another surface, which does not have a boundary. This method, a geometric version of the barrier method from constrained optimization (Geiger and Kanzow, 2002), is described and illustrated in (Hofer and Pottmann, 2004). 5.2. Applications Energy minimizing splines on surfaces possess a variety of applications. Partially, these go beyond standard CAGD applications. Some major applications are the following ones (Hofer and Pottmann, 2004): Curve design on the curved manifold, which represents the Lie group of Euclidean rigid body motions, leads to variational motion design for animation in Computer Graphics and motion planning problems in robotics. Applying the concept to so-called image manifolds, as introduced in (Kimmel et al., 2000), one can use splines for image sensitive drawing and for interactive image segmentation.
Acknowledgement This work was supported by the Austrian Science Fund (FWF) under grant P16002-N05.
References Blake, A., Isard, M., 1998. Active Contours. Springer, Berlin. Bohl, H., 1999. Kurven minimaler Energie auf getrimmten Flächen. PhD thesis, Univ. Stuttgart, Germany. Brunnett, G., Crouch, P.E., 1994. Elastic curves on the sphere. Adv. Comp. Math. 2, 23–40.
708
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
Brunnett, G., Crouch, P.E., Silva Leite, F., 1995. Spline elements on spheres. In: Daehlen, M., Lyche, T., Schumaker, L. (Eds.), Mathematical Methods for Curves and Surfaces. Vanderbilt Univ. Press, pp. 49–54. Brunnett, G., Hagen, H., Santarelli, P., 1993. Variational design of curves and surfaces. Surv. Math. Ind. 3, 1–27. Camarinha, M., Silva Leite, F., Crouch, P.S., 2001. On the geometry of Riemannian cubic polynomials. Differential Geom. Appl. 15, 107–135. Cheng, B.-T., Burchard, P., Merriman, B., Osher, S., 2002. Motion of curves constrained on surfaces using a level set approach. J. Comput. Phys. 175 (2), 604–644. Dietz, R., Hoschek, J., Jüttler, B., 1993. An algebraic approach to curves and surfaces on the sphere and on other quadrics. Computer Aided Geometric Design 10, 211–229. Geiger, C., Kanzow, C., 2002. Theorie und Numerik restringierter Optimierungsaufgaben. Springer, Heidelberg. Gelfand, I.M., Fomin, S.V., 1963. Calculus of Variations. Prentice-Hall, Englewood Cliffs, NJ. Hofer, M., Pottmann, H., Ravani, B., 2003. Geometric design of motions constrained by a contacting surface pair. Computer Aided Geometric Design 20, 523–547. Hofer, M., Pottmann, H., 2004. Energy-minimizing splines in manifolds. In: Proceedings of ACM SIGGRAPH 2004, Trans. Graph. 23 (3), 284–293. Hoschek, J., Lasser, D., 1993. Fundamentals of Computer Aided Geometric Design. A.K. Peters, Wellesley, MA. Jüttler, B., Wagner, M., 2002. Kinematics and animation. In: Farin, G., Kim, M.-S., Hoschek, J. (Eds.), Handbook of Computer Aided Geometric Design. North Holland, Amsterdam, pp. 723–748. Kelley, C.T., 1999. Iterative Methods for Optimization. SIAM, Philadelphia. Kim, M.-J., Kim, M.-S., Shin, S., 1995. A general construction scheme for unit quaternion curves with simple high order derivatives. In: SIGGRAPH ’95, Comput. Graph. 29, 369–376. Kimmel, R., Malladi, R., Sochen, N., 2000. Images as embedded maps and minimal surfaces: movies, color, texture and volumetric medical images. Internat. J. Computer Vision 39, 111–129. Memoli, F., Sapiro, G., 2001. Fast computation of weighted distance functions and geodesics on implicit hyper-surfaces. J. Comp. Phys. 173, 730–764. Moreton, H., Sequin, C., 1993. Scale-invariant minimum-cost curves: fair and robust design implements. Comput. Graph. Forum 12, 473–484. Neuberber, J.W., 1997. Sobolev Gradients and Differential Equations. Springer, Heidelberg. Noakes, L., Heinzinger, G., Paden, B., 1989. Cubic splines on curved spaces. IMA J. Math. Control Inform. 6, 465–473. Noakes, L., 1998. Non-linear corner cutting. Adv. Comp. Math. 8, 165–177. Noakes, L., 1999. Accelerations of Riemannian quadratics. Proc. Amer. Math. Soc. 127, 1827–1836. Noakes, L., 2003. Null cubics and Lie quadratics. J. Math. Phys. 44 (3), 1436–1448. Osher, S., Paragios, N. (Eds.), 2003. Geometric Level Set Methods. Springer, New York. Park, F.C., Ravani, B., 1995. Bézier curves on Riemannian manifolds and Lie groups with kinematic applications. ASME J. Mech. Design 115, 36–40. Park, F.C., Ravani, B., 1997. Smooth invariant interpolation of rotations. ACM Trans. Graph. 16, 277–295. Pottmann, H., Wallner, J., 2001. Computational Line Geometry. Springer, Heidelberg. Ramamoorthi, R., Barr, A., 1997. Fast construction of accurate quaternion splines. In: SIGGRAPH ’97, Comput. Graph. 31, 287–292. Reinsch, C.H., 1967. Smoothing by spline functions. Numer. Math. 10, 177–183. Renner, G., Weiss, V., 2004. Exact and approximate computation of B-spline curves on surfaces. Computer-Aided Design 36, 351–362. Röschel, O., 1998. Rational motion design—a survey. Computer-Aided Design 30, 169–178. Sapiro, G., 2001. Geometric Partial Differential Equations and Image Analysis. Cambridge Univ. Press, Cambridge. Schweikert, D., 1966. An interpolation curve using a spline in tension. J. Math. Phys. 45, 312–317. Shoemake, K., 1985. Animating rotations with quaternion curves. In: SIGGRAPH ’85, Comput. Graph. 19, 245–254. Sprott, K., Ravani, B., 1997. Ruled surfaces, Lie groups and mesh generation. In: Proc. ASME Design Eng. Techn. Conf. No DETC97/DAC-3966. Strubecker, K., 1969. Differentialgeometrie II: Theorie der Flächenmetrik. Walter de Gruyter, Berlin. Várady, T., Martin, R., 2002. Reverse Engineering. In: Farin, G., Hoschek, J., Kim, M.-S. (Eds.), Handbook of Computer Aided Geometric Design. Elsevier, Amsterdam, pp. 651–681. Wahba, G., 1990. Spline Models for Observational Data. SIAM, Philadelphia.
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
709
Wallner, J., 2004. Existence of set-interpolating and energy-minimizing curves. Computer Aided Geometric Design 21 (9), 837–857. Wallner, J., Dyn, N., 2005. Convergence and C 1 analysis of subdivision schemes on manifolds by proximity. Computer Aided Geometric Design 22. In press.
A variational approach to spline curves on surfaces Helmut Pottmann ∗ , Michael Hofer Research Group ‘Geometric Modeling and Industrial Geometry’, Vienna University of Technology, Wiedner Hauptstr. 8-10, A-1040 Wien, Austria Received 4 March 2004; accepted 13 January 2005 Available online 15 July 2005
Abstract Given an m-dimensional surface Φ in Rn , we characterize parametric curves in Φ, which interpolate or approximate a sequence of given points pi ∈ Φ and minimize a given energy functional. As energy functionals we study familiar functionals from spline theory, which are linear combinations of L2 norms of certain derivatives. The characterization of the solution curves is similar to the well-known unrestricted case. The counterparts to cubic splines on a given surface, defined as interpolating curves minimizing the L2 norm of the second derivative, are C 2 ; their segments possess fourth derivative vectors, which are orthogonal to Φ; at an end point, the second derivative is orthogonal to Φ. Analogously, we characterize counterparts to splines in tension, quintic C 4 splines and smoothing splines. On very special surfaces, some spline segments can be determined explicitly. In general, the computation has to be based on numerical optimization. Applications of splines on surfaces go beyond geometric modeling, and arise for example also in motion planning and image processing. 2005 Elsevier B.V. All rights reserved. Keywords: Spline interpolation; Spline approximation; Variational curve design; Curves in manifolds; Shape optimization
1. Introduction Curve design using splines is one of the most fundamental topics in CAGD. B-spline curves possess a beautiful shape preserving connection to their control polygon. They allow us the formulation of efficient algorithms for processing, especially subdivision algorithms. Moreover, at least the curves of odd degree * Corresponding author.
E-mail addresses: [email protected] (H. Pottmann), [email protected] (M. Hofer). 0167-8396/$ – see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cagd.2005.06.006
694
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
and maximal smoothness also arise as solutions of variational problems. There is a huge body of literature on these curves, and they received many generalizations (Hoschek and Lasser, 1993). In extension of the standard spline methods, variational curve design has been investigated in a large number of contributions (see (Brunnett et al., 1993; Moreton and Sequin, 1993) and the references therein). It is quite surprising that there seem to be relatively few contributions on the design of spline curves which are restricted to surfaces or more general surfaces (manifolds) in arbitrary dimensions. Of course, these curves can no longer be expressed as linear combinations of control points with suitable basis functions. Moreover, shape preserving schemes on manifolds, e.g., based on subdivision, are different from the solutions of variational problems. In the present paper, we will investigate the latter topic. We will study the functionals which are minimized by splines in affine spaces, but restrict the candidate curves to a given surface. The solution is not as simple and explicit as in the unrestricted case. However, we will find remarkable counterparts to known results. 1.1. Previous work Let us briefly sketch the literature dealing with splines on surfaces. We are focussing here only on curves which have some shape design handles or are generated as solutions of variational problems. Shoemake (1985) introduced spherical counterparts of Bézier curves. He extended de Casteljau’s algorithm to the sphere, replacing straight line segments by geodesic arcs and replacing ratios of Euclidean distances by ratios of geodesic distances. Using this on the 3-sphere in R4 , he generated the spherical component of rigid body motions and applied it to Computer Animation. Soon, it became apparent that the transition to the sphere has the following problem: One can use the spherical de Casteljau algorithm for evaluation of the curve, but one looses the powerful and important subdivision property. Thus, alternative constructions on the sphere have been proposed, see, e.g., (Kim et al., 1995; Jüttler and Wagner, 2002; Ramamoorthi and Barr, 1997). Again motivated by the application of motion design, Park and Ravani (1995) defined Bézier curves on Riemannian manifolds. Moreover, Sprott and Ravani (1997) extended B-spline algorithms to Lie groups and applied them to motion design. However, smoothness properties have not been proved. An algebraic approach to NURBS curves on quadrics, which is also rooted in line geometry and kinematics, has been proposed by Dietz et al. (1993). It has been shown that this method delivers a useful control structure and is suitable for interactive design (Pottmann and Wallner, 2001). The algebraic approach is extendable to low degree algebraic surfaces. In particular, it can be used for motion design (see e.g. (Jüttler and Wagner, 2002; Röschel, 1998)), but it is not suitable for curve design on general surfaces and manifolds. Among the few contributions dealing with subdivision schemes on manifolds we point to papers by L. Noakes (1998, 1999). His results concern C 1 smoothness of extensions of familiar subdivision schemes to manifolds. For example, he proved that the counterpart to the Lane–Riesenfeld algorithm for cubic B-spline curves on Riemannian manifolds yields a limit curve, which is differentiable and its derivative Lipschitz. Very recently, Wallner and Dyn (2005) found a technique to analyze nonlinear schemes based on certain ‘proximity’ conditions to known linear schemes. Their universally applicable results give remarkable extensions of known work; for example, it turns out that the Riemannian extension of cubic subdivision gives actually C 2 curves. We will not pursue subdivision in this paper, but concentrate on variational design. The ‘intersecting topic’, variational subdivision for curves on surfaces, has been computationally addressed, but not analyzed in (Hofer et al., 2003).
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
695
A number of papers dealing with variational design on surfaces is considering the sphere (Brunnett and Crouch, 1994; Brunnett et al., 1995), partially in view of the application to motion design (Park and Ravani, 1997; Ramamoorthi and Barr, 1997). An early contribution to variational curve design on more general surfaces is due to Noakes et al. (1989). The authors characterize the minimizers of an intrinsic geometric counterpart to the L2 norm of the second derivative. This is the integral of the squared covariant derivative of the first derivative with respect to arc length. For recent results on these intrinsic splines, see (Noakes, 2003; Camarinha et al., 2001). We will not pursue the intrinsic formulation in the present paper. The most closely related work to our contribution is the PhD thesis of H. Bohl (1999). It deals with the usual energy functionals used in CAGD, but restricts the curves to a given surface in R3 . Bohl proves the existence of a solution and gives computed examples on parametric surfaces, based on a quasi Newton optimization algorithm. A large part of his work deals with patch boundaries. If we consider on a surface the minimizers of the curve length, we obtain the well-studied geodesics. Their geometric properties have been investigated in classical differential geometry. Geodesics in a scaled arc length parameterization also arise as minimizers of the L2 norm of the first derivative. A variety of applications of geodesics has been described within Computer Vision and Image Processing (see e.g. (Osher and Paragios, 2003; Sapiro, 2001)). The present work has also been inspired by research on active contours (Blake and Isard, 1998), especially by work on active curves and geometric flows of curves on surfaces (see, e.g. (Cheng et al., 2002; Memoli and Sapiro, 2001; Osher and Paragios, 2003)). 1.2. Contributions and outline of the present paper We will study the minimizers of standard energy functionals of spline theory, which interpolate or approximate given data points, and are restricted to a given m-dimensional surface Φ ⊂ Rn . The restriction to Φ is the new aspect, and it is also the source of arising complications. However, we will find nice characterizations of the minimizers. In Section 2, we characterize the counterparts on surfaces to C 2 cubic splines, to splines in tension and to C 4 quintic splines. It will turn out that differentiability at the interpolation points is the familiar one, but the characterization of spline segments is slightly different. We give an example: Whereas the minimizers of the L2 norm of the second derivative have cubic segments (vanishing fourth derivative) in the unrestricted case, the corresponding splines on surfaces have segments with vanishing tangential component of the fourth derivative; we call such segments ‘tangentially cubic’. Hence, the differential equation c(4) = 0 changes to tpr c(4) = 0, with tpr denoting the tangential projection (orthogonal projection into the corresponding tangent space of Φ). Analogous changes are found in the characterizing differential equations for the other spline schemes. Section 3 is devoted to approximating curves, in particular to the counterparts of smoothing splines. In view of the importance of tangentially cubic curves, we give in Section 4 some explicit representations of such curves on special surfaces, namely certain cylinder surfaces. In particular, we address special tangentially cubic parametrizations of a curve; these are parametrizations whose fourth derivative is orthogonal to the curve. Computational algorithms are not the focus of the present paper. We address the computation just briefly in Section 5, where we also conclude the paper with examples and pointers to a number of interesting applications in Geometric Modeling, Computer Graphics, Robotics and Image Processing.
696
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
2. Interpolating spline curves on surfaces In view of the applications we have in mind, it is necessary to work on surfaces of arbitrary dimension and codimension. Thus, we consider an m-dimensional surface Φ in Euclidean Rn , m < n. Moreover, a sequence of points pi ∈ Φ, i = 1, . . . , N , and real numbers u1 < · · · < uN shall be given. We are seeking interpolating splines on the surface Φ. Sometimes we will also call Φ a manifold; if not stated otherwise, we work with a manifold without boundary. Thus, we have closed or unbounded surfaces. As we will point out in Section 5, the computation of energy-minimizing splines which respect patch boundaries or holes in a surface, can be reduced to the computation of a spline on (another) surface without boundary. 2.1. The counterparts to cubic splines Let us recall the situation, where we are not confined to a surface: Among all curves x(u) ⊂ Rn , whose first derivative is absolutely continuous, x˙ ∈ AC(I ), and whose second derivative satisfies x¨ ∈ L2 (I ) on I = [u1 , uN ], and which interpolate the given data, x(ui ) = pi , the unique minimizer of uN 2 (1) E2 (x) = x¨ (u) du, u1
is the interpolating C 2 cubic spline c(u). In the following we would like to extend this well-known result to the case where the admissible curves x(u) are restricted to the given surface Φ. We are not changing the functional E2 , whose interpretation requires an embedding space. We are considering the restriction to the surface as a constraint, rather than formulating the problem in terms of the intrinsic geometry of the manifold. As we will see later, this is a suitable formulation for the problems we would like to solve. Theorem 1. Consider real numbers u1 < · · · < uN and points p1 , . . . , pN on an m-dimensional C 4 surface Φ in Euclidean Rn . We let I = [u1 , uN ]. Then among all C 1 curves x : I → Φ ⊂ Rn , which interpolate the given data, i.e., x(ui ) = pi , i = 1, . . . , N , and whose restrictions to the intervals [ui , ui+1 ], i = 1, . . . , N − 1, are C 4 , a curve c which minimizes the functional E2 of Eq. (1) is C 2 and possesses segments c|[ui , ui+1 ], whose fourth derivative vectors are orthogonal to Φ. Moreover, at the end points p1 = c(u1 ) and pN = c(uN ) of the solution curve, the second derivative vector is orthogonal to Φ. Proof. If a solution curve c exists, the first variation of the energy functional must vanish there (Gelfand and Fomin, 1963). To express this condition, we consider neighboring curves x(u) ⊂ Φ, written as x(u, ε) = c(u) + h(u, ε).
(2)
For any fixed u˜ ∈ I , the curve h(u, ˜ ε) is a curve in Φ with h(u, ˜ 0) = 0. Its Taylor expansion at ε = 0 reads ε2 ˜ 0) + hεε (u, ˜ 0) + · · · , h(u, ˜ ε) = εhε (u, 2 ˜ 0) =: t(u) ˜ is a tangent vector of Φ at c(u). ˜ where the subscripts indicate differentiation. Note that hε (u, The displacement curves h to the given interpolation points vanish, h(ui , ε) = 0, for all ε. In particular,
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
697
Fig. 1. Displacement of the solution curve for deriving the first variation of the energy.
we have hε (ui , 0) = t(ui ) = 0. For the following, it is important to see that the mixed partial derivative vector hεu (ui , 0) = tu (ui ) is a tangent vector of Φ at c(ui ) = pi . Geometrically, this follows from the fact that the curve c(u) + t(u) is a curve on the ruled surface r(u, v) = c(u) + vt(u), which touches Φ along c. This curve passes through each interpolation point pi = c(ui ) + t(ui ), and thus its derivative vector cu (ui ) + tu (ui ) is a tangent vector of Φ. Together with the tangency of cu (ui ) this implies tangency of tu (ui ). Whatever field of displacement curves h(u, ε) we have chosen, the energy functional must assume a stationary value at ε = 0, d E2 x(u, ε) ε=0 = 0. (3) dε In view of 2 2 c¨ (u) + huu (u, ε) du = c¨ + εtuu + ε2 (· · ·) du, E2 x(u, ε) = I
I
this is equivalent to c¨ · tuu du = 0.
(4)
I
Integration by parts on each segment yields ui+1 ui+1 ui+1 (3) ui+1 c¨ · tuu du = c¨ · tu |ui − c · t|ui + c(4) · t du. ui
ui
We first note that the middle term vanishes, since t(ui ) = 0. Summing up over all intervals and denoting left and right derivatives at the knots ui by a subscript − and +, respectively, the condition for a solution c reads, N −1 (5) c¨ − (ui ) − c¨ + (ui ) · tu (ui ) + c¨ − (uN ) · tu (uN ) + c(4) · t du. 0 = −¨c+ (u1 ) · tu (u1 ) + i=2
I
698
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
Since t(u) is an arbitrary tangent vector field along c, and tu (ui ) are arbitrary tangent vectors at the given interpolation points, a solution curve c(u) must satisfy the following orthogonality conditions to the surface Φ. To formulate them, we denote by tpr the orthogonal projection of a vector at a point p ∈ Φ onto the tangent space τ of Φ at p, tpr c¨ + (u1 ) = tpr c¨ − (uN ) = 0, tpr c¨ − (ui ) − c¨ + (ui ) = 0, i = 2, . . . , N − 1,
(6)
tpr c(4) (u) = 0 on each interval [ui , ui+1 ].
(8)
(7)
Exactly these conditions are stated in Theorem 1. Eq. (6) means that the second derivative vector at the end points is orthogonal to the surface; this implies (but is not equivalent to) vanishing geodesic curvature at the end points. Eq. (7) shows that the tangential component of the second derivative is continuous along c; we call this behaviour tangentially C 2 or briefly T C 2 . It will be shown later that in view of the C 4 manifold Φ this implies C 2 . Finally, by Eq. (8), the fourth derivative vectors of the curve (which may be discontinuous at the interpolation points) are orthogonal to the manifold Φ. Vanishing fourth derivative characterizes a cubic curve; since here only the tangential component vanishes, we speak of a tangentially cubic curve in the following. To complete the proof, we show the following result: Lemma 2. On an m-dimensional C k manifold Φ ⊂ Rn , m < n, we consider a curve c(u), which is C k (k 2) everywhere except at a point c(u0 ). At c(u0 ), all derivatives up to order k have a continuous tangential component, i.e., (i) i = 1, . . . , k. (9) tpr c(i) − (u0 ) − c+ (u0 ) = 0, Then, the curve is also C k at c(u0 ). Proof. We may represent the manifold at least in a neighborhood of c(u0 ) as an intersection of n − m hypersurfaces, whose implicit representations shall be Fj (x) = 0,
j = 1, . . . , n − m.
Since the curve lies in all these hypersurfaces, we have the identity Fj (c(u)) = 0 in u. By differentiation, we obtain ∇Fj (c) · c˙ = 0, and for the ith derivative we find an expression of the form Gi−1 (c, . . . , c(i−1) ) + ∇Fj (c) · c(i) = 0.
(10)
Here, Gi−1 (c, . . . , c(i−1) ) abbreviates a function, which depends on Fj , up to differentiation order i, and on c, but only up to differentiation order i − 1. At c(u0 ) this holds for the left and right derivatives. The first derivative is tangential, and thus c is C 1 at u0 . This is the basis of an induction on the differentiation order i. Assuming continuity of the derivative c(i−1) at u0 , Eq. (10) expresses that the components ∇Fj (c) · c(i) of the ith derivative are continuous at u0 , since they equal Gi−1 (c, . . . , c(i−1) ), which is continuous by the induction hypothesis. Since the vectors ∇Fj (c(u0 )), j = 1, . . . , n − m, span the normal space of Φ at c(u0 ), we have shown that the normal component of the ith derivative is continuous. The tangential component is continuous by (9) and thus we have shown continuity of c(i) . This completes the proof of Lemma 2 and thus also the proof of Theorem 1. Note that the case k = 2 follows also immediately by Meusnier’s theorem. 2
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
699
Fig. 2. For a minimizer of (left) E1 and (right) E2 , the second respectively fourth derivative vectors (yellow) are orthogonal to the surface. (For interpretation of the references in color in this figure legend, the reader is referred to the web version of this article.)
The minimizers possess a characterization which is very similar to the familiar cubic C 2 splines. They are C 2 and tangentially cubic. However, the problem is nonlinear. Tangentially cubic curves can in general not be computed explicitly, but only by a numerical algorithm. Existence of the solution follows from the work of Bohl (1999) and Wallner (2004). Uniqueness cannot be expected, as will be clear from the following considerations. 2.2. Geodesics Let us replace the energy in (1) by the L2 norm of the first derivative, uN 2 E1 (x) = x˙ (u) du.
(11)
u1
Moreover, we just prescribe the two end points. Then we find with the same approach as in the proof of Theorem 1 that the curve’s second derivative vectors c¨ are orthogonal to Φ; we may say the curve is tangentially linear. In particular, c˙ and c¨ are orthogonal, and hence ˙c2 = const. This proves, that the curve is parameterized by a constant multiple of its arc length. Moreover, it shows that c¨ represents the principal normal, and orthogonality of the principal normal to Φ characterizes a geodesic. Thus we find the known result that the minimizers of E1 are geodesics in a scaled arc length parameterization. Note that the functionals we are considering are also optimizing the parameterization. Since even for geodesics we do not have simple results on uniqueness, the more involved case of splines will hardly allow us a characterization of situations with a unique solution.
700
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
Fig. 3. Spline curves c1 , ct , c2 on a surface Φ, interpolating the same points p1 , . . . , p4 and minimizing energies E1 , Et , and E2 , respectively.
2.3. Splines in tension on surfaces A linear combination of energies (11) and (1) leads in the unrestricted case to the well known splines in tension as minimizers (Schweikert, 1966). The counterpart on manifolds is characterized in the following theorem, whose proof can be omitted since it is completely analogous to that of Theorem 1. Theorem 3. Consider a sequence of data points, parameter values and admissible curves on a surface Φ as in Theorem 1. Then, a minimizer of the functional uN (12) Et (x) = (¨x2 + wx˙ 2 ) du, w = const > 0, u1
is a C 2 curve which satisfies tpr c(4) (u) − w¨c(u) = 0
(13)
on all segments. The end conditions are tpr c¨ + (u1 ) = tpr c¨ − (uN ) = 0. Increasing w brings the solution closer to the curve which minimizes E1 , i.e., the curve composed by geodesic segments between the interpolation points. Hence, the parameter w controls the tension of the curve, see Fig. 3. 2.4. Counterparts to higher order splines on surfaces It is known from spline theory, that the minimizers of the L2 norm of the third derivative, uN 2 E3 (x) = x(3) (u) du, u1
(14)
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
701
under interpolation conditions are quintic C 4 splines. It is not difficult to extend this result to manifolds. Theorem 4. Consider real numbers u1 < · · · < uN and points p1 , . . . , pN on an m-dimensional C 6 manifold Φ in Euclidean Rn . Then among all C 2 curves x : [u1 , uN ] → Φ ⊂ Rn , which interpolate the given data and whose restrictions to the intervals [ui , ui+1 ], i = 1, . . . , N − 1, are C 6 , a curve c which minimizes the L2 norm E3 of the third derivative is C 4 and possesses segments c|[ui , ui+1 ], whose sixth derivative vectors are orthogonal to Φ. Moreover, at the end points p1 = c(u1 ) and pN = c(uN ) of the solution curve, the third derivative vanishes and the fourth derivative vector is orthogonal to Φ. Proof. The proof resembles the one of Theorem 1. We use the same notation and just mention the few essential differences. A solution c(u) must satisfy c(3) · tuuu du = 0. (15) I
Integration by parts on each segment yields ui+1 ui+1 (3) (3) ui+1 (4) ui+1 (5) ui+1 c · tuuu du = c · tuu |ui − c · tu |ui + c · t|ui − c(6) · t du. ui
ui
From the four terms on the right hand side, the third one vanishes because of the interpolation conditions, expressed by t(ui ) = 0. Since t is an arbitrary tangent vector field along c, the last term vanishes iff tpr c(6) = 0. Now we can sum up over all intervals and collect the contributions from the first two terms c(3) · tuu and c(4) · tu . Since tu (ui ) is tangential, we find a T C 4 condition at the inner knots and orthogonality of the fourth derivative to the manifold at the end points. It is not hard to see that tuu (ui ) can be an arbitrary vector and this shows that a solution must have vanishing third derivative at the ends and it must be C 3 at the inner knots. Application of Lemma 2 shows that the curve is even C 4 at the inner knots, which completes the proof. 2 Remark 5. So far we always considered an open curve and natural end conditions. It is clear from the derivation of our results how to handle other cases. Let us discuss this at hand of the minimization of E3 . The segments are of course always tangentially quintic, i.e., tpr c(6) = 0. A closed spline is C 4 at all knots. Moreover, instead of working for an open curve with natural end conditions, we could prescribe first and second derivative vectors at the end points. The extension of the results to the minimization of the L2 norm of even higher derivatives, or combinations of L2 norms of different derivatives, is straightforward and thus it can be omitted.
3. Approximating spline curves on surfaces Reinsch (1967) relaxed the interpolation conditions x(ui ) − pi = 0 by adding the sum of squared interpolation errors as a penalty term to E2 , uN N 2 s x(ui ) − pi 2 , λ, µ > 0. ¨ x(u) du + µ (16) E2 (x) = λ u1
i=1
702
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
The minimizers, called smoothing splines, found a variety of applications in the analysis of observational data (Wahba, 1990). The choice of the smoothing parameter λ : µ in this method for data approximation is not a simple problem. Often one uses statistical methods for this critical task (Wahba, 1990). We should also mention that an analogous functional is used for freeform surface fitting in reverse engineering applications (Várady and Martin, 2002). We can use the same method as in the proof of Theorem 1 to derive the following result on the counterparts to cubic smoothing splines on surfaces. Theorem 6. Consider a sequence of data points, parameter values and differentiability conditions on an admissible curve on a surface Φ as in Theorem 1. Then, a minimizer of the smoothing spline functional E2s has all the properties of the minimizer of E2 in Theorem 1. However, the interpolation conditions c(ui ) = pi are replaced by the following transition conditions between segments, (3) = 0. (17) tpr λ c(3) + (ui ) − c− (ui ) + µ c(ui ) − pi In an analogous way, one can investigate the minimizers of smoothing spline counterparts Ets , E3s , . . . , to other energy functionals. For example, the case of E3s yields in modification of Theorem 4 curves, where the interpolation conditions are replaced by the smoothness conditions (5) = 0. (18) tpr λ c(5) − (ui ) − c+ (ui ) + µ c(ui ) − pi 4. Examples of tangentially cubic curves 4.1. Tangentially cubic curves on cylinder surfaces in R3 Let us consider a general cylinder surface Φ ⊂ R3 . We may choose a directrix curve l(u) = (x1 (u), x2 (u), 0) and rulings parallel to (0, 0, 1), i.e., a parameterization of Φ in the form Φ : x(u, v) = (x1 (u), x2 (u), v). For a tangentially cubic spline curve c(t) = (c1 (t), c2 (t), c3 (t)) ⊂ Φ, the fourth derivative c(4) must be orthogonal to Φ, and thus also orthogonal to the corresponding ruling of Φ. This requires a vanishing third component c3(4) (t) = 0 and shows that the function c3 is cubic, c3 (t) = a0 + a1 t + · · · + a3 t 3 . The other two coordinate functions (c1 (t), c2 (t)) have its fourth derivative orthogonal to a parameterization (c1 (t), c2 (t), 0) of the (planar) orthogonal cross section l of Φ. This is therefore a tangentially cubic parameterization of the cross section curve. Choosing the data points for a spline on a cross section, we find a result on energy minimizing parameterizations of a given curve. Corollary 7. Given a C 4 curve l in Rn , and a sequence of points pi with parameters ui , i = 1, . . . , N, on it. Then, among all C 1 , piecewise C 4 parameterizations of the curve which satisfy the interpolation conditions c(ui ) = pi , a minimizer of the L2 norm of the second derivative is tangentially cubic, i.e., a parameterization whose fourth derivative is orthogonal to l. Moreover, c is C 2 and the second derivative at the end points is orthogonal to l.
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
703
The derivation in Rn can be done with a cylinder surface in Rn+1 . We may also note that the case of a manifold of dimension 1 is actually included in the formulation and in the proof of Theorem 1. Since the study of splines on cylinders is reduced to tangentially cubic parameterizations of curves, we will study those in the following in more detail. 4.2. Tangentially cubic arc length parameterizations In view of the important role of arc length parameterizations for curves, we answer the question on which curves in Rn the arc length parameterization c(s) (or a scaled version of it), is tangentially cubic. Denoting the Frenet frame of c by e1 (s), . . . , en (s), and the curvatures by κ1 , . . . , κn−1 , we have c = e1 ,
e1 = κ1 e2 ,
e1 = κ1 e2 − κ 2 e1 + κ1 κ2 e3 .
Hence, the condition for a tangentially cubic parameterization becomes 0 = c(4) · c = e(3) 1 · e1 = −3κ1 κ1 .
Proposition 8. The arc length parameterization c(s) of a curve in Rn is tangentially cubic iff the curve possesses constant curvature κ1 . In the plane, we get only circles and straight lines. However, already in R3 this family is relatively rich, since the torsion τ = κ2 can be arbitrary. The most prominent examples are helices which even have constant torsion. 4.3. Tangentially cubic parameterizations of the circle; spline curves on a right circular cylinder Consider the unit circle x2 = 1 and its parameterization c(t) = cos φ(t), sin φ(t) . We want to determine φ(t) such that c(t) is a tangentially cubic parameterization of the circle. A simple computation shows that orthogonality of first and fourth derivative is equivalent to (19) φ (4) − 6φ˙ 2 φ¨ = 0. ˙ and thus have to solve the differential equation We set ψ(t) := φ(t) ˙ ψ (3) = 6ψ 2 ψ. We can immediately perform a first integration, ψ¨ = 2ψ 3 + A0 ,
(20)
with a constant A0 ∈ R. Since the particular solution ψ˙ = 0, i.e., the scaled arc length parameterization φ = at + b is already known from the previous considerations, we may assume ψ˙ = 0 from now on. Multiplying Eq. (20) by 2ψ˙ and integrating again, we arrive with A := 2A0 at ψ˙ 2 = ψ 4 + Aψ + B. Separation of the variables and integration leads to dψ . t +C =± ψ 4 + Aψ + B
(21)
704
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
Fig. 4. Examples of tangentially cubic curves with parameterization (23), (a) and (b), and parameterization (24), (c) and (d).
The integral on the right hand side requires Legendre normal integrals. We omit its discussion in the general case and just look at the special case A = B = 0. There, we obtain ψ = ±1/(t + C) and φ(t) = ln |C ± t| + D.
(22)
Thus, even in the very simple case of a right circular cylinder, just two very special families of tangentially cubic curves can be written down explicitly. On the unit cylinder x12 + x22 = 1, these are the curves derived from a scaled arc length parameterization c(t) = cos(at + b), sin(at + b), a0 + a1 t + · · · + a3 t 3 , (23) and the curves to the circle parameterization (22), c(t) = cos ln |C ± t| + D , sin ln |C ± t| + D , a0 + a1 t + · · · + a3 t 3 .
(24)
A few examples of curves (23) and (24) are depicted in Fig. 4. 4.4. Tangentially cubic parameterization by the direction angle of the curve tangent It turns out that nice explicit representations of tangentially cubic curves on cylinder surfaces arise as follows. We would like to determine those planar curves, whose parameterization c(φ) with respect to the direction angle φ of the curve tangent is tangentially cubic. Therefore, we write the family of tangent lines of such a curve in the form −x1 sin φ + x2 cos φ + h(φ) = 0.
(25)
Here, φ is the angle between the tangent, oriented by the vector (cos φ, sin φ), and the x1 -axis; h denotes the signed distance of the origin to the oriented tangent. The function h(φ) is known as support function. The curve parameterization c(φ) is computed as envelope of the family of tangents, i.e., from Eq. (25) and its derivative with respect to φ. We find c(φ) = (h sin φ + h˙ cos φ, −h cos φ + h˙ sin φ). (26) A tangentially cubic parameterization requires c˙ · c(4) = 0, which is equivalent to h(5) − 2h(3) − 3h˙ = 0.
(27)
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
The solution of this differential equation is, with real constants bi , √ h(φ) = b0 + b1 cos φ + b2 sin φ + b3 cosh( 3φ + b4 ).
705
(28)
To understand these curves, we note that an appropriate translation of the origin can eliminate the term b1 cos φ + b2 sin φ, and b4 just accounts for a rotation. Thus it is sufficient to study the case b1 = b2 = b4 = 0. Addition of the constant b0 to the support function means construction of the offset at signed distance b0 . The curve with support function √ (29) h(φ) = cosh 3φ, is a special instance of a so-called hypercycloid. Hypercycloids are a certain counterpart to the familiar epicycloids or hypocycloids, which are obtained as locus of a point on a circle, which rolls on another circle. A hypercycloid possesses such a kinematic generation, but the generating circles are not real: the moving circle has a complex radius and the fixed circle a purely imaginary radius. Hypercycloids arise for example as orthogonal projections of geodesics on a paraboloid of revolution, when projecting parallel to the rotational axis (see (Strubecker, 1969, pp. 227–232)). All other nontrivial solutions of our problem arise from the special hypercycloid with support function (29) by offsetting and a similarity transformation. The trivial solutions belong to b3 = 0 and are circles. On a cylinder surface with the hypercycloid to (29) as cross section, the curves √ √ √ cosh 3φ sin φ + 3 sinh 3φ cos φ √ √ √ c(φ) = − cosh 3φ cos φ + 3 sinh 3φ sin φ (30) a0 + a1 φ + a2 φ 2 + a3 φ 3 are tangentially cubic. Fig. 5 shows for φ ∈ [−2.5, 2.5] tangentially cubic curves (30) on a cylinder surface with the hypercycloid (29) as cross section, with (a) a0 = a1 = a3 = 0, a2 = 5 and (b) a0 = 0, a1 = 1, a2 = 2, a3 = 1.
Fig. 5. Tangentially cubic curves (30) on a cylinder surface with the hypercycloid (29) as cross section.
706
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
5. Computation based on an optimization algorithm The presented spline curves possess nice geometric characterizations. However, we found out in Section 4 that the problem is not even explicitly solvable on very simple surfaces. Hence, we have to use numerical algorithms based on a discretization. Geometrically this means that we replace the curve c(t) ⊂ Φ by a sufficiently dense polygon Pc with vertices ci ∈ Φ. Discretizing a quadratic functional, such as E1 , E2 , Et , results in a quadratic function in the coordinates of the vertices ci of the approximating polygon Pc . Since the polygon vertices are restricted to Φ, we end up with the constrained minimization of a quadratic function in a rather high dimensional space. From a geometric viewpoint, we have to construct the closest point p∗ on some manifold S to a given point p. Here, p represents the unrestricted (discretized) spline curve in Rn . Basically, one can use any algorithm for constrained optimization (Geiger and Kanzow, 2002). If the surface is given in parametric form, we can of course eliminate the constraints by the use of a polygon in the parameter domain. Then, an appropriate algorithm for unconstrained optimization over the parameter domain can be employed (Kelley, 1999). For example, H. Bohl (1999) uses the standard BFGS algorithm. However, several other efficient algorithms (quasi-Newton methods, Newton method with improvement of global behaviour, etc. (Kelley, 1999)) are available. We have developed a geometrically motivated optimization algorithm for the solution of the present optimization problem (Hofer and Pottmann, 2004). It extends and stabilizes the simple algorithm for variational curve design in (Hofer et al., 2003), and can also be interpreted within the Sobolev gradient method (Neuberber, 1997). This algorithm can be used for the computation of splines on parametric surfaces, level sets, triangle meshes and point set surfaces, see (Hofer and Pottmann, 2004). Fig. 3 depicts energy minimizing curves on a parametric surface. Fig. 6 shows energy minimizing curves on a point set surface and on a level set. All of them have been computed with the optimization algorithm presented in (Hofer and Pottmann, 2004).
Fig. 6. Curve minimizing E2 on a point set surface (left), and level set (right).
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
707
Another computational approach could directly use the characterization of the splines in a discretized way. For example, a spline minimizing E2 has 4th central difference vectors, which are orthogonal to Φ. Formulation of this property and the other characterizations gives a system of nonlinear equations, which can be solved with a Newton iteration. One will not apply the pure Newton algorithm, but use one of the known methods for improving the global convergence behaviour (Kelley, 1999). Remark 9. For applications in CAD, the surface Φ and the preimage of the desired spline curve in the parameter domain, have to be in NURBS form. To achieve a NURBS representation of the energy minimizing curve, we can use numerical optimization of the energy over a finite-dimensional linear space of B-spline curves in the parameter domain. Another approach separates the computation of the shape from the computation of a B-spline representation. One first computes the energy minimizing curve, represented by a polygon, with the numerical methods indicated above; then the resulting polygon vertices with their parameter values are approximated by one of the known schemes from CAGD. For the state of the art on the latter topic, we refer to (Renner and Weiss, 2004). 5.1. Respecting boundaries and holes Throughout the paper we assumed surfaces without boundary. This might appear as a severe limitation, especially in view of the frequent presence of boundaries in CAD related applications. Holes and boundaries can be treated as obstacles in a curve design problem. Fortunately, variational curve design in the presence of obstacles can be reformulated as curve design on another surface, which does not have a boundary. This method, a geometric version of the barrier method from constrained optimization (Geiger and Kanzow, 2002), is described and illustrated in (Hofer and Pottmann, 2004). 5.2. Applications Energy minimizing splines on surfaces possess a variety of applications. Partially, these go beyond standard CAGD applications. Some major applications are the following ones (Hofer and Pottmann, 2004): Curve design on the curved manifold, which represents the Lie group of Euclidean rigid body motions, leads to variational motion design for animation in Computer Graphics and motion planning problems in robotics. Applying the concept to so-called image manifolds, as introduced in (Kimmel et al., 2000), one can use splines for image sensitive drawing and for interactive image segmentation.
Acknowledgement This work was supported by the Austrian Science Fund (FWF) under grant P16002-N05.
References Blake, A., Isard, M., 1998. Active Contours. Springer, Berlin. Bohl, H., 1999. Kurven minimaler Energie auf getrimmten Flächen. PhD thesis, Univ. Stuttgart, Germany. Brunnett, G., Crouch, P.E., 1994. Elastic curves on the sphere. Adv. Comp. Math. 2, 23–40.
708
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
Brunnett, G., Crouch, P.E., Silva Leite, F., 1995. Spline elements on spheres. In: Daehlen, M., Lyche, T., Schumaker, L. (Eds.), Mathematical Methods for Curves and Surfaces. Vanderbilt Univ. Press, pp. 49–54. Brunnett, G., Hagen, H., Santarelli, P., 1993. Variational design of curves and surfaces. Surv. Math. Ind. 3, 1–27. Camarinha, M., Silva Leite, F., Crouch, P.S., 2001. On the geometry of Riemannian cubic polynomials. Differential Geom. Appl. 15, 107–135. Cheng, B.-T., Burchard, P., Merriman, B., Osher, S., 2002. Motion of curves constrained on surfaces using a level set approach. J. Comput. Phys. 175 (2), 604–644. Dietz, R., Hoschek, J., Jüttler, B., 1993. An algebraic approach to curves and surfaces on the sphere and on other quadrics. Computer Aided Geometric Design 10, 211–229. Geiger, C., Kanzow, C., 2002. Theorie und Numerik restringierter Optimierungsaufgaben. Springer, Heidelberg. Gelfand, I.M., Fomin, S.V., 1963. Calculus of Variations. Prentice-Hall, Englewood Cliffs, NJ. Hofer, M., Pottmann, H., Ravani, B., 2003. Geometric design of motions constrained by a contacting surface pair. Computer Aided Geometric Design 20, 523–547. Hofer, M., Pottmann, H., 2004. Energy-minimizing splines in manifolds. In: Proceedings of ACM SIGGRAPH 2004, Trans. Graph. 23 (3), 284–293. Hoschek, J., Lasser, D., 1993. Fundamentals of Computer Aided Geometric Design. A.K. Peters, Wellesley, MA. Jüttler, B., Wagner, M., 2002. Kinematics and animation. In: Farin, G., Kim, M.-S., Hoschek, J. (Eds.), Handbook of Computer Aided Geometric Design. North Holland, Amsterdam, pp. 723–748. Kelley, C.T., 1999. Iterative Methods for Optimization. SIAM, Philadelphia. Kim, M.-J., Kim, M.-S., Shin, S., 1995. A general construction scheme for unit quaternion curves with simple high order derivatives. In: SIGGRAPH ’95, Comput. Graph. 29, 369–376. Kimmel, R., Malladi, R., Sochen, N., 2000. Images as embedded maps and minimal surfaces: movies, color, texture and volumetric medical images. Internat. J. Computer Vision 39, 111–129. Memoli, F., Sapiro, G., 2001. Fast computation of weighted distance functions and geodesics on implicit hyper-surfaces. J. Comp. Phys. 173, 730–764. Moreton, H., Sequin, C., 1993. Scale-invariant minimum-cost curves: fair and robust design implements. Comput. Graph. Forum 12, 473–484. Neuberber, J.W., 1997. Sobolev Gradients and Differential Equations. Springer, Heidelberg. Noakes, L., Heinzinger, G., Paden, B., 1989. Cubic splines on curved spaces. IMA J. Math. Control Inform. 6, 465–473. Noakes, L., 1998. Non-linear corner cutting. Adv. Comp. Math. 8, 165–177. Noakes, L., 1999. Accelerations of Riemannian quadratics. Proc. Amer. Math. Soc. 127, 1827–1836. Noakes, L., 2003. Null cubics and Lie quadratics. J. Math. Phys. 44 (3), 1436–1448. Osher, S., Paragios, N. (Eds.), 2003. Geometric Level Set Methods. Springer, New York. Park, F.C., Ravani, B., 1995. Bézier curves on Riemannian manifolds and Lie groups with kinematic applications. ASME J. Mech. Design 115, 36–40. Park, F.C., Ravani, B., 1997. Smooth invariant interpolation of rotations. ACM Trans. Graph. 16, 277–295. Pottmann, H., Wallner, J., 2001. Computational Line Geometry. Springer, Heidelberg. Ramamoorthi, R., Barr, A., 1997. Fast construction of accurate quaternion splines. In: SIGGRAPH ’97, Comput. Graph. 31, 287–292. Reinsch, C.H., 1967. Smoothing by spline functions. Numer. Math. 10, 177–183. Renner, G., Weiss, V., 2004. Exact and approximate computation of B-spline curves on surfaces. Computer-Aided Design 36, 351–362. Röschel, O., 1998. Rational motion design—a survey. Computer-Aided Design 30, 169–178. Sapiro, G., 2001. Geometric Partial Differential Equations and Image Analysis. Cambridge Univ. Press, Cambridge. Schweikert, D., 1966. An interpolation curve using a spline in tension. J. Math. Phys. 45, 312–317. Shoemake, K., 1985. Animating rotations with quaternion curves. In: SIGGRAPH ’85, Comput. Graph. 19, 245–254. Sprott, K., Ravani, B., 1997. Ruled surfaces, Lie groups and mesh generation. In: Proc. ASME Design Eng. Techn. Conf. No DETC97/DAC-3966. Strubecker, K., 1969. Differentialgeometrie II: Theorie der Flächenmetrik. Walter de Gruyter, Berlin. Várady, T., Martin, R., 2002. Reverse Engineering. In: Farin, G., Hoschek, J., Kim, M.-S. (Eds.), Handbook of Computer Aided Geometric Design. Elsevier, Amsterdam, pp. 651–681. Wahba, G., 1990. Spline Models for Observational Data. SIAM, Philadelphia.
H. Pottmann, M. Hofer / Computer Aided Geometric Design 22 (2005) 693–709
709
Wallner, J., 2004. Existence of set-interpolating and energy-minimizing curves. Computer Aided Geometric Design 21 (9), 837–857. Wallner, J., Dyn, N., 2005. Convergence and C 1 analysis of subdivision schemes on manifolds by proximity. Computer Aided Geometric Design 22. In press.