An explicit formula for the control points of periodic uniform spline interpolants and its application

Computer Aided Geometric Design 30 (2013) 389–397

Contents lists available at SciVerse ScienceDirect

Computer Aided Geometric Design www.elsevier.com/locate/cagd

An explicit formula for the control points of periodic uniform spline interpolants and its application ✩ Chongyang Deng School of Science, Hangzhou Dianzi University, Hangzhou 310018, China

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 10 February 2011 Received in revised form 8 February 2013 Accepted 8 February 2013 Available online 13 February 2013 Keywords: Explicit formula Uniform cubic spline Curve interpolation Data polygon Periodic spline

Usually we obtain the B-spline control points of a cubic uniform spline interpolant by solving a system of linear equations. In this paper, we derive an explicit formula for the control points of periodic spline interpolants. As an application of this formula, we derive the least possible bound for the deviation of the periodic cubic spline interpolant to the data polygon. © 2013 Elsevier B.V. All rights reserved.

1. Introduction Data interpolation by cubic splines is a typical problem in Computer Aided Geometric Design (CAGD). In this paper we consider the problem of using periodic uniform cubic spline s : [0, k − 1] → Rd to interpolate a sequence of given points p 0 , p 1 , . . . , pk−1 (k  3) in Rd , where d  2, p i , p i +1 are distinct, and pk = p 0 , such that s(i ) = p i , i = 0, 1, . . . , k − 1. Let

{ p i }ki =−01 be the control points of s, then by the interpolation conditions (Prautzsch et al., 2002),

⎡

4 1 ⎢1 4 1

⎢

·

1⎢ ⎢

⎢ ⎢ ⎣

6⎢ 1

1

⎤⎡

⎤

p0 ⎥ ⎢ p1 ⎥

⎤

⎡

p0 ⎢ p1 ⎥

⎥⎢ ⎥ ⎢ ⎥ ⎥⎢ · ⎥ ⎢ · ⎥ ⎥⎢ ⎥ ⎢ ⎥ · ⎥⎢ · ⎥ = ⎢ · ⎥. ⎥⎢ ⎥ ⎢ ⎥ · ⎥⎢ · ⎥ ⎢ · ⎥ ⎦ ⎦ ⎦ ⎣ ⎣ 1 4 1 p k −2 p k −2 1 4

p k −1

(1)

p k −1

Eq. (1) can be abbreviated as N P = P . The sparse matrix N in (1) is called collocation matrix (Prautzsch et al., 2002, Section 6.8). It is well known that the linear system can be solved using direct methods (Prautzsch et al., 2002) or iterative methods (Yamaguchi, 1988; Lin et al., 2004; Lin, 2010). In this paper, we present an explicit formula for the control points of the interpolation curve. With the explicit formula, some theoretical analysis about the uniform cubic spline interpolation are easier than before. For example, while Floater (2008) derived a bound of the deviation of a uniform cubic spline interpolants from its data polygon, in this paper, we derive the least possible bound by the explicit formula, and find the conditions for achieving the least possible bound. ✩

This paper has been recommended for acceptance by H. Prautzsch. E-mail address: [email protected].

0167-8396/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cagd.2013.02.003

390

C. Deng / Computer Aided Geometric Design 30 (2013) 389–397

The rest of this paper is organized as follows. In Section 2 we give the explicit formula. As the application of the explicit formula, in Sections 3 we derive the least possible bound of the deviation from the uniform cubic spline interpolant to its data polygon. Some numerical examples are presented in Section 4. Finally we conclude the paper with discussions of future work in Section 5. 2. The explicit formula −1 Theorem 2.1. Given a sequence of points { p i }ki = in Rd , where d  2 and p i and p i +1 are distinct. Let { p i } be the control points of the 0 periodic uniform cubic spline which interpolates { p i }, then:

(i) If k is an odd number, k −1

p i = a0 p i +

2 

an ( p i −n + p i +n );

(2)

an ( p i −n + p i +n ) + a k p i − k , 2 2

(3)

n =1

(ii) If k is an even number, k −2

p i = a0 p i +

2 

n =1

where p i ±k = p i and an is defined as

√

an =

3(−1)n

(uk−n − v k−n ) + (−1)k (un − v n ) , uk + v k + 2(−1)k−1

(4)

and

u=2+

√

v =2−

3,

√

3.

To prove this theorem, we first present a lemma about the coefficients {an }. Lemma 2.2. For the coefficients {an } defined by (4),

2a0 + a1 = 3,

(5)

4an + an−1 + an+1 = 0 (1  n  k − 1),

(6)

an = ak−n

(7)

(0  n  k),

2a k + a k−2 = 0 (k is an even number),

(8)

5a k−1 + a k−3 = 0 (k is an odd number).

(9)

2

2

2

2

√

Proof. Because u = 2 +

2u − 1 = u+

1 u

3, v = 2 −

√

√

3, it is easy to verify that

√

3u ,

−4= v +

1 v

2v − 1 = − 3v ,

(10)

− 4 = 0.

(11)

By the definition of an and (10) we have

2a0 + a1 3

So (5) holds.

√

√ uk−1 − v k−1 + (−1)k (u − v ) 3 + 3[uk + v k + 2(−1)k−1 ] √ (2u − 1)uk−1 − (2v − 1) v k−1 + 2 3(−1)k−1 = = 1. √ 3[uk + v k + 2(−1)k−1 ]

=

2 3( u k − v k )

3[ u k

vk

+ 2(−1)k−1 ]

−

C. Deng / Computer Aided Geometric Design 30 (2013) 389–397

391

By the definition of an and (11), for 1  n  k − 1 we have

4an + an−1 + an+1 6

√ (uk−n − v k−n ) + (−1)k (un − v n ) = 4 3(−1)n 6[uk + v k + 2(−1)k−1 ] √ (uk−(n−1) − v k−(n−1) ) + (−1)k (un−1 − v n−1 ) + 3(−1)n−1 6[uk + v k + 2(−1)k−1 ] k −( n + 1 ) √ (u − v k−(n+1) ) + (−1)k (un+1 − v n+1 ) + 3(−1)n+1 6[uk + v k + 2(−1)k−1 ] k − n k n √ [u + (−1) u ](4 − u − u1 ) − [ v k−n + (−1)k v n ](4 − v − 1v ) = 3(−1)n 6[uk + v k + 2(−1)k−1 ] = 0.

So (6) holds. By the definition of an , for 0  n  k,

√

an =

√ = √ =

3(−1)n [uk−n − v k−n + (−1)k (un − v n )] uk + v k + 2(−1)k−1 k−n

3(−1)

[(−1)2n−k (uk−n − v k−n ) + (−1)k+2n−k (un − v n )] uk + v k + 2(−1)k−1

3(−1)k−n [uk−(k−n) − v k−(k−n) + (−1)k (uk−n − v k−n )] uk + v k + 2(−1)k−1

= ak−n .

So (7) holds. If k is an even number, by (7) we have a k−2 = a k+2 . Then by (6), 2

2a k + a k−2 = 2

4a k + a k−2 + a k+2 2

2

2

2

2

2

= 0.

If k is an odd number, by (7) we have a k−1 = a k+1 . Then by (6), 2

2

5a k−1 + a k−3 = 4a k−1 + a k−3 + a k+1 = 0. 2

2

2

2

2

These two equations mean that (8) and (9) hold.

2

Now we prove Theorem 2.1. Proof Theorem 2.1. We classify two cases, i.e., k is an odd or even number. 1 If k is an odd number, let m = k− , combining with (2), Lemma 2.2, p i −m = p i +m+1 and p i +m = p i −m−1 we have 2





m  pki−1 + 4pki + pki+1 4 = a0 p i + an ( p i −n + p i +n ) 6 6

+ =

1



n =1

a 0 ( p i −1 + p i +1 ) +

6

m 

an ( p i −1−n + p i −1+n + p i +1−n + p i +1+n )

n =1

4a0 + 2a1 6

pi +

m −1 

4an + an−1 + an+1

n =1

= pi .

6

( p i −n + p i +n ) +

5am + am−1 6

( p i −m + p i +m )

If k is an even number, let m = 2k , combining with (3), Lemma 2.2 and p i −m = p i +m we have



m −1  pki−1 + 4pki + pki+1 4 = a0 p i + an ( p i −n + p i +n ) + am p i −m 6 6

+

1 6





n =1

a 0 ( p i −1 + p i +1 ) +

m −1  n =1

an ( p i −1−n + p i −1+n + p i +1−n + p i +1+n )

392

C. Deng / Computer Aided Geometric Design 30 (2013) 389–397

+ am ( p i −1−m + p i +1−m ) =

4a0 + 2a1 6

pi +

m −1 

4an + an−1 + an+1 6

n =1

= pi .

( p i −n + p i +n ) +

4am + 2am−1 6

p i −m

These two cases mean that the periodic uniform cubic spline with control points p i (0  i  k − 1) interpolates the points p i (0  i  k − 1). 2 Remark 2.3. We get this explicit formula by solving the inverse matrix of N, which is defined by (1). Remark 2.4. For different k, the coefficients {an } are different. Another choice for the coefficients are {ank }, but they seem to be complex. So we adopt {an } as the form presented in Eq. (4). 3. Application According to Floater (2008), the global deviation of the interpolanting spline is measured by

max0i k−1 dist(s|[t i ,t i+1 ] , [ p i , p i +1 ])

μ(s) =

max0i k−1  p i +1 − p i 

Floater derives

3 4

μ

.

for uniform cubic splines. In this paper, we show that for uniform cubic spline interpolant, the √

3( 3−1)

≈ 0.27452 < 34 . Furthermore, for each k, we derive the least possible bound and find bound can be improved to be 8 the conditions to achieving the least possible bound with the aid of the explicit formula. −1 Theorem 3.1. Letting s be the periodic uniform cubic spline interpolates { p i }ki = 0 , we have

√

μ

3( 3 − 1) 8

(12)

.

Before the proof of this theorem, we first prove a lemma on the coefficients {an }. Lemma 3.2. For the coefficients {an } defined by (4),



0 < (−1)n an

0n

 |an | > |an+1 | √ √

0n

k 2 k



(13)

,

2

(14)

,

3 > a0

(k is an odd number), 2(−1)m am 3 > a0 − (k is an even number).

(15) (16)

3

Proof. It is easy to verify that (13) and (14) hold. For odd k,

√ a0 =

3( u k − v k )

uk

+

vk

+2

<

√

3.

For even k  4, letting m = 2k , and noting that um v m = 1 we have

a0 −

2(−1)m am 3

√

=

3( u k − v k )

√

−

2 2 3( u m − v m )

3 uk + v k − 2 uk + v k − 2 √ 3(um + v m )(um − v m ) − 4(um − v m ) = 3 3(um − v m )2

=

√

3 1−

4 − 6v m 3( u m

−

vm)





√

3 1−

4 − 6v 2 3( u m

−

vm)

 <

√

3.

2

C. Deng / Computer Aided Geometric Design 30 (2013) 389–397

393

In the following we give the proof of Theorem 3.1. Proof Theorem 3.1. The line segment L i (t ) := [ p i , p i +1 ] has the cubic Bézier representation



L i (t ) = (1 − t )3 p i + 3(1 − t )2 t

 + 3(1 − t )t 2

1 2

−

a0

1 2



+

a0



pi +

6



pi +

6



1

+

2

a0

2



6

1

−

a0





6



p i +1

p i +1 + t 3 p i +1 ,

0t 1

(17)

where a0 is defined by (4) By the tetrahedral algorithm to convert a B-spline representation into a Bézier representation (Prautzsch et al., 2002, Section 5.10), the curve segment s|[i ,i +1] has the cubic Bézier representation

 3

2

s|[i ,i +1] := (1 − t ) p i + 3(1 − t ) t

2

pi +

3

1 3



 + 3(1 − t )t

p i +1

2

1 3

pi +

2 3

p i +1

+ t 3 p i +1 ,

0  t  1.

(18)

Letting  p i = p i +1 − p i ,  p i = p i +1 − p i , by (17), (18) and 2a0 + a1 = 3, we have

Φ i (t ) : = L i (t ) − s|[i ,i +1]  

  1 a0 1 a0 2 1 = 3(1 − t )t (1 − t ) + pi + − p i +1 − p i − p i +1

 +t

1 2

−

a0 6

2





pi +

6

1 2

+

a0 6

2



p i +1 −

6

1 3

3

pi −

2 3



3

p i +1

   = (1 − t )t (1 − t ) 2( p i − p i ) + ( p i +1 − p i +1 ) − (a0 − 1) p i /2   + t ( p i − p i ) + 2( p i +1 − p i +1 ) + (a0 − 1) p i /2

  4p i + p i −1 + p i +1 4p i +1 + p i + p i +2 − p i + (1 + t ) − p i +1 = (1 − t )t (2 − t ) 6

− (1 − 2t ) =

3a0 − (2a0 + a1 ) 6

6

  pi

t (1 − t ) 

 (1 + t ) p i +1 − (2 − t ) p i −1 + (1 − 2t ) p i − (1 − 2t )(a0 − a1 ) p i .

6

In the following, we classify two cases, i.e., k is odd or even. 1 , by (2), we have If k is an odd number, letting m = k− 2

 p i = p i +1 − p i = a 0  p i +

m  



an ( p i −n +  p i +n ) .

n =1

Then for 0  t  1, by (19), (6), (7), (13)–(15) and p i +k = p i we have

 m   t (1 − t )     Φ i (t ) = (1 − 2t )an − (2 − t )an−1 + (1 + t )an+1  p i −n   6 n =1      + (1 − 2t )an − (2 − t )an+1 + (1 + t )an−1  p i +n    m  1    (1 − 2t )an − (2 − t )(−4an − an+1 ) + (1 + t )an+1  p i −n   24  n =1      + (1 − 2t )an − (2 − t )(−an−1 − 4an ) + (1 + t )an−1  p i +n    m        1    n n n −1 n (−1) (3 − 2t )an + an+1 (−1)  p i −n + (−1) (3 − 2t )an + an−1 (−1)  p i +n  =   8 n =1

(19)

394

C. Deng / Computer Aided Geometric Design 30 (2013) 389–397

 = =

1

m        p i  (−1)n (3 − 2t )an + an+1 + (−1)n−1 (3 − 2t )an + an−1



max

8 0i k−1 1

n =1 m     p i  a0 − a1 + (−1)m−1 am + (−1)m am+1



max

8 0i k−1

n =1

3(a0 − 1) 8

So for odd k we have

√  3( 3 − 1)  p i  <



max

8

0i k−1

μ

max

0i k−1



  p i  .

(20)

√

3( 3−1) . 8

If k is an even number, letting m = 2k , by (3) we have m −1 



 p i = p i +1 − p i = a 0  p i +



an ( p i −n +  p i +n ) + am  p i −m .

n =1

Then for 0  t  1, by (19), (6), (7) (13), (14), (16) and p i +k = p i we have

m   t (1 − t )    Φ i (t ) = (1 − 2t )an − (2 − t )an−1 + (1 + t )an+1  p i −n   6 n =1  m − 1     (1 − 2t )an − (2 − t )an+1 + (1 + t )an−1  p i +n  +  n =1 m  m −1      t (1 − t )    (3 − 2t )an + an+1  p i −n + (3 − 2t )an + an−1  p i +n  =    2 n =1 n =1 m    t (1 − t )   =  (−1)n (3 − 2t )an + an+1 (−1)n  p i −n  2 n =1  m − 1      n −1 n −1 (3 − 2t )an + an−1 (−1)  p i +n  (−1) +  n =1  m −1    m     t (1 − t ) n n −1  p i  (3 − 2t )an + an−1 (−1) (3 − 2t )an + an+1 + (−1) max  2

= = 

0i k−1

t (1 − t ) 2

n =1

max

   p i  a0 − a1 + (−1)m−2 am−1 + (−1)m−1 am + (−1)m am |1 − 2t |

0i k−1

t (1 − t ) 2



max

0i k−1



1

   p i  3(a0 − 1) − 3(−1)m am + (−1)m am |1 − 2t | 

 p i  3(a0 − 1) − 2(−1)m am max 8 0i k−1 √

<

n =1



3( 3 − 1) 8

max

0i k−1

(21)





  p i  .

(22)

√

So for even k, we have μ  3( 83−1) . By (20) and (22) we conclude that Theorem 3.1 holds. Remark 3.3. For odd k, by (20) we can see that n−1

(−1)

 p i +n (1  n 

to show that

3(a0 −1) 8

k−1 ) 2

2 3(a0 −1) if and only if (−1)n  p i −n 8 1 = 2 . In the next section we design

μ achieves the maximum

have the same directions and lengths and t

(1  n 

k−1 ), 2

such examples

is the least possible bound for odd k. √

Remark 3.4. For even k, by (22) we can see that least possible bound for even k.

3( 3−1) 8

is not the least possible bound. In the following we deduce the

C. Deng / Computer Aided Geometric Design 30 (2013) 389–397

Fig. 1. The initial data. (a) k is an odd number; (b) k is an even number and

k 2

is an odd number; (c) both k and

395

k 2

are even numbers.

For 0  t  12 , by (21) we have

μ =

t (1 − t )  2

t (1 − t )  2



3(a0 − 1) − 3(−1)m am + (−1)m am (1 − 2t ) 3(a0 − 1) − 2(−1)m am − 2(−1)m am t



= ( A − Bt )t (1 − t ) := f (t ), where A =

3(a0 −1)−2(−1)m am ,B 2  2

(23)

= (−1) am . m

Because f (t ) = 3Bt − 2( A + B )t + A, for 0  t 

1 2

we have

f (t )  f (tk ), √ ( A + B )− A 2 − A B + B 2 where tk = is a root of f  (t ) = 0 and 0 < tk < 12 . 3B 1 For 2  t  1, we have similar conclusions. So the least possible value of μ for even k is f (tk ), and μ achieves f (tk ) if 2 and only if both the directions and lengths of (−1)n  p i −n (1  n  2k ), (−1)n−1  p i +n (1  n  k− ) are equal. In the next 2 section we design such examples to show that f (tk ) is the least possible bound for even k. Combining Remarks 3.3, 3.4 and the numerical examples of the next section we have Theorem 3.5. 3(a −1)

0 (i) For odd k, the least possible value for μ is ; 8 (ii) For even k, the least possible value for μ is f (tk ).

4. Numerical examples In this section, we give some numerical examples to demonstrate the validity of the explicit formula, and show that for 3(a0 −1) is the least possible bound; for even k, the constant f (tk ) is the least possible bound. odd k, the constant 8 Let ε be a positive real number, we select the initial data as follows: For odd k, { p i } (0  i  k − 1) are defined as (see Fig. 1(a))



p i = i ε , (−1)i



(0  i  k − 1).

For even k, { p i } (0  i  k − 1) are defined as (see Fig. 1(b), (c))



p i = iε, −

1 + (−1)i +1 2



0i

k 2

−1 ,

  1 + (−1)(k−i ) k p i = (k − 1 − i )ε , 1 − i k−1 . 2

2

The initial data and interpolation curves for some k, ε are displayed in Fig. 2. In all the figures, the initial points and polygon are displayed by blue “◦” and blue lines, respectively; the interpolation curve, its control points and its control polygon are displayed by green lines, green “◦”, green lines, respectively. From Fig. 2 we can see the curves derived by the explicit formula interpolate the initial data. These numerical examples mean that the explicit formula is correct. We list the corresponding μ for ε and some odd ks, some even ks in Table 1 and Table 2, respectively. From Table 1 and 3(a0 −1) (k is an odd number) Table 2, we can see that in all these examples μ converges to the least possible bound, i.e., 8 and f (tk ) (k is an even number) when ε → 0. These results are consistent with the theoretical analysis in Section 3.

396

C. Deng / Computer Aided Geometric Design 30 (2013) 389–397

Fig. 2. The initial data and the interpolation curves. (a) k = 3, the reader is referred to the web version of this article.)

ε = 0.5; (b) k = 4, ε = 0.5; (c) k = 5, ε = 0.3. (For interpretation of the references to color,

Table 1

μ for some ε and odd k. k

3 5 7 9 11

Table 2 μ for some k

4 6 8 10 12

ε 0.5

0.2

0.1

0.05

0.02

0.01

3(a0 −1) 8

0.242536 0.264584 0.182927 0.143818 0.132895

0.248759 0.271374 0.273028 0.273147 0.273156

0.249688 0.272387 0.274048 0.274167 0.274176

0.249922 0.272642 0.274305 0.274424 0.274433

0.249988 0.272714 0.274377 0.274496 0.274505

0.249997 0.272724 0.274387 0.274506 0.274515

0.250000 0.272727 0.274390 0.274510 0.274518

0.5

0.2

0.1

0.05

0.02

0.01

f (tk )

0.188328 0.233507 0.250131 0.253949 0.255170

0.188689 0.246970 0.264544 0.269218 0.270512

0.188760 0.249275 0.267012 0.271828 0.273135

0.188778 0.249870 0.267648 0.272501 0.273811

0.188784 0.250037 0.267828 0.272691 0.274002

0.188784 0.250061 0.267853 0.272718 0.274029

0.188785 0.250069 0.267862 0.272728 0.274038

ε and even k. ε

5. Conclusion An explicit formula for uniform cubic spline interpolation is introduced in this paper. Using the explicit formula we derive the least possible bound for the deviation from the periodic uniform cubic spline interpolant to its data polygon. The conditions to achieving the least possible bound and such examples are also presented in this paper. In this paper we only consider the explicit formula for periodic uniform cubic spline interpolation. It would be interesting to know whether other spline interpolants, such as those with different parameterizations and higher degrees, have explicit formula. Acknowledgements This work is supported by NSFC (61003194). We would like to thank all anonymous reviewers for their comments and suggestions. We also wish to thank Dr. Hongwei Lin for his help in improving the language and Prof. Hartmut Prautzsch for his suggestions on the terminology.

C. Deng / Computer Aided Geometric Design 30 (2013) 389–397

397

References Floater, M.S., 2008. On the deviation of a parametric cubic spline interpolant from its data polygon. Computer Aided Geometric Design 25, 148–156. Lin, H., 2010. Local progressive-iterative approximation format for blending curves and patches. Computer Aided Geometric Design 27, 322–339. Lin, H., Wang, G., Dong, C., 2004. Constructing iterative non-uniform B-spline curves and surface to fit data points. Science in China, Series F 47, 315–331. Prautzsch, H., Boehm, W., Paluszny, M., 2002. Bézier and B-spline Techniques. Springer-Verlag. Yamaguchi, F., 1988. Curves and Surfaces in Computer Aided Geometric Design. Springer-Verlag.