Fast degree elevation and knot insertion for B-spline curves

Computer Aided Geometric Design 22 (2005) 183–197 www.elsevier.com/locate/cagd Fast degree elevation and knot insertion for B-spline curves Qi-Xing Huang a , Shi-Min Hu a,∗ , Ralph R. Martin b a Department of Computer Science and Technology, Tsinghua University, Beijing 100084, PR China b School of Computer Science, Cardiff University, Cardiff, CF24 3AA, United Kingdom Received 30 April 2004; received in revised form 15 October 2004; accepted 10 November 2004 Available online 30 November 2004 Abstract We give a new, simple algorithm for simultaneous degree elevation and knot insertion for B-spline curves. The method is based on the simple approach of computing derivatives using the control points, resampling the knot vector, and then computing the new control points from the derivatives. We compare our approach with previous algorithms and illustrate it with examples.  2004 Elsevier B.V. All rights reserved. Keywords: B-splines; Degree elevation; Knot insertion 1. Introduction Several methods have been given for degree elevation of B-spline curves (Cohen et al., 1985; Liu and Wayne, 1997; Prautzsch, 1984; Prautzsch and Piper, 1991; Pigel and Tiller, 1994), the fastest of which is Prautzsch and Piper’s algorithm (Prautzsch and Piper, 1991). Unfortunately, their algorithm suffers from being complicated and hard to implement. On the other hand, Piegl and Tiller’s (1994) algorithm is more straightforward, and easier to understand. It splits the B-spline curve into Bézier curve pieces, raises the degree of each piece, and then recombines the degree-elevated Bézier curves to produce the new B-spline curve. Liu’s algorithm (Liu and Wayne, 1997) has the benefits of being both simple to implement and fast. It takes the approach of computing the new control points using a series of knot insertions followed by a series of knot deletions. * Corresponding author. E-mail address: shimin@tsinghua.edu.cn (S.-M. Hu). 0167-8396/$ – see front matter  2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cagd.2004.11.001 184 Q.-X. Huang et al. / Computer Aided Geometric Design 22 (2005) 183–197 We give here a new fast, simple method for degree elevation of B-spline curves, which can also simultaneously insert knots. Our approach is based on the well known fact that a polynomial curve is uniquely determined by its value at a given point together with the values of all its derivatives at that point. This observation can be generalised to B-splines which are piecewise polynomials, and hence B-spline curves: each curve segment over the knot interval [ti , ti+1 ] is determined by the value of the curve and its derivatives at the knot ti . We have already used this observation to devise a knot adjustment algorithm for B-splines (Tai et al., 2003). The derivatives of a B-spline can be computed efficiently using a variety of approaches (Ferrari et al., 1994; Wang et al., 1996). It has been shown that using the inverse approach leads to a knot refinement algorithm (Sankar et al., 1994) which is significantly faster than the usual Oslo algorithm (Cohen et al., 1980). In a similar vein to (Sankar et al., 1994), this paper gives a new degree elevation algorithm based on representing a B-spline using derivatives, and converting it back to the usual piecewise polynomial form. Our new algorithm has the following advantages: • our algorithm is more efficient than existing algorithms, • our algorithm is simple to implement, with a readily understood conceptual basis, • our algorithm can handle unclamped B-spline curves, whereas previous algorithms only work in the clamped case, • our algorithm can be generalized to perform knot insertion. Section 2 outlines various B-spline formulae we need later. Section 3 describes our new algorithms. Section 4 provides examples and a comparison with other approaches. We close the paper with some conclusions and a discussion. 2. B-spline formulae Here we summarize various relevant B-spline formulae. A B-spline curve of order k is defined by a linear combination of B-spline basis functions as: P (t) = n
Pi Ni,k (t), tk−1
t
tn+1 , (1) i=0 where the Pi are control points, and Ni,k (t), i = 0, . . . , n, are B-spline basis functions defined recursively on the knot vector T = [t0 , . . . , tn+k ] as (Piegl and Tiller, 1997):  1 if ti
t < ti+1 , Ni,1 (t) = (2) 0 otherwise, Ni,k (t) = ti+k − t t − ti Ni,k−1 (t) + Ni+1,k−1 (t). ti+k−1 − ti ti+k − ti+1 (3) (By convention, if 0/0 appears in the formula, we replace it by 0.) We could explicitly require that ti+k > ti : if ti+k = ti , this leads to Ni,k (t) = 0, resulting in a B-spline curve which splits into two separate B-spline curves. However, there is no need to impose this condition, Q.-X. Huang et al. / Computer Aided Geometric Design 22 (2005) 183–197 185 and Eq. (1) is still valid with this choice of knot vector. This fact is important because the (k − p)th derivative of a B-spline curve may not satisfy the condition that ti+p > ti , but it is still convenient to treat it as a single B-spline curve. As some knots with consecutive subscripts may be equal, for the sake of convenience, we rewrite the knot vector in another form as follows: T = [t0 , . . . , tk−2 , u0 , u1 , . . . , u1 , u2 , . . . , u2 , . . . , uS−1 , . . . , uS−1 , uS , tn+2 , . . . , tn+k ],          z1 z2 (4) zS−1 where t0
· · ·
tk−2
u0 , uS
tn+2
· · ·
tn+k , and {ui }i=0,...,S is a strictly increasing sequence, with {zi }i=1,...,S−1 being a positive integer sequence giving the multiplicities of each of the knots: 1
zi
k; i = 1, 2, . . . , S − 1. The multiplicity of each ui is zi . Let P (l) (t) denote the lth derivative of P (t). Then P (l) (t) = n−l
Pil Ni+l,k−l (t), (5) i=0 where Ni+l,k−l (t) are the B-spline basis functions defined over the knot vector given by Eq. (4), and the Pil are defined recursively by:  if l = 0,  Pi l−1 l−1 k−l l ) if l > 0 and ti+k > ti+l , (P − P (6) Pi = ti+k −ti+l i+1 i  0 if l > 0 and ti+k = ti+l . l−1 Alternatively, we can compute Pi+1 from Pil−1 and Pil by a rearrangement of Eq. (6): ti+k − ti+l l Pi . k−l l−1 Pi+1 = Pil−1 + (7) j We call Pi the derivative coefficients of the B-spline P (t). When a B-spline curve has only simple knots, Wang (Wang et al., 1996) gives the following formula to compute the (k − 1)th derivatives at the knots j using the Pi as follows: P (k−1)(ui ) = Pik−1 . (8) For curves having multiple knots, we now give a similar formula: Theorem 1. j P (j ) (ui ) = Pβi , where βi = k − zi
j
k − 1, 1
i
S − 1, i l=1 zl . Proof. P (j ) (ui ) = n−j
j Pi Ni+j,k−j (ui ) = n−j
j Pi Ni+j,k−j (tβi +k−1 ) i=0 i=0 βi +k−j −1 =
i=βi j Pi Ni+j,k−j (tβi +k−1 ) (9) 186 Q.-X. Huang et al. / Computer Aided Geometric Design 22 (2005) 183–197 j = Pβi Nβi +j,k−j (tβi +k−1 ) (10) j = Pβi . (11) Eq. (10) follows from ui = tβi +k−1 = ti+j and Eq. (11) follows from βi +k−j −1 Nβi +j,k−j (tβi +k−1 ) =
Ni+j,k−j (tβi +k−1 ) = 1. ✷ i=βi Knot vectors of B-splines can be classified as clamped or unclamped (Piegl and Tiller, 1997). The knot vector of a clamped B-spline curve satisfies t0 = t1 = · · · = tk−2 = u0 and uS = tn+2 = · · · = tn+k . We may also say that P (t) is left-clamped if t0 = t1 = · · · = tk−1 . It is well known that a left-clamped B-spline curve P (t) satisfies j P (j ) (u0 ) = P0 , 0
j
k − 1. (12) Existing degree elevation algorithms (Liu and Wayne, 1997; Prautzsch and Piper, 1991; Pigel and Tiller, 1994) can only handle clamped B-spline curves, so existing approaches to raising the degree of an unclamped B-spline curve first clamp its knot vector by means of a clamping algorithm (Piegl and Tiller, 1997). As an alternative, we use a knot adjustment algorithm (Tai et al., 2003) for this purpose; it can easily be combined with our new degree elevation algorithm to obtain greater overall efficiency. 3. Degree elevation 3.1. Degree elevation of a clamped B-spline curve Since a B-spline curve is a piecewise polynomial curve, it is possible to raise its degree from k to k + m, where m is an integer greater than or equal to 1. Thus, there must exist control points Qi and a new knot vector T = [t¯0 , . . . , t¯n̄+k+m ] such that P (t) = Q(t) = n̄
i,k+m (t)Qi , N (13) i=0 i,k+m (t), i = 0, . . . , n̄, are the B-spline basis where n̄ is the number of control points of Q(t), and N functions of order k + m defined on the knot vector T. The curves P (t) and Q(t) have the same geometry and parameterization. The computation of n̄, Qi , and T is referred to as raising the degree of the curve (Pigel and Tiller, 1994). The knot vector T and n̄ can be computed as follows. Assume that T takes the form given in Eq. (4). Since degree elevation preserves continuity, Q(t) has continuity of order C k−zi at ui , and the new knot vector must take the form T = [u0 , . . . , u0 , u1 , . . . , u1 , u2 , . . . , u2 , . . . , uS−1 , . . . , uS−1 , uS , . . . , uS ],                k+m z1 +m z2 +m so that n̄ = n + S × m. We now consider how to find the Qi . zS−1 +m k+m (14) Q.-X. Huang et al. / Computer Aided Geometric Design 22 (2005) 183–197 187 3.1.1. The new degree elevation algorithm Theorem 2. The derivative coefficients of P (t) and Q(t) are related as follows j j Q0 = P0 , Qtβp +pm = 0
j
k − 1, Pβip , (15) 1
p
S − 1, k − zp
i
k − 1, k−1 Qk−1 βp +pm+j = Qβp +pm , (16) 1
p
S − 1, 1
j
m. (17) Proof. Theorem 1 gives that, P (i) (u0 ) = P0i , (i) Q (u0 ) = Qi0 , P (i) (up ) = Pβip , 0
i
k − 1, 0
i
k − 1, Q(i) (up ) = Qiβp , 1
p
S − 1, k − zp
i
k − 1. As P (t) and Q(t) have the same geometry and parametrization, so do their derivatives, which proves Eqs. (15) and (16). Consider one segment of the knot vector, t ∈ [up , up+1 ). It is well known that at most k + m of the i,k+m (t) are nonzero in this knot segment; more precisely, N i+k,m (t) is nonzero B-spline basis functions N on [up , up+1 ) when βp + (p − 1)m − k + 1
i
βp + pm − k. Consider the kth derivatives of P (t) and Q(t). As the degree of P (t) is k − 1, its kth derivative equals zero, and thus so does the kth derivative of Q(t), so that: 0=P (k) (k) (t) = Q (t) = n+S·m
i=0 βp +pm−k Qki n̄i+k,m (t) =
Qki n̄i+k,m (t). (18) i=βp +(p−1)m−k+1 Eq. (18) allows us to deduce that Qki = 0, βp + (p − 1)m − k + 1
i
βp + pm − k, and as a result we can deduce Eq. (17): t¯i+k+m − t¯i+k k k−1 Q. ✷ Qk−1 + i+1 = Qi (k + m) − k i Remark. It is obvious that Theorem 2 holds as long as P (t) and Q(t) are just left-clamped; it does not matter if the curve is right-unclamped. To do degree elevation for an unclamped curve, we can turn it into a left-clamped curve using the knot adjustment algorithm mentioned earlier. For a given l in Eqs. (6) and (7), we can see that division by a common factor of (k − l) is needed for all i, so from a software engineering point of view, it simplifies matters if we define instead j j j (k − l), Pi = Pi l=1 j j j Qi = Qi (k + m − l), l=1 which lets us rewrite Eqs. (6) and (7) in simpler form: j −1 j Pi = j −1 Pi+1 − Pi , ti+k − ti+l (19) 188 Q.-X. Huang et al. / Computer Aided Geometric Design 22 (2005) 183–197 j −1 j −1 Pi+1 = Pi j (20) 0
j
k − 1, (21) + (ti+k − ti+l )Pi . Eqs. (15)–(17) now become j j Q0 = l=1   k−l j P , k+m−l 0 j j Qβp +pm = l=1   k−l j Pβp , k+m−l k−1 Qk−1 βp +pm+i = Qβp +pm , 1
p
S − 1, k − zp
j
k − 1, 1
p
S − 1, 1
i
m. (22) (23) This leads to greater efficiency. For example, in the case where the knots of P (t) are not repeated, then j k−l l=1 ( k+m−l ) can be computed a priori, so while Eqs. (15) and (16) add a further n multiplications, Eqs. (19) and (20) save a total of n(k − 1) multiplications and mn(k − 1) divisions respectively. Based on the equations developed above, we now give a procedural method for degree elevation of a clamped B-spline curve as follows: Algorithm 1. Raise a clamped B-spline curve from degree k to degree k + m. • • • • j Use Eq. (19) to compute P0 , 0
j
k − 1 and Pβip , 1
p
S − 1, k − zp
i
k − 1. Use Eq. (14) to compute T and set ñ to n + S × m. j Use Eqs. (21)–(23) to get Q0 , 0
j
k − 1, Qiβp +pm , Qk−1 βp +pm+j . 0 Use Eq. (20) to compute new control points Qi . 3.1.2. Example We now give a numerical example showing how our algorithm works. Fig. 1 shows a clamped B-spline curve. The original curve P (t) is an order 4 B-spline curve with knot vector {0, 0, 0, 0, 0.5, 0.5, 1, 1, 1, 1} and control points {(260, 100), (100, 260), (260, 420), (420, 420), (580, 260), (420, 100)}. The curve Q(t) after degree elevation is an order 5 B-spline with knots {0, 0, 0, 0, 0, 0.5, 0.5, 0.5, 1, 1, 1, 1, 1} and with control points {(260, 100), (140, 220), (180, 340), (280, 420), (400, 420), (500, 340), (540, 220), (420, 100)}. The differential coefficients of the two curves are presented below together with the relations from Theorem 2: Q00 = P00 , 1 Q32 = P22 , 2 3 Q01 = P10 , 4 1 Q33 = P23 , 4 1 Q02 = P20 , 2 Q13 = Q03 , 1 Q03 = P30 , 4 Q43 = Q33 . Q.-X. Huang et al. / Computer Aided Geometric Design 22 (2005) 183–197 Fig. 1. Degree elevation from order 4 to order 5. Differential coefficients of P (t): P00 (260, 100) P10 (100, 260) P20 (260, 420) P01 (−320, 320) P11 (320, 320) P21 (160, 0) P02 (1280, 0) P12 (−320, −640) P22 (320, −640) P03 (−3200, −1280) P13 (320, −640) P23 (−3200, 1280) P30 (420, 420) P40 (580, 260) P50 (420, 100) P31 (320, −320) P41 (−320, −320) P32 (−1280, 0). Differential coefficients of Q(t): Q00 (260, 100) Q01 (140, 220) Q02 (180, 340) Q03 (280, 420) Q10 (−240, 240) Q11 (80, 240) Q12 (200, 160) Q13 (120, 0) Q20 (640, 0) Q21 (240, −160) Q22 (−160, −320) Q23 (160, −320) Q30 (−800, −320) Q31 (−800, −320) Q32 (160, −320) Q33 (−800, 320) Q40 (0, 0) Q41 (0, 0) Q42 (0, 0) Q43 (0, 0) Q04 (400, 420) Q05 (500, 340) Q06 (540, 220) Q07 (420, 100) Q14 (200, −160) Q15 (80, −240) Q16 (−240, −240) Q24 (−240, −160) Q25 (−640, 0) Q34 (−800, 320). 189 190 Q.-X. Huang et al. / Computer Aided Geometric Design 22 (2005) 183–197 3.2. Degree elevation of an unclamped B-spline curve The above method can readily be generalised to the case of an unclamped B-spline curve. This is done by considering the appropriate knot adjustment and degree raising results in turn, and combining them. 3.2.1. Problem statement Let P (t) = ni=0 Ni,k (t)Pi be a B-spline curve defined on the knot vector T = [t0 , . . . , tk−2 , u0 , u1 , . . . , u1 , . . . , uS−1 , . . . , uS−1 , uS , tn+2 , . . . , tn+k ].       z1 zS−1 We wish to raise the degree of this B-spline from k to k + m, and compute an appropriate new knot vector. We denote the final curve by Q(t), with knot vector of the form: T = [t¯0 , . . . , t¯k+m−2 , u0 , u1 , . . . , u1 , . . . , uS−1 , . . . , uS−1 , uS , t¯n+Sm+2 , . . . , t¯n+(S+1)m+k ].       z1 +m zS−1 As before {t¯0
t¯1
· · ·
t¯k+m−2
u0 } and {uS
t¯n+Sm+2
. . .
t¯n+(S+1)m+k } are input values which may be chosen by the user. 3.2.2. Degree elevation algorithm First, we change the knot vector of P (t) to its left-clamped form T1 = [u0 , . . . , u0 , u1 , . . . , u1 , . . . , uS−1 , . . . , uS−1 , uS , tn+1 , . . . , tn+k ].          z1 k zS−1 The new curve R(t) defined on T1 has control points Ri , i = 0, . . . , n. By considering the derivative coefficients of P (t) and R(t) arising in the knot adjustment algorithm, we have l l R0k−1 = P0k−1 Rk−2−l = Pk−2−l , Rβi p = Pβip , 0
l
k − 2. (24) 1
p
S − 1, k − zp
i
k − 1, l − (ti+k − tk−1 )Ril+1, Ril = Ri+1 (25) i = 0, . . . , k − 3 − l, l = k − 3, . . . , 0. (26) We now use the degree elevation algorithm for a clamped B-spline curve to raise the degree of R(t) as desired. Let the result be U (t) with control points Ui and knot vector T2 = [u0 , . . . , u0 , u1 , . . . , u1 , . . . , uS−1 , . . . , uS−1 , uS , t˜n+Sm+2 , . . . , t˜n+(S+1)m+k ].          z1 +m k+m j zS−1 j The relations between the Ri and the Ui are given by Eqs. (21)–(23) where we replace P by R and Q by U . Finally, we change the knot vector T2 to T and find Q(t). By considering the derivative coefficients of U (t) and Q(t), and as Qli = Uil = 0, l  k from the knot adjustment algorithm, we have j j Qβp +pm = Uβp +pm , 1
p
S − 1, k − zp
j
k − 1, k−1 k−1 k−1 Qk−1 βp +pm+i = Qβp +pm = Uβp +pm+i = Uβp +pm , j j Qm+k−1−j = Um+k−1−j , 1
j
k − 1, (27) 1
p
S − 1, 1
i
m, Qk−1 = Uik−1 , i 1
l
m, (28) (29) Q.-X. Huang et al. / Computer Aided Geometric Design 22 (2005) 183–197 j j +1 j Qi = Qi+1 − (t˜i+k+m − t˜i+j +1 )Qi , 191 i = 0, . . . , k + m − 2 − j, j = 0, . . . , k − 2. (30) These may readily be cast into algorithmic form as before. 3.3. Combining knot insertion and degree elevation 3.3.1. Problem statement In this section, we only consider the case of a clamped B-spline curve. Let P (t) = a B-spline curve defined over the knot vector n i=0 Pi Ni,k (t) be T = [u0 , . . . , u0 , u1 , . . . , u1 , u2 , . . . , u2 , . . . , uS−1 , . . . , uS−1 , uS , . . . , uS ]                z1 k z2 k zS−1 as before. We now wish to raise its degree from k to k + m, and also to insert a set of new knots si each with multiplicity yi : [s0 , . . . , s0 , s1 , . . . , s1 , . . . , sl , . . . , sl ].          y1 y0 (31) yl We denote the final curve by Q(t). We may express the final knot vector as: T = [ū0 , . . . , ū0 , ū1 , . . . , ū1 , . . . , ūl[1] , . . . , ul[1] , . . . , ūl[S]−1 , . . . , ūl[S]−1 , ūl[S] , . . . , ūl[S] ],                z̄1 k+m z̄l[1] zl[S]−1 (32) k+m where ūl[i] = ui , 0
i
S, are the knots of original curve and the knots [ū1 , . . . , ū1 , . . . , ūl[1]−1 , . . . , ūl[1]−1 , ūl[1] , . . . , ūl[1] , . . . , ūl[S]−1 , . . . , ūl[S]−1 ]             z̄1 z̄l[1]−1 z̄l[1] −z1 −m (33) z̄l[S]−1 are the ones inserted. Thus, the number of new knots is n̄ = n + Sm + li=0 yi , where number of knots being inserted. Here Eq. (33) gives another way of expressing Eq. (31). l i=0 yi is the 3.3.2. Combining knot insertion and degree elevation Existing degree elevation and knot insertion algorithms use different approaches which are hard to combine. However, our degree elevation algorithm and the knot insertion algorithm proposed by (Sankar et al., 1994) use the same idea, i.e., computing derivatives from control points, resampling the knot vector, and computing new control points from derivatives. It is thus possible to combine the two algorithms. The resulting algorithm is more efficient than performing degree elevation and knot insertion separately. The following theorem describes the relations between derivatives of the curve before and after degree elevation and knot insertion. Theorem 3. j j Q0 = P0 , j Qβ̃ 0
j
k − 1, j l[i] Qk−1 β̃ = Pβi , l[i]+h +j (34) k − zi
j
k − 1, 1
i
S − 1, (35) 1
i
S − 1, 0
h
l[i + 1] − l[i] − 1, 0
j
min(z̃l[i]+h − 1, m), (36) = Pβk−1 , i 192 Q.-X. Huang et al. / Computer Aided Geometric Design 22 (2005) 183–197 j j Qβ̃ = Qβ̃ h h  , +(k−m−1−z̃ −j ) i k + m − z̃h
j
k − 2 if h = l[i], k + m − z̃h
j
k − zi − 1 if h = l[i]. Proof. The ideas of the proof follow those of Theorem 2 and we omit them for brevity. (37) ✷ We may now give an algorithm for simultaneously raising the degree and inserting new knots, which follows by analogy with Algorithm 1. Algorithm 2. Simultaneous degree elevation and knot insertion for a clamped B-spline curve. The order is being raised from k to k + m. j • Use Eq. (19) to compute P0 , 0
j
k − 1, and Pβip , 1
p
S − 1, k − zp
i
k − 1. • Set T by Eq. (32) and set ñ = n + Sm + li=0 yi as above. j j • Use Theorem 3 to get Q0 , 0
j
k − 1, Qβ̃ , k − zi
j
k − 1, 1
i
S − 1, and l[i] Qk−1 β̃ l[i]+h , +j 1
i
S − 1, 0
h
l[i + 1] − l[i] − 1, 0
j
min(z̃l[i]+h − 1, m). • Use Eqs. (20) and (37) to compute the new control points Q0i . Our algorithm is very efficient. For a 2D B-spline curve with unique knots, we need only 3(k − 1) additions and 2(k −1) multiplications per inserted knot, while Böhm’s knot insertion algorithm (Goldman and Lyche, 1993) takes 3(k + m − 1) additions and 2(k + m − 1) multiplications. 4. Comparison and discussion We now give the results of comparing our new algorithm with existing methods. 4.1. Degree elevation Here we only consider the case of a clamped B-spline curve, as the algorithms proposed by Piegl (Pigel and Tiller, 1994) and Prautzsch (Prautzsch and Piper, 1991) cover this case. Firstly, we counted the number of operations needed by various algorithms. Prautzsch’s algorithm is faster than Piegl’s for raising the order by one, but Piegl’s algorithm is more efficient when raising the order by an arbitrary degree m. Thus we compared our algorithm with Prautzsch’s in the former case and with Piegl’s in the latter. The following tables give the number of operations required, assuming all knots are of multiplicity one; we consider knots with higher multiplicity later. Tables 1 and 2 clearly show that for arbitrary n, k, and m, our algorithm takes less arithmetic operations. Secondly, we experimentally tested and compared our algorithm to the algorithms considered above. The following tests were performed: • The growth rate as a function of m using different starting degrees (k = 2, 3, . . .). Q.-X. Huang et al. / Computer Aided Geometric Design 22 (2005) 183–197 193 Table 1 Raising the degree by one Operations Prautzsch’s algorithm Our algorithm +, − × / (n − k + 1)(10k − 8) + k (n − k + 1)(5k − 2) + 2k (n − k + 1)(2k − 2) 6(n − k + 1)(k − 1) + 2k(k − 1) (n − k + 1)(2k − 1) + k(k − 1)/2 (n − k + 1)(k − 1) + k(k − 1)/2 Table 2 Raising the degree by m Operations Pigel’s algorithm Our algorithm +, − × / ((2m + 1)(k − 1) + 2k ∗ k)(n − k + 1) ((m + 1)(k − 1) + 3k(k − 1)/2)(n − k + 1) k 2 (n − k + 1)/2 2(m + 2)(n − k + 1)(k − 1) + 2k(k − 1) (n − k + 1)(km + (k − m)) + k(k − 1)/2 (n − k + 1)(k − 1) + k(k − 1)/2 Table 3 Times (in seconds) taken to elevate the degree starting at order 3 Elevation Pigel’s algorithm Prautzsch’s algorithm Our algorithm 3 to 4 3 to 5 3 to 6 3 to 7 3 to 8 3 to 9 0.9432 0.9432 1.5162 1.8116 2.108 2.3834 0.7916 1.8412 3.1578 4.7107 6.5777 8.6342 0.5338 0.689 0.8352 0.9734 1.1446 1.3 • The growth rate as a function of the starting order k using different elevations of order (k goes to k + 1, k goes to k + 2, . . .). In an attempt to be as fair as possible to previous methods, we used the previous authors’ own code given in (Prautzsch and Piper, 1991) and (Pigel and Tiller, 1994). The hardware used was an Intel Pentium IV, 1.4 GHZ computer. We ran each program 10000 times, and counted the total time taken in seconds. Each test curve was a 2D B-spline curve with 20 randomly chosen control points and randomly distributed knots. These tests all led to similar conclusions; representative results for order 3 curves are shown in Table 3 and Fig. 2. These tests showed that: • Our algorithm is best in terms of absolute time taken. • Prautzsch and Piper’s algorithm has the worst growth characteristics; our algorithm has a slightly slower growth rate than Piegl and Tiller’s algorithm. These examples used B-spline curves with interior knots of multiplicity one. We also investigated timings in the order 4 case when the interior knots were all double or triple knots, again using a B-spline 194 Q.-X. Huang et al. / Computer Aided Geometric Design 22 (2005) 183–197 Fig. 2. Times taken to elevate the degree starting at order 3. Table 4 Times (in seconds) taken to raise the order by 3 starting at various orders Elevation Pigel’s algorithm Prautzsch’s algorithm Our algorithm 2 to 4 3 to 5 4 to 6 5 to 7 6 to 8 7 to 9 8 to 10 0.8344 1.2218 1.701 2.3283 2.8171 3.381 4.057 1.3525 1.8412 2.3516 2.9312 3.1965 3.5331 3.8652 0.4135 0.689 0.9441 1.1355 1.3819 1.564 1.7143 with 20 control points. Similar results were again obtained, again showing that our algorithm to be best both in absolute time taken, and in growth rate, when the splines have multiple knots. Next, we studied the performance of the various algorithms with respect to different starting orders, and differing amounts of degree elevation. Again, as a representative sample, we illustrate the results obtained in the order k to order k + 3 case in Table 4 and Fig. 3. The results follow a similar pattern to the previous tests, and our algorithm is the clear winner. The results of our practical experiments validated our theoretical comparisons in terms of the number of operations used. Our new algorithm is clearly more efficient for degree elevation than either of the existing algorithms used as benchmarks. 4.2. Degree elevation and knot insertion We also experimentally tested and compared our combined degree elevation and knot insertion algorithm to the use of a separate degree elevation algorithm followed by a knot insertion algorithm. For the same reasons noted above, we used Prautzsch’s degree elevation algorithm when raising the degree by one, and Piegl’s algorithm when raising the degree by more than one. We used Böhm’s (1980) knot insertion algorithm, being the fastest. Other test conditions were as described in the previous Section. Q.-X. Huang et al. / Computer Aided Geometric Design 22 (2005) 183–197 195 Fig. 3. Times taken to raise the order by 3 starting at various orders. Table 5 Times (in seconds) taken to raise the order by 1 and insert varying numbers of knots Number of knots inserted 10 30 50 70 90 110 Separate algorithms Our algorithm 0.874 1.042 1.203 1.375 1.542 1.709 0.5894 0.7106 0.8118 0.9231 1.045 1.153 Fig. 4. Times (in seconds) taken to raise the order by 1 and insert varying numbers of knots. Two tests were carried out. In the first, we always started with an order 3 curve, and raised its order to 4, while at the same time inserting some number of knots, with the number being inserted increasing from 10 to 110. Times taken are shown in Table 5, and graphed in Fig. 4. 196 Q.-X. Huang et al. / Computer Aided Geometric Design 22 (2005) 183–197 In the second test, we started with an order 4 curve, and raised its order by varying amounts to order 5, 6, . . . , 10, while at the same time always inserting 100 knots. In this case, the relative advantage of our method increased with the amount of degree elevation. Overall, our combined algorithm is faster than using separate algorithms, as expected. 4.3. Handling rounding errors Real arithmetic on a computer is not performed to perfect accuracy. Normally, the resulting errors (rounding errors) can be ignored because the noise they represent is extremely small compared to the signal. However, in certain ill-conditioned cases, the errors can be significant. Usually, normal floating point arithmetic works well for our algorithms, but the results can be unstable when the knots are distributed in an extremely non-uniform manner, e.g., if (ui+1 − ui )/(ui − ui−1 ) < 10−7 for some i. In this case, a modification to our new method may be used. We insert ui+1 and ui , k − zi+1 and k − zi times respectively, using Böhm’s approach (Böhm, 1980), so that the curve is segmented into three separate B-spline curves. Degree elevation can then be performed for each curve using our algorithm, and then Eck’s method (Eck and Hadenfeld, 1995) can be used to remove the knots ui and ui+1 , kzi and kzi+1 times respectively, to obtain the overall B-spline curve after degree elevation. Stability is achieved by dividing the B-spline curve into segments such the that knots are relatively uniformly distributed within each segment. Note that this process only takes O(k 2 ) time for each knot interval [ui , ui+1 ], which has to be processed in this way. If we do this procedure for every knot, we obtain Piegl and Tiller’s algorithm. 5. Conclusion In this paper, we have given an efficient algorithm to elevate the degree of a B-spline curve, and we have also shown how the process can be combined with a complementary algorithm working on similar principles for knot insertion, to give an efficient algorithm which can do degree elevation and knot insertion simultaneously. These methods are computationally superior to existing approaches. The new method presented in this paper can clearly also be extended to the case of degree reduction. Acknowledgement The work was supported by the Natural Science Foundation of China (Project Number 60225016, 60273012, 60321002) and the National Basic Research Project of China (Project Number 2002CB312101). References Böhm, W., 1980. Inserting new knots into B-spline curves. Computer Aided Design 12 (4), 199–201. Cohen, E., Lyche, T., Riesenfeld, R.F., 1980. Discrete B-spline subdivision techniques in computer aided geometric design and computer graphics. Computer Graphics and Image Processing 14 (2), 87–111. Q.-X. Huang et al. / Computer Aided Geometric Design 22 (2005) 183–197 197 Cohen, E., Lyche, T., Schumaker, L., 1985. Algorithms for degree raising of splines. ACM Trans. Graph. 4 (3), 171–181. Eck, M., Hadenfeld, J., 1995. Knot removal for B-spline curves. Computer Aided Geometric Design 12 (3), 259–282. Goldman, R.N., Lyche, T., 1993. Knot Insertion Deletion Algorithms for B-Spline Curves and Surfaces. Society for Industrial and Applied Mathematics. Ferrari, L.A., Sankar, P.V., Silbermann, M.J., 1994. Efficient algorithms for the implementations of general B-splines. Computer Vision, Graphics and Image Processing 56 (1), 102–105. Liu, W., Wayne, 1997. A simple, efficient degree raising algorithm for B-spline curves. Computer Aided Geometric Design 14 (7), 693–698. Prautzsch, H., 1984. Degree elevation of B-spline curves. Computer Aided Geometric Design 18 (12), 193–198. Prautzsch, H., Piper, B., 1991. A fast algorithm to raise the degree of B-spline curves. Computer Aided Geometric Design 8 (4), 253–266. Pigel, L., Tiller, W., 1994. Software-engineering approach to degree elevation of B-spline curves. Computer-Aided Design 26 (1), 17–28. Piegl, L., Tiller, W., 1997. The NURBS Book, second ed. Springer-Verlag. Sankar, P.V., Silbermann, M.J., Ferrari, L.A., 1994. Curve and surface generation and refinement based on a high speed derivative algorithm. Computer Vision, Graphics and Image Processing 56 (1), 94–101. Tai, C.-L., Hu, S.-M., Huang, Q.-X., 2003. Approximate merging of B-spline curves via knot adjustment and constrained optimization. Computer-Aided Design 35 (10), 893–899. Wang, S.Y., Ferrari, L., Silbermann, M.J., 1996. High speed computation of spline functions and applications. Internat. J. Imag. Syst. Technol. 7, 1–15.