C1 and C2 -continuous polynomial parametric Lp splines (p⩾1)

Computer Aided Geometric Design 24 (2007) 373–394 www.elsevier.com/locate/cagd

C 1 and C 2 -continuous polynomial parametric Lp splines (p  1) P. Auquiert a , O. Gibaru b,∗ , E. Nyiri b a (LAMAV) Université de Valenciennes et du Hainaut-Cambrésis, Le Mont Houy, F-59313 Valenciennes cedex 9, France b (L2MA) C.E.R. Ecole Nationale Supérieure d’Arts et Métiers de Lille, 8 Bld Louis XIV, F-59046 Lille cedex, France

Received 20 October 2006; received in revised form 1 February 2007; accepted 16 April 2007 Available online 17 May 2007

Abstract We introduce Lp (1  p  ∞) piecewise polynomial parametric splines of degree 3 or 5 which smoothly interpolate data points. The coefficients of these polynomial splines are calculated by minimizing the Lp norm of the second derivative. It is demonstrated that the C 1 -continuous cubic (resp. G1 -continuous cubic and C 2 -continuous quintic) Lp polynomial splines are unique for 1 < p < ∞. The focus is mainly on L1 polynomial splines which preserve the shape of data even with abrupt changes in direction or spacing. When there are several candidates we introduce a tension parameter which selects from the set of solutions the smallest in norm. The optimisation uses a nonlinear functional; therefore we discretize it so as to obtain a linear problem and we give an upper bound to the error. Computational examples using a primal affine (interior point) algorithm are presented. We will also show how interesting it can be to be able to change the partition associated to the data points in order to obtain convexity properties when possible. Moreover, we will show that the solution is not unvarying by rotation in the L1 -case and we propose an alternate approach with a local change of coordinates on each segment [Pi , Pi+1 ]. We show that this method can be applied to any kind of norm for the minimization problem. © 2007 Elsevier B.V. All rights reserved. MSC: 65D17; 14B05; 65D10; 65D05; 68U05 Keywords: Interpolation; Lp polynomial parametric splines; L1 polynomial parametric splines

1. Introduction Spline curves have been intensively investigated (for instance see (De Boor, 1978; Farin, 1993)). They have various useful properties making them a standard in schemes in computer aided geometric design and in NC tool path generation. In this article we have chosen to study Lp polynomial splines because they are more easily exported on computer. Let (a = x0 < x1 < · · · < xn = b) be a strictly monotonic partition of a finite real interval [a, b]. In (Lavery, 2000), Lavery defines cubic Lp splines which interpolate the data {xi , yi = f (xi )}i∈{0,...,n} with C 1 piecewise cubic polynomial, the coefficients of which are calculated by minimizing the Lp -norm of their second derivative. Lavery shows that cubic L2 splines, which coincide with the conventional C 2 -cubic splines, and cubic L∞ splines do not preserve * Corresponding author.

E-mail address: [email protected] (O. Gibaru). 0167-8396/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.cagd.2007.04.007

374

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

shape well. In contrast, cubic L1 splines provide C 1 -smooth continuity, shape preserving, multiscale interpolation of arbitrary data, including data with abrupt changes along the y-axis in spacing and magnitude, with no need for monotonicity or convexity constraints, node adjustment or other user input. When the given adjacent points show abrupt changes, the C 1 -continuous cubic L1 splines avoid the Gibbs phenomena (see Figs. 1 and 3) in contrast to the C 1 -continuous cubic L2 splines which are actually C 2 . Indeed, Zang and Martin (1997) demonstrate that for a function f with isolated discontinuity, the C 2 -cubic spline which interpolates it at {xi , yi = f (xi )}i∈{0,...,n} , oscillates near the discontinuous points with a maximal overshoot about 4% in the limit when h = maxi∈{0,...,n−1} |xi+1 − xi | tends to zero. They also demonstrate that the oscillation decreases exponentially in the region away from the discontinuity (see Fig. 1). Moreover Saff and Tashev (1999) show that the overshoot and oscillation of the best Lp -polygonal line approximation over [−1, 1] with knots at { ni }ni=−n of functions having a jump discontinuity, decrease as p tends to 1 and disappear completely when p = 1. In (Cheng et al., 2002, 2005), Cheng, Fang and Lavery develop a geometric programming approach for the univariate cubic L1 splines, which gives some results on the ability of cubic L1 spline to preserve the shape of nonparametric data. They demonstrate that when four consecutive data points lie on a straight line, the cubic L1 spline is linear in the interval between the second and the third data points. They also demonstrate that the cubic L1 spline of convex/concave data preserves convexity/concavity if the first divided differences of the data do not increase/decrease too rapidly. We propose to extend this work to RN (N  2), where Lp splines of degree 3 or 5 are defined such that they smoothly interpolate the data points. The interpolation problem we are tackling is as follows. Let {Pi ∈ RN }i=0,...,n be a set of n + 1 points and let a = u0 < u1 < · · · < un = b be an arbitrary and strictly monotonic partition of the finite real interval [a, b]. We define the Lp spline γ of degree k (k ∈ N∗ ) to be the piecewise Bézier curves of degree k such that γ (ui ) = Pi , for i = 0, . . . , n. The coefficients of these polynomial splines are calculated by minimizing the Lp -norm of the second derivative of these piecewise Bézier curves. The paper is organised as follows. After introducing the basic notation for the interpolation by Lp polynomial splines of degree k, 1  p  ∞, we will see in detail interpolation by C 1 -continuous cubic, G1 -continuous cubic and C 2 -continuous quintic Lp polynomial splines in RN . We will demonstrate that they are unique for 1 < p < ∞. The focus is mainly on L1 polynomial splines which preserve the shape of the data points with abrupt changes along an axis direction. When there are several candidates we will introduce a tension parameter which selects the smallest solution in norm from the set of solutions. In the L1 -case, minimization gives a nonlinear functional, so we will discretize it by the midpoint rule with equal subintervals. The approximation thus obtained is a linear problem that can be solved by a primal affine algorithm (Vanderbei et al., 1986). We will give an upper bound for the error. Examples in the plane are given to illustrate the case of C 1 -continuous cubic, G1 -continuous cubic and C 2 -continuous quintic L1 polynomial splines and we compare them with υ-splines and τ -splines. We show that the ability to modify the partition ui is very important so as to search a convexity property according to an axis direction. Moreover, we show that whatever the order of the derivative used in the minimization norm (see (Lavery, 2006)) the solutions are not unvarying by rotation in the L1 -case. To do so, we create an alternate approach with a local change of coordinates on each segment [Pi , Pi+1 ] (see Figs. 14, 15a, 15b and 15c). 2. Interpolation by parametric Lp splines of degree k Definition 1. Let us give a set of n + 1 points Pi ∈ RN (N  2) for i = 0, . . . , n. Let a = u0 < u1 < · · · < un = b be an arbitrary and strictly monotonic partition of the finite real interval [a, b]. Then we define the Lp spline γ of degree k (k ∈ N∗ ) to be the piecewise Bézier curve of degree k such that:   i for u ∈ [ui , ui+1 ], Qj ∈ RN (a) γ (u) = Bjk [ui , ui+1 ](u)Qij ,   i+1 −u k−j u−ui j where Bjk [ui , ui+1 ](u) = jk uui+1 −ui ui+1 −ui , (b) γ (ui ) = Pi ,

for i = 0, . . . , n.

(c) spline γ minimizes  2 p d γ     du2  du when 1  p < ∞ p

(2.1)

(2.2)

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

or

 2  d γ    0in−1 ui uui+1  du2 ∞ max

sup

when p = ∞.

375

(2.3)

3. C 1 -continuous cubic Lp splines (1  p  +∞) (k = 3) In that case we assume that we have continuous first derivatives at each node ui , i.e., lim

u→u− i

dγ dγ (u) = lim (u)(= Ti ) + du du u→ui

for i = 1, . . . , n − 1,

(3.1)

moreover, we note: T0 =

dγ (u0 ) du

and Tn =

dγ (un ). du

(3.2)

Let T be the set of the n + 1 vectors Ti ∈ RN denoted by T = (T0 , T1 , . . . , Tn ). Theorem 2. The C 1 -continuous cubic Lp splines exist for 1  p  +∞. Proof. In the interval [ui , ui+1 [ a C 1 -continuous cubic Lp spline γ can be expressed by γ (u) =

3 

Bj3 [ui , ui+1 ](u)Qij

for i = 0, . . . , n − 1.

j =0

The constraints (2.1), (3.1) and (3.2) give: ui+1 − ui ui+1 − ui Ti and Qi2 = Pi+1 − Ti+1 . 3 3 It leads to the Hermite or Fergusson formula for a cubic spline curve: Qi0 = Pi ,

Qi3 = Pi+1 ,

γ (u) = Pi + (u − ui )Ti +

Qi1 = Pi +

(u − ui )2 (u − ui )3 (3Pi − 2Ti − Ti+1 ) + (Ti+1 + Ti − 2Pi ) hi h2i

for u ∈ [ui , ui+1 [ and i = 0, . . . , n − 1, where hi = ui+1 − ui and Pi = Pi+1hi−Pi . The second derivative can be calculated at all points of the open interval ]ui , ui+1 [ by     d 2γ 6 1 hi T + u − u (Ti+1 + Ti − 2Pi ). (u) = − T − i+1 i i hi 2 du2 h2i For 1  p < ∞ we substitute on each ]ui , ui+1 [ the suitable value of changing the variable of integration on the interval [ui , ui+1 ] from u to

d2γ (u) (see du2 h u−ui − 2i t= hi

(3.3)

(3.3)) into the functional (2.2), gives the following functional

1

E(T )

n−1  2  i=0

− 12

 1  (1 + 6t)Ti+1 + (−1 + 6t)Ti − 12tPi p dt.

(3.4)

p

p−1 hi

Therefore, we assume that the coordinates of a vector Ti are written Ti,k for k = 1, . . . , N . For instance Ti = (Ti,1 , . . . , Ti,N ). The functional (3.4) can be written as follows E(T ) = N k=1 Ek (T ) where 1

Ek (T ) =

n−1  2  i=0

− 12

p 1

(1 + 6t)Ti+1,k + (−1 + 6t)Ti,k − 12tPi,k dt

p−1 hi

for 1  k  N.

376

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

As we can see, each Ek (T ) is positive, then we minimise each Ek (T ) so as to obtain the minimal value of E(T ). Ek (T ) =

n−1 

1

ϕ p−1 i,k h i=0 i

where for i = 0, . . . , n − 1,

1

2 ϕi,k (T ) =



(1 + 6t)Ti+1,k + (−1 + 6t)Ti,k − 12tPi,k p dt.

− 12

The functions ϕi,k are convex, continuous and bounded below by 0. This implies that Ek is convex, continuous and bounded below by 0. Moreover limT 2 →+∞ Ek (T ) = +∞ which, via a classical result ((Brezis, 1999) Corollary III.20) allows us to conclude that the minimum of E(T ) exists for all p ∈ [1, +∞[. For p = ∞ we substitute (3.3) into the functional (2.3). Changing the variable from u to t = functional:  1 (1 + 6t)Ti+1 + (−1 + 6t)Ti − 12tPi  F (T ) = max sup ∞ h 0in−1 1 − t 1 i 2

h u−ui − 2i hi

gives the (3.5)

2

which can be written: F (T ) =

max

0in−1

sup − 12 t 21

1 max φi,k,t (T ) hi 1kN

where for k = 1, . . . , N , i = 0, . . . , n − 1 and t ∈ [− 12 , 12 ] the functions φi,k,t : (RN )n+1 → R are defined by



φi,k,t (T ) = (1 + 6t)Ti+1,k + (−1 + 6t)Ti,k − 12tPi,k . Functions φi,k,t are continuous and convex so functions φi,t ∞ = max1kN φi,k,t are lower semi-continuous, convex and moreover bounded below by 0. Consequently the functions φi,t ∞ , h1i φi,t ∞ , sup− 1 t 1 h1i φi,t ∞ 2

2

and F (T ) = max0in−1 sup− 1 t 1 h1i φi,t ∞ are lower semi-continuous ((Choquet, 1969) Theorem 8-6), convex 2 2 and bounded below by 0. Moreover limT 2 →+∞ F (T ) = +∞ and via ((Brezis, 1999) Corollary III.20) we conclude that the minimum of F (T ) exists. Hence, for any given p, there exists a set of vectors Ti that minimizes E(T ) and F (T ), which implies that the corresponding C 1 -continuous cubic Lp splines exist for all p ∈ [1, +∞]. 2 Theorem 3. The C 1 -continuous cubic Lp splines are unique for 1 < p < +∞. Proof. When 1 < p < +∞ functional (3.4) is strictly convex, hence, there exists a unique set of vectors Ti that minimizes this expression, which implies that the corresponding cubic Lp splines are unique. 2 Remark 4. As in (Lavery, 2000) (Figs. 4 and 5), Fig. 5 shows that the L1 parametric splines are not necessarily unique. For p = 1 (respectively p = ∞), functional (3.4) (respectively (3.5)) are not necessarily strictly convex. To show this, let Θ be defined by Θ = max1kN max0in−1 Pi,k and the hi be such that hi = h for i = 0, . . . , n − 1. Let α be an arbitrary value in ]0, 1[ and B (respectively C) be a vector of (RN )n+1 such that its coordinate values are defined by Bi,k = Θ + 1 (respectively Ci,k = Θ + 2) for i = 0, . . . , n and k = 1, . . . , N . For p = 1 after some calculations we obtain that: Ek (αB + (1 − α)C) = 24

n−1  

1 2

(Θ − Pi,k + 2 − α)t dt,

i=0 0

Ek (B) = 24

n−1   i=0 0

1 2

(Θ − Pi,k + 1)t dt,

Ek (C) = 24

n−1   i=0 0

1 2

(Θ − Pi,k + 2)t dt.

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

377

Hence, αEk (B) + (1 − α)Ek (C) = Ek (αB + (1 − α)C), which allows us to conclude. For p = +∞, by setting i = min1kN Pi,k for i = 0, . . . , n − 1 and Θ  = max0in−1 i we obtain after some calculations that: 6 (Θ − α + 2 − Θ  ), h 6 6 F (B) = (Θ + 1 − Θ  ), F (C) = (Θ + 2 − Θ  ). h h F (αB + (1 − α)C) =

Hence, αF (B) + (1 − α)F (C) = F (αB + (1 − α)C), which allows us to conclude. The choice of p influences the form of the interpolating curves. In the following examples, the data points show abrupt changes. Fig. 1 shows that L1 splines do not show overshoot in contrast to the L2 splines (see (Lavery, 2000) and (Zang and Martin, 1997)). So we will now focus on the L1 case and will give a practical method to generate such splines. The curvatures of the C 1 -continuous cubic L1 splines of Fig. 1 are not continuous, they show abrupt changes induced by discontinuity at the Pi points (see calculus in Remark 5 for Fig. 4).

Fig. 1. Interpolation of points with cubic splines using a uniform parametrization. - - - - C 1 -continuous cubic L2 spline, —— C 1 -continuous cubic i , 0)i=6 , ( 1 , 1 ), ( i , 1)i=14 . L1 spline. Data points on the left: (0, 0), ( 16 , 0), ( 13 , 0), ( 12 , 12 ), ( 23 , 1), ( 56 , 1), (1, 1). Data points on the right: ( 14 i=0 2 2 i=8 14

Fig. 2. The curvatures of Fig. 1 curves.

Fig. 3. Interpolation of points with cubic splines using a uniform parametrization. - - - C 1 -continuous cubic L2 spline, —— C 1 -continuous cubic 3 7 L1 spline. Data points: (0, 0), ( 12 , 0), (1, 0), (1, 1), ( 32 , 1), (2, 1), ( 13 4 , − 2 ), ( 2 , 0), (4, 0), (5, 0).

378

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

Fig. 4. The curvatures of the Fig. 3 curves.

Fig. 5. - - - the Ti are respectively: (9.898, 0), (9.898, 0), (7, 0), (7, 7), (7, 7), (7, 0), (9.898, 0), (9.898, 0). —– the Ti are respectively: (9.898, 0), (9.898, 0), (7, 7), (7, 7), (7, 7), (7, 7), (9.898, 0), (9.898, 0).

The curvature of the C 1 -continuous cubic L1 spline is not continuous. It shows discontinuity at the Pi points shown by the calculus in Remark 5 below. Remark 5. The C 1 -continuous cubic L2 splines are always C 2 , but the C 1 -continuous cubic L1 splines are generally 2 only C 1 . A C 1 -continuous cubic L1 spline is C 2 if h2i (3Pi − 2Ti − Ti+1 ) = hi−1 (−3Pi−1 + Ti−1 + 2Ti ) for i = 1, . . . , n − 1, which is not satisfied in general. In the above example (see Fig. 3) we calculate

Ti+1 ) −

2

hi−1 (−3Pi−1

2 hi (3Pi

− 2Ti −

+ Ti−1 + 2Ti ) for each value of i. The values are given in the following table: Number of the point

x

y

1 2 3 4 5 6 7 8

6.225383 × 10−8

−4.60425 × 10−9 486 486 −2.80464 × 10−9 −1215 −486 729 −2.99167 × 10−9

−243 −243 3.63898 × 10−8 405 445.5 −5.986988 × 10−7 −121.5

3.1. C 1 -continuous cubic L1 splines The functional (3.4) can be written as follows: E(T ) =

N

k=1 Ek (T )

where

1 2

n−1  



6t (Ti+1,k + Ti,k − 2Pi,k ) + Ti+1,k − Ti,k dt Ek (T ) =

for 1  k  N.

i=0 1 −2

As we can see Ek (T ) is a nonlinear nonsmooth convex function. Moreover E(T ) is not strictly convex, so its minima need not be unique. We can have several solutions, for example with the data set (0, 0), (1.414, 0), (2.818, 0), (3.818, 1), (4.818, 2), (5.818, 3), (7.242, 3), (8.656, 3), the two curves shown in Fig. 5 are candidates.

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

379

To reduce the set of solutions we propose as in (Lavery, 2000) to select derivative vectors Ti which are the shortest (in norm) possible. Consequently, we create a C 1 -continuous cubic L1 spline minimizing the new functional G(T ) = Here ε is a strictly positive real that acts as a tension parameter. This functional can be written as E(T ) + εT 1 . follows G(T ) = N k=1 Gk (T ) where Gk (T ) = Ek (T ) + ε

n 

|Ti,k |.

(3.6)

i=0

A discretization of the integral in (3.6) by the midpoint rule with d equal subintervals gives the following approximation of Gk (T ): k (T ) = G

n−1  d−1 

1



(−2d + 6j + 3)T + (−4d + 6j + 3)T − (−6d + 12j + 6)P i+1,k i,k i,k d2 i=0 j =0 n 

+ε

|Ti,k |.

(3.7)

i=0

Proposition 6. An upper bound for the discretization error of functional (3.6) is given by n−1

 1 3|Ti+1,k + Ti,k − 2Pi,k |, if |Ti+1,k − Ti,k |< 3|Ti+1,k + Ti,k − 2Pi,k |

Gk (T ) − G k (T )  d 2 0, otherwise. i=0

k (T )| = O( 12 ) for all T ∈ (RN )n . Hence |Gk (T ) − G d To demonstrate this result we use the following trivial lemma: Lemma 7. Let f : [c, d] → R be an affine function such that f (t) = at + b.

d (a) If |f (t)| > 0 for all t ∈ ]c, d[ then c |f (t)| dt = (d − c)|f ( d+c 2 )|, (b) If there exists ξ ∈ ]c, d[ such that f (ξ ) = 0 then

d

















d + c



2 a





 (d − c) .

f (t) dt − (d − c) f



2 2 c

(c) Moreover if c = − 12 and d = 12 , |f (t)| > 0 for all t ∈ ]− 12 , 12 [ if and only if 2|b|  |a|. Proof of Proposition 6. If |Ti+1,k − Ti,k |  3|Ti+1,k + Ti,k − 2Pi,k | then from Lemma 7 there is no error between

12 1 |f (t)| dt and its discretization. Otherwise, there exists an index l (0  l  d − 1) such that f takes the value 0 in −2

the interval ]− 12 + dl , − 12 +

l+1 d [.

The statement of Proposition 6 follows directly.

2

k (T ) is a linear problem which can be solved by the primal affine algorithm For k = 1, . . . , N , minimizing G developed by Vanderbei, Meketon and Freedman and outlined in (Vanderbei et al., 1986) and (Lavery, 2001). The k (T ) can be written functionals G     T0,k     − Ck  . A.   Tn,k 1 where the matrix A = (ai,j ) 1id×n+n+1 and the vector Ck = (ck,i )1id×n+n+1 are defined by: 1j n+1

for 0  i  d − 1 and for 0  j  n − 1 −4d + 6i + 3 −2d + 6i + 3 , aj ×d+i+1,j +2 = , aj ×d+i+1,j +1 = d2 d2

380

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

(−6d + 12i + 6)Pj,k , d2 for 1  i  n + 1 and for 1  j  n + 1

ck,j ×d+i+1 =

if i = j then an×d+j,i = 0

else an×d+j,j = ε

ck,n×d+j = 0. 4. G1 -continuous cubic Lp splines (1  p  +∞) (k = 3) Let T be a set of n + 1 given vectors Ti ∈ RN denoted by T = (T0 , T1 , . . . , Tn ). Let ε be an arbitrary value of ]0, +∞[. Let α = (α1 , α2 , . . . , α2n ) be a set of 2n real numbers which are bounded below by ε. We now assume that we have derivatives at each node ui such that lim

dγ (u) = α2i Ti du

and

lim

dγ (u) = α1 T0 du

and

u→u− i u→u+ 0

lim

dγ (u) = α2i+1 Ti du

lim

dγ (u) = α2n Tn . du

u→u+ i u→u− n

for i = 1, . . . , n − 1,

Theorem 8. Let us give a set of n + 1 vectors Ti ∈ RN and let ε be a positive real. Then the G1 -continuous cubic Lp splines exist for all 1  p  +∞. Proof. In the interval [ui , ui+1 [ a G1 -continuous cubic Lp spline γ can be expressed by γ (u) =

3 

Bj3 [ui , ui+1 ](u)Qij .

j =0

The previous constraints give: Qi0 = Pi ,

Qi3 = Pi+1 ,

Qi2 = Pi+1 −

Qi1 = Pi +

(ui+1 − ui ) α2i+2 Ti+1 . 3

(ui+1 − ui ) α2i+1 Ti 3

and

It follows that γi (u) = Pi + (u − ui )α2i+1 Ti + +

(u − ui )2 (3Pi − 2α2i+1 Ti − α2i+2 Ti+1 ) hi

(u − ui )3 (α2i+2 Ti+1 + α2i+1 Ti − 2Pi ) h2i

for u ∈ [ui , ui+1 [ and i = 0, . . . , n − 1, where hi = ui+1 − ui and Pi = Pi+1hi−Pi . The second derivative can be calculated at all points of the open interval ]ui , ui+1 [ by   d 2 γi 6 hi 1 u − u (α2i+2 Ti+1 + α2i+1 Ti − 2Pi ) + (α2i+2 Ti+1 − α2i+1 Ti ). (u) = − i 2 2 2 hi du hi For 1  p < ∞ we substitute on each ]ui , ui+1 [ the expression for Changing the variable of integration on the interval [ui , ui+1 ] from u to

d2γ (u) given in (4.1) into functional du2 h u−ui − 2i t= gives the functional hi

(4.1) (2.2).

1

E(α, T ) =

n−1  2  i=0

− 12

 1  (1 + 6t)α2i+2 Ti+1 + (−1 + 6t)α2i+1 Ti − 12tPi p dt

p−1 hi

where α = (α1 , . . . , α2n ).

p

(4.2)

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

Now we set the vectors Ti for i = 0, . . . , n, then E(α) =

N k=1

n−1

1 i=0 hp−1 ψi,k i

381

where

1

2 ψi,k (α) =



(1 + 6t)α2i+2 Ti+1,k + (−1 + 6t)α2i+1 Ti,k − 12tPi,k p dt

− 12

for i = 0, . . . , n − 1 and k = 1, . . . , N . The functions ψi,k are convex, continuous and bounded below by 0. Consequently,

1 p−1

hi

ψi,k ,

n−1

1 i=0 hp−1 ψi,k , E(α) i

are convex, continuous and bounded below by 0. Moreover limα2 →+∞ E(α) = +∞ which allows us to conclude that the minimum of E(α) exists for all p ∈ [1, +∞[. For p = ∞ we substitute (4.1) into functional (2.3). By changing the variable from u to t = setting the vectors Ti for i = 0, . . . , n we have

1

(1 + 6t)α2i+2 Ti+1,k + (−1 + 6t)α2i+1 Ti,k − 12tPi,k

F (α) = max sup max 0in−1 1 1kN hi − t 1 2

h u−ui − 2i hi

and now (4.3)

2

or F (α) =

max

0in−1

sup − 12 t 21

1 max χi,k,t (α) hi 1kN

where the functions χi,k,t : R2n+1 → R are defined by



χi,k,t (α) = (1 + 6t)α2i+2 Ti+1,k + (−1 + 6t)α2i+1 Ti,k − 12tPi,k

for k = 1, . . . , N , i = 0, . . . , n − 1 and t ∈ [− 12 , 12 ]. Functions χi,k,t are continuous and convex so functions χi,t ∞ = max1kN φi,k,t are continuous, convex and moreover bounded below by 0. Consequently the functions χi,t ∞ , h1i χi,t ∞ , sup− 1 t 1 h1i χi,t ∞ and 2

F (α) =

max

0in−1

sup − 12 t 21

2

1 χi,t ∞ hi

are lower semi-continuous ((Choquet, 1969) Theorem 8-6), convex and bounded below by 0. Moreover lim

α2 →+∞

F (α) = +∞

and via ((Brezis, 1999) Corollary III.20) we conclude that the minimum of F (α) exists. Hence, for any given p, there exists a set (α1 , . . . , α2n ) ∈ R2n which minimizes E(α) and F (α). This implies that the corresponding G1 -continuous cubic Lp spline exists for all p ∈ [1, +∞]. 2 Theorem 9. The G1 -continuous cubic Lp splines are unique for 1 < p < +∞. Proof. When 1 < p < +∞ functional (4.2) is strictly convex, hence, there exists a unique set of reals αi that minimizes this expression, which implies that the corresponding cubic Lp splines are unique. 2 4.1. G1 -continuous cubic L1 splines We have set the Ti vectors such that the minimization problem of the discretization of E(α) remains linear. To this end, we choose for the Ti vectors the ones obtained for the C 1 -continuous cubic L1 spline case. Then for the G1 -continuous cubic L1 spline case, the problem consists in finding the αi for i = 1, . . . , 2n which minimize   N n−1  2  





6t (α2i+2 Ti+1,k + α2i+1 Ti,k − 2Pi,k ) + α2i+2 Ti+1,k − α2i+1 Ti,k dt. E(α) = 1

i=0

− 12

k=1

(4.4)

382

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

The integrals in (4.4) are discretized by the midpoint rule with d equal subintervals  = E(α)

n−1  d−1  N



(−4d + 6j + 3)

(−2d + 6j + 3) (−6d + 12j + 6)

α2i+1 Ti,k + α2i+2 Ti+1,k − Pi,k

.

2 2 2 d d d i=0 j =0 k=1

The function g(t) = (1 + 6t)α2i+2 Ti+1,k + (−1 + 6t)α2i+1 Ti,k − 12tPi,k is affine, so via Lemma 4, an upper bound to discretization error of the functional (4.4) is given by n−1  N 



E(α) − E(α)   1 ej,k d2

where

i=0 k=1

ej,k = 3|α2i+2 Ti+1,k + α2i+1 Ti,k − 2Pi,k | if |α2i+2 Ti+1,k − α2i+1 Ti,k | < 3|α2i+2 Ti+1,k + α2i+1 Ti,k − 2Pi,k | and 0 otherwise. We write     α1      E(α) = A · − C .   α2n 1 where the matrix A = (ai,j ) 1id×N ×n and the vector C = (ci )1id×N ×n are defined by: 1j 2n

for 0  i  d − 1, for 0  j  n − 1 and for 0  k  N − 1 −4d + 6i + 3 Tj,k+1 , aj ×N ×d+k×d+i+1,j +1 = d2 −2d + 6i + 3 Tj +1,k+1 , aj ×N ×d+k×d+i+1,j +2 = d2 −6d + 12i + 6 Pj,k+1 . cj ×N ×d+k×d+i+1 = d2 5. Applications In the L1 case to generate the computational results presented below, the partition of the finite real interval [a, b], a = u0 < u1 < · · · < un = b is chosen uniform and the integrals are discretized by the midpoint rule with 50 equal subintervals for each integral. The tension parameter ε is set to 10−3 . The resulting linear systems are solved by the Vanderbei, Meketon and Freedman primal affine algorithm (Vanderbei et al., 1986) outlined in (Lavery, 2001). In each case, we give the number of iterations required by the primal affine method to compute the matrix of coefficients. We also give for each integral, the discretized value and the upper bound for the error of the minimal value of the discretized functional. In the following example (see Fig. 6) the data points are: (0, 0), (1, 0), (1, 1), (2, 1), (2, 0). We calculate the C 1 continuous L1 spline and the G1 -continuous L1 spline solutions. They are compared with the υ-spline (L2 ) solution (see (Hoschek and Lasser, 1993) Section 3.6.2). In this example the G1 -spline is more stretched than the C 1 -spline and the υ-spline. Our method favours abrupt changes along the axis due to the L1 -norm. In fact, contrary to (Lavery, 2000), our method separates each axis and finally the data are the couples (ui , xi ) and (ui , yi ). In next example (see Fig. 7) the C 1 -continuous cubic L1 spline and the G1 -continuous cubic L1 spline give a better solution than the C 1 -continuous cubic L1 spline of kind y = f (x) defined in (Lavery, 2000) in terms of specific objective of smoothness, overshoot and shape preservation. In (Lavery, 2006), for two-dimension data, in an example Lavery compares five types of parametric L1 and L2 splines calculated by minimizing expressions involving L1 norms, L2 norms and squares of L2 norms of second derivatives, and five types of parametric L1 and L2 splines calculated by minimizing analogous expressions involving first derivatives minus first differences among themselves. He notes that the interactive-component first-derivativebased L1 spline, which is a parametric spline (x(u), z(u)) that is obtained by minimizing

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

Number of iterations for C 1  ) Discretized values G(T Upper bound for the error

x

y

18 17.7271 0.0121

15 13.766 0.0030

383

Fig. 6. · · · G1 -continuous cubic L1 spline, —– C 1 -continuous cubic L1 spline, – – – υ-spline (L2 ). Number of iterations for G1 spline: 12.

Number of iterations for C 1  ) Discretized values G(T Upper bound for the error

x

y

18 59.4 0.0476

21 95.981 0.0428

Fig. 7. —— C 1 -continuous cubic L1 spline, – – – υ-spline (L2 ), · · · G1 -continuous cubic L1 spline, - - - C 1 -continuous cubic L1 spline of kind y = f (x) as defined in (Lavery, 2000). Data points: (0, 0), (1, 0), (2, 0), (3, 0), (3.2, 2), (4.8, 2), (5, 0), (6, 0), (7, 0), (8, 0). Number of iterations for the G1 spline: 12. Number of iterations for C 1 -continuous cubic L1 spline of kind y = f (x): 22.

Fig. 8. —– interactive-component first-derivative-based L1 spline, - - - - C 1 -continuous cubic L1 spline, for the data points of Fig. 3. With uniform parametrization (left) and chordal parametrization (right).



⎫

dx

dz ⎪ ⎪ − xi

+

− xi + − zi

⎬ 1 du

du du



du



⎪ hi ⎪ ⎪ + dx − xi − dz + zi + dz − zi ⎪ i=0 ui ⎩

du

du

⎭ du



 n



dx





(ui ) + dz (ui )

+ε

du

du

n−1 

⎧

⎨

dx u ⎪

i+1⎪

(5.1)

i=0

does an excellent job. In Figs. 8 and 9, we calculate the C 1 -continuous L1 splines which are compared with the interactive-component first-derivative-based L1 splines.

384

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

Fig. 9. —– Interactive-component first-derivative-based L1 spline, - - - - C 1 -continuous cubic L1 spline, for the data points of Figs. 1–12 of (Lavery, 2006) with chordal parametrization.

 ) according to the rotation angle applied to the data points of Fig. 6. Fig. 10. Values of G(T

Fig. 11. C 1 -continuous cubic L1 splines interpolating the same data point configuration. · · · initial data point, - - - rotation of π8 , – – – rotation of π and —– rotation of 3π . 4 8

5.1. Invariant curve under a rotation of the data points  ) for each L1 spline solution associated with the In Fig. 10, we calculate the discretized minimal value of G(T data points issued from the rotations of the initial data points of Fig. 6. As in our method the x-axis and y-axis have symmetric properties, we limit our study to rotation angles belonging to [0, π2 ]. For the data points of Fig. 6, in Fig. 11, we give the C 1 -continuous cubic L1 spline solutions for three different rotations ( π8 , π4 and 3π 8 ) and in Fig. 12, the interactive-component first-derivative-based L1 spline solutions for two π 1 different rotations ( 8 and 3π 8 ). Using a chordal parametrization, in Fig. 13, we give the C -continuous cubic L1 spline solutions and the interactive-component first-derivative-based L1 spline solutions for a rotation of π8 for the data points of Figs. 1–12 of (Lavery, 2006). The curves obtained are then superimposed and we can see (Figs. 11–13) that they are not invariant with respect to rotations. We have modified our method so as to be invariant under rotations. We create an alternate approach with a local change of coordinates on each segment [Pi , Pi+1 ].

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

385

Fig. 12. Interactive-component first-derivative-based L1 splines interpolating the same data point configuration. · · · initial data point, - - - rotation of π8 and —– rotation of 3π 8 .

Fig. 13. C 1 -continuous cubic L1 splines (left). Interactive-component first-derivative-based L1 splines (right). ——- initial data point configuration of Figs. 1–12 of (Lavery, 2006), - - - rotation of π8 . −−−−−→

−−−−−→

Pi Pi+1  Pi Pi+1  → − → On each interval [ui , ui+1 [, we choose (Pi , − u , √ ) (see i , vi ) such that the coordinates of Pi+1 are ( √ 2 2 the following graphic).

Now the functional (2.2) for the C 1 -continuous cubic L1 spline is written by 1

n−2  2

 √ −−−−−→



(1 + 6t)(βi Ti+1,1 + δi Ti+1,2 ) + (1 − 6t)Ti,1 − 6 2t Pi Pi+1  dt

u −u

i=0

i+1

i

− 12

1

−−−−→

n−2  2

 √ −

Pi Pi+1 



+

(1 + 6t)(βi Ti+1,2 − δi Ti+1,1 ) + (1 − 6t)Ti,2 − 6 2t u − u dt i=0

i+1

− 12

1

2

−−−−−→

√ −

Pn−1 Pn 

+

(1 + 6t)Tn,1 + (1 − 6t)Tn−1,1 − 6 2t

dt u −u n

n−1

− 12

1

2

−−−−−→

√ −

Pn−1 Pn 

+

(1 + 6t)Tn,2 + (1 − 6t)Tn−1,2 − 6 2t

dt u −u n

n−1

− 12

→ − → where (Ti,1 , Ti,2 ) are the coordinates of Ti in (Pi , − u i , vi ) for i = 0, . . . , n − 1.

i

386

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394 −→ −−−→ The coordinates (Tn,1 , Tn,2 ) of Tn are then expressed in the basis (Pn−1 , − u−n−1 , vn−1 ). For i = 0, . . . , n − 2, we have

βi =

−−−−−→ −−−−−−−−→ Pi Pi+1 · Pi+1 Pi+2 −−−−−→ −−−−−−−−→ Pi Pi+1 .Pi+1 Pi+2 

−−−−−−−−→ −−−−−→

Det(Pi+1 Pi+2 , Pi Pi+1 ) and δi = −−−−−→ −−−−−−−−→ . Pi Pi+1 .Pi+1 Pi+2 

We use the same algorithm as in the case of the C 1 -continuous cubic L1 spline in order to minimize this function. For the independent-component first-derivative-based L1 spline (see (Lavery, 2006)), which is a parametric spline (x(u), z(u)) that is calculated by minimizing u



i+1

n−1 n

 

dx

dx

1





du + ε − x (u ) i

i



hi du du i=0

ui

i=0

and

u



i+1

n−1 n

 

dz

dz

1





du + ε − z (u ) i

i .



hi du du i=0

ui

i=0

The functional u u



i+1

i+1

n−1 n−1  

dx

dz

1 1



− xi du + − zi

du



hi du hi du i=0

ui

i=0

ui

is written with the same notation by:   −−−−−→



 1

1 1 Pi Pi+1 

2 2 2



√ n−2  

3t − t − 4 Ti,1 + −6t + 2 h 2 i

dt

 



1 2



i=0 1

(Ti+1,1 βi + Ti+1,2 δi ) + 3t + t −

−2 4   −−−−−→ 

 1

1 1 Pi Pi+1  2 Ti,2 + −6t 2 + √ n−2  2 3t − t − 

4 2 h 2 i

  +

1 i=0 1

(Ti+1,2 βi − Ti+1,1 δi ) + 3t 2 + t − −2 4







dt





1

   −−−−−−→

   2 

1 1 1 Pn−1 Pn 

2 2 2

+ 3t − t − dt Tn−1,1 + 3t + t − Tn,1 + −6t + √ 4 4 2 hn−1 2

− 12

1

   −−−−−−→

   2 

1 1 1 Pn−1 Pn 

+

3t 2 − t − Tn−1,2 + 3t 2 + t − Tn,2 + −6t 2 + √ dt. 4 4 2 hn−1 2

− 12

The functional



⎫ ⎧

dx

dx

dz ⎪ ⎪ u





i+1 ⎪ ⎪  + − x − x − z + n−1 ⎬ ⎨

i

i i

du  1 du

du



du



⎪ hi ⎪ ⎪ + dx − xi − dz + zi + dz − zi ⎪ i=0 ui ⎩

du



⎭ du du for the interactive-component first-derivative-based L1 spline is written with the same notations by:   −−−−−→



 1

1 1 Pi Pi+1 

2 2 2



√ n−2  

3t − t − 4 Ti,1 + −6t + 2 h 2 i

dt

 



1 2



i=0 1

(Ti+1,1 βi + Ti+1,2 δi ) + 3t + t −

−2 4    −−−−−→ 

 1

1 1 √ Pi Pi+1 

2 2 3t 2 − t −

 (T + T ) −6t + 2 n−2 i,1 i,2 



4 2 hi

dt

  +



  1 2

i=0 1

Ti+1,1 (βi − δi ) + Ti+1,2 (βi + δi )

+ 3t + t − −2 4

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

387

Fig. 14. Interactive-component first-derivative-based L1 splines (right) and C 1 -continuous cubic L1 splines interpolating (left) the same data point configuration after a rotation with the invariant method.



 1

1

 2 3t 2 − t − (T − T ) i,1 i,2



4

dt

  +



  1 2

i=0 1

T + 3t + t − (β + δ ) + T (δ − β ) i+1,1 i i i+1,2 i i

−2 4   −−−−−→



 1

1 1 Pi Pi+1 

2 2 3t 2 − t −

 T + −6t + √ n−2 i,2 



4 2 2 h i

dt

  +



1 2



i=0 1

(T + t − β − T δ ) + 3t i+1,2 i i+1,1 i

−2 4 n−2 

1

   −−−−−−→

   2 

1 1 1 Pn−1 Pn 

+

3t 2 − t − Tn−1,1 + 3t 2 + t − Tn,1 + −6t 2 + √ dt 4 4 2 hn−1 2

− 12

  





1 1 2 2

 3t − t − (T (T + T ) + 3t + t − + T ) n−1,1 n−1,2 n,1 n,2



4 4

dt   −−−−−−→  +



√ 1 Pn−1 Pn  2



+ −6t + 2 1

−2 2 h 1 2

n−1

1 2

    



1 1 2 2

+ 3t − t − (Tn−1,1 − Tn−1,2 ) + 3t + t − (Tn,1 − Tn,2 )

dt 4 4 − 12 1

   −−−−−−→

   2 

1 1 1 Pn−1 Pn 

+

3t 2 − t − Tn−1,2 + 3t 2 + t − Tn,2 + −6t 2 + √ dt. 4 4 2 hn−1 2

− 12

As we can see in Figs. 14, 15a, 15b and 15c, the solutions are then invariant under rotations 5.2. Choice of the parametrization for the convexity properties Let us recall two theorems proved in (Cheng et al., 2005): Theorem 10. If four consecutive data points (xi , yi ), (xi+1 , yi+1 ), (xi+2 , yi+2 ) and (xi+3 , yi+3 ) lie on a straight line, then a cubic L1 -spline y = S(x) preserves linearity over the middle interval [xi+1 , xi+2 ].

388

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

(a)

(b)

(c) Fig. 15. (a) C 1 -continuous cubic L1 splines interpolating the same data point configuration after a rotation with the invariant method. (b) Independent-component first-derivative-based L1 splines interpolating the same data point configuration after a rotation with the invariant method. (c) Interactive-component first-derivative-based L1 splines interpolating the same data point configuration after a rotation with the invariant method.

Theorem 11. If four consecutive data points (xi , yi ), (xi+1 , yi+1 ), (xi+2 , yi+2 ) and (xi+3 , yi+3 ) take up yi+2 −yi+1 xi+2 −xi+&

 interpolant.

yi+3 −yi+2 xi+3 −xi+2 ,

yi+1 −yi xi+1 −xi



then, on the interval [xi+1 , xi+2 ], the cubic L1 -spline y = S(x) is bounded above by the linear

As we can see in Fig. 16, the L1 spline does not preserve the shape of the data points and does not satisfy linearity over the middle of three consecutive intervals. Moreover, the interpolant crosses the piecewise linear interpolant of the points even if these data points are convex according to the y-axis direction. If we change the ui to the chordal partition (see (Hoschek and Lasser, 1993) 4.4.1 page 201), then the (ui , xi ) and (ui , yi ) coordinates satisfy the hypotheses of Theorems 9 and 10 (see Fig. 19). Hence, the resulting L1 spline in the (x, y) plane does not cross the piecewise linear interpolant and is bounded above by it according to the y-axis direction (see Fig. 18).

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

Curves of Fig. 16

x

y

Number of iterations for C1

19 26.0309 0.0189

16 25.8112 0.0056

 ) Discretized values G(T Upper bound for the error

389

Fig. 16. —– C 1 -continuous cubic L1 spline, - - - - G1 -continuous cubic L1 spline. Data points: (0, 3), ( 12 , 52 ), ( 52 , 12 ), (3, 0), (4, 1), (5, 2). With a uniform parametrization.

i=n Fig. 17. The couples {(ui , xi )}i=n i=0 (left) and {(ui , yi )}i=0 (right) for the data points of Fig. 16.

Curves of Fig. 18

x

y

Number of iterations for C 1

3 4.77573 × 10−6 3.7089 × 10−9

16 16.664 0.0012

 ) Discretized values G(T Upper bound for the error

Fig. 18. —– C 1 -continuous cubic L1 spline, - - - - G1 -continuous cubic L1 spline. Data points: (0, 3), ( 12 , 52 ), ( 52 , 12 ), (3, 0), (4, 1), (5, 2) with a chordal partition a = u0 < u1 < · · · < un = b.

i=n Fig. 19. The couples {(ui , xi )}i=n i=0 (left) and {(ui , yi )}i=0 (right) for the data points of Fig. 18.

390

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

6. C 2 -continuous quintic Lp splines (k = 5) In this section, we study the C 2 case where in Definition 1 the minimal value of k is 5. We now assume that we have continuous second derivatives at each node ui , i.e., ⎧ dγ dγ ⎪ ⎪ lim (u) = lim (u) (= Ti ) ⎪ + du ⎨ u→u− du u→ui i for i = 1, . . . , n − 1, ⎪ d 2γ d 2γ ⎪ ⎪ (u) = lim (u) (= M ) lim i ⎩ du2 du2 u→u− u→u+ i i moreover we note: dγ (u0 ) = T0 , du

d 2γ (u0 ) = M0 du2

dγ (un ) = Tn , du

and

d 2γ (un ) = Mn . du2

Let T (respectively M) be the set of the n+1 vectors Ti ∈ RN (respectively Mi ∈ RN ) denoted by T = (T0 , T1 , . . . , Tn ) (respectively M = (M0 , M1 , . . . , Mn )). By a similar analysis as in Sections 3 and 4 it can easily be shown that: Theorem 12. The C 2 -continuous quintic Lp splines exist for 1  p  +∞ . Theorem 13. The C 2 -continuous quintic Lp splines are unique for 1 < p < +∞. 6.1. C 2 -continuous quintic L1 splines In the L1 case, the functional (2.2) can be written as follows: E(T , M) =

N 

Ek (T , M)

where for 1  k  N :

k=1





3

2 3



+ 15t − 6t − 60t Ti+1,k 1

2 2

     n−1 



3 3 1

+ − + 15t + 6t 2 − 60t 3 Ti,k + − − t + 3t 2 + 10t 3 hi Mi+1,k dt. Ek (T , M) =



4 2

 2  i=0 1



−2

+ − 1 + 3 t + 3t 2 − 10t 3 h M + −30t + 120t 3 P



i i,k i,k 4 2 We have to minimize each Ek (T , M). Similarly to the previous section, E(T , M) is not strictly convex, so its minima need not be unique. To this end, we propose to reduce the set of solutions and focus on vectors Ti and Mi which are the shortest (in norm) possible. Consequently, we create a C 2 -continuous cubic L1 spline minimizing the ε2 are strictly positive reals which act as new functional G(T , M) = E(T , M) + ε1 T 1 + ε2 M1 . Here ε1 and tension parameters. This functional can be written as follows G(T , M) = N k=1 Gk (T , M) where Gk (T , M) = Ek (T , M) + ε1

n 

|Ti,k | + ε2

i=0

n 

|Mi,k |.

(6.1)

i=0

For k = 1, . . . , n, we discretize the integrals in (6.1) by the midpoint rule with d equal subintervals and we obtain the following functional  ⎛  ⎞ T0,k      ⎜ . ⎟   ⎜  ⎟ T ⎜   ⎟ n,k k (T , M) = A · ⎜ G ⎟ − Ck   ⎜ M0,k ⎟   ⎝  ⎠ .     Mn,k 1 where the matrix A = (ai,j ) 1id×n+n+1 and the vector Ck = (ck,i )1id×n+n+1 are defined by: 1j 2(n+1)

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

391

Fig. 20. —– C 2 -continuous quintic L1 spline, - - - C 2 -continuous quintic L2 spline, ...... τ -spline (L2 ). Data points: (0, 0.5), (0, 0), (1, 1), (1, 0), (0, 1), (0, 0.5).

for 0  i  d − 1 and for 0  j  n − 1         2i + 1 2i + 1 2 1 2i + 1 3 −36 + 96 − 60 aj ×d+i+1,j +1 = d 2d 2d 2d     2    2i + 1 2i + 1 1 2i + 1 3 −24 + 84 aj ×d+i+1,j +2 = − 60 d 2d 2d 2d     2    2i + 1 2i + 1 1 2i + 1 3 hj 1−9 + 18 aj ×d+i+1,n+j +2 = − 10 d 2d 2d 2d         2i + 1 2i + 1 2 1 2i + 1 3 hj 3 − 12 aj ×d+i+1,n+j +3 = + 10 d 2d 2d 2d         2i + 1 2i + 1 2 1 2i + 1 3 Pi,k −60 + 180 ck,j ×d+i+1 = − 120 d 2d 2d 2d for 1  i  n + 1, and for 1  j  n + 1 if i = j then an×d+j,i = an×d+j,n+1+i = 0 else an×d+j,j = ε1 and an×d+j,n+1+j = ε2 ck,n×d+j = 0. 6.2. Examples In these examples we present C 2 -continuous quintic L1 spline and C 2 -continuous quintic L2 spline with ε = ε1 = ε2 = 10−3 . In the L1 case to generate the computational results presented below, the integrals are discretized by the midpoint rule with 50 equal subintervals for each integral. In Fig. 20 the C 2 -continuous quintic L1 spline is more stretched than the C 2 -continuous quintic L2 spline and the τ -spline (see (Hoschek and Lasser, 1993) Section 3.6.2). Remark 14. The C 2 -continuous quintic L1 splines are generally only C 2 . A C 2 -continuous quintic L1 spline is C 3 if   1  1 20Pi − 12Ti − 8Ti+1 + hi (Mi+1 − 3Mi ) = 2 20Pi−1 − 8Ti−1 − 12Ti + hi−1 (3Mi − Mi−1 ) 2 hi hi−1 for i = 1, . . . , n − 1, which is not satisfied in general. In the above example (see Fig. 20) we calculate   1 1  20Pi − 12Ti − 8Ti+1 + hi (Mi+1 − 3Mi ) − 2 20Pi−1 − 8Ti−1 − 12Ti + hi−1 (3Mi − Mi−1 ) 2 hi hi−1 for each value of i. The values are given in the following table:

392

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

Fig. 21. Curvatures of the Fig. 20 curves.

Fig. 22. —— C 2 -continuous quintic L1 spline, - - - C 2 -continuous quintic L2 spline.

Fig. 23. Curvatures of Fig. 22 curves. Number of the point

x

y

1 2 3 4

−1402.43 −63.2836 −63.2836 −1402.43

−3062.91 5230.84 −5230.84 3062.91

In Fig. 22, we use the data points of Fig. 16 with a chordal partition of the ui . As we can see the curve of the C 2 -continuous quintic L1 spline is more “stretched” than that of the C 2 -continuous quintic L2 as indicated by the curvature profiles in Fig. 23. 6.3. Invariant curve under a rotation of the data points 2 In Fig. 24, we give the L1 splines solutions for three different rotations ( π8 , π4 and 3π 8 ). The C -continuous quintic L1 splines obtained are then superimposed and we can see (Fig. 24) that these curves, as in the case of the C 1 continuous cubic L1 splines, are not invariant with respect to rotations. We use the same alternate approach of the C 1 -continuous cubic L1 splines with a local change of coordinates on each segment [Pi , Pi+1 ].

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

393

Fig. 24. C 2 -continuous quintic L1 splines interpolating the same data point configuration. · · · initial data points configuration, - - - rotation of π8 , – – – rotation of π4 and —— rotation of 3π 8 .

→ − → We recall that: on each interval [ui , ui+1 [, we choose (Pi , − u i , vi ) such that the coordinates of Pi+1 are −−−−−→ → − → The coordinates of Ti = (Ti,1 , Ti,2 ) and Mi = (Mi,1 , Mi,2 ) are written in (Pi , − u i , vi ) for i = 0, . . . , n − 1. The coordinates of Tn = (Tn,1 , Tn,2 ) and Mn = (Mn,1 , Mn,2 ) are expressed in the basis −→ −−−→ (Pn−1 , − u−n−1 , vn−1 ). For i = 0, . . . , n − 2, we have −−−−−→

Pi+1  Pi Pi+1  ( Pi√ , √ ). 2 2

βi =

−−−−−→ −−−−−−−−→ Pi Pi+1 · Pi+1 Pi+2 −−−−−→ −−−−−−−−→ Pi Pi+1 .Pi+1 Pi+2 

−−−−−−−−→ −−−−−→

and δi =

Det(Pi+1 Pi+2 , Pi Pi+1 ) −−−−−→ −−−−−−−−→ . Pi Pi+1 .Pi+1 Pi+2 

In this case, the functional (2.2) is written by:   



3

3

2 3 2 3

+ 15t − 6t − 60t (βi Ti+1,1 + δi Ti+1,2 ) + − + 15t + 6t − 60t Ti,1

1

2 2

  n−2  2





1 3

+ hi − − t + 3t 2 + 10t 3 (βi Mi+1,1 + δi Mi+1,2 )

dt



 4 2  i=0 1

−−−−−→

√ −2

+ h − 1 + 3 t + 3t 2 − 10t 3 M + 15 24t 3 − t  Pi Pi+1 



i i,1 4 2 ui+1 − ui    



3

3



+ 15t − 6t 2 − 60t 3 (βi Ti+1,2 − δi Ti+1,1 ) + − + 15t + 6t 2 − 60t 3 Ti,2

1

2 2

  n−2  2





+ hi − 1 − 3 t + 3t 2 + 10t 3 (βi Mi+1,2 − δi Mi+1,1 )

dt +



 4 2  i=0 1

−−−−−→

−2

+ h − 1 + 3 t + 3t 2 − 10t 3 M + 15√24t 3 − t  Pi Pi+1 



i i,2 4 2 ui+1 − ui   



3

3 2 3 2 3



T T + 15t − 6t + 15t + 6t − 60t + − − 60t n,1 n−1,1 1

2

2 2    



1 3 1 3 +

+ hi − − t + 3t 2 + 10t 3 Mn,1 + hi − + t + 3t 2 − 10t 3 Mn−1,1

dt 4 2 4 2



−−−−−−→

√ − 12

  P P   n−1 n

+ 15 2 4t 3 − t



un − un−1



  

3

3



+ 15t − 6t 2 − 60t 3 Tn,2 + − + 15t + 6t 2 − 60t 3 Tn−1,2 1

2 2 2



    



1 3 1 3 2 3 2 3

+ + hi − − t + 3t + 10t Mn,2 + hi − + t + 3t − 10t Mn−1,21

dt. 4 2 4 2



−−−−−−→

√ − 12

+ 15 24t 3 − t  Pn−1 Pn 



u −u n

n−1

We use the same algorithm as in the case of the C 2 -continuous quintic L1 spline in order to minimize the functional. As we can see in Fig. 25 the solutions are then invariant under rotations.

394

P. Auquiert et al. / Computer Aided Geometric Design 24 (2007) 373–394

Fig. 25. C 2 -continuous quintic L1 splines interpolating the same data point configuration after a rotation with the modified method.

References Brezis, H., 1999. Analyse fonctionnelle : Théorie et applications. Dunod, Paris. Cheng, H., Fang, S.-C., Lavery, J.E., 2002. Univariate cubic L1 splines—A geometric programming approach. Mathematical Methods of Operations Research 56, 197–229. Cheng, H., Fang, S.-C., Lavery, J.E., 2005. Shape-preserving properties of univariate cubic L1 splines. Journal of Computational and Applied Mathematics 174, 361–382. Choquet, G., 1969. Cours d’analyse, tome 2, Topologie. Masson et Cie, Paris. De Boor, C., 1978. A Practical Guide to Splines. Springer, Berlin. Farin, G., 1993. Curves and Surfaces for Computer Aided Geometric Design: a Practical Guide. Academic press, New-York. Hoschek, J., Lasser, D., 1993. Fundamentals of Computer Aided Geometric Design. A.K. Peters, Wellesley. Lavery, J.E., 2000. Univariate cubic Lp splines and shape-preserving, multiscale interpolation by univariate cubic L1 splines. Computer Aided Geometric Design 17, 319–336. Lavery, J.E., 2001. Shape-preserving, multiscale interpolation by bi- and multivariate cubic L1 splines. Computer Aided Geometric Design 18, 321–343. Lavery, J.E., 2006. Shape-preserving, first-derivative-based parametric and nonparametric cubic L1 spline curves. Computer Aided Geometric Design 23, 276–296. Saff, E.B., Tashev, F., 1999. Gibbs phenomenon for best Lp approximation by polygonal lines. East Journal on Approximations 5 (2), 235–251. Vanderbei, R.J., Meketon, M.J., Freedman, B.A., 1986. A modification of Karmarkar’s linear programming algorithm. Algorithmica 1, 395–407. Zang, Z., Martin, C.F., 1997. Convergence and Gibb’s phenomenon in cubic spline interpolation of discontinuous functions. Journal of Computational and Applied Mathematics 87, 359–371.