A note on cubic polynomial interpolation

“The NURBS Book” [L. Piegl, W. Tiller, The NURBS Book, second edn, Springer, 1997] is very popular in the fields of computer aided geometric design (CAGD) and geometric modeling. In Section 9.5.2 of the book, the well-known problem of the local cubic spline approximation is discussed. The key in local cubic spline approximation is cubic polynomial interpolation. In this short paper, we present the concept of single-side/double-side cubic curves and obtain the necessary and sufficient condition of a cubic curve being a singleside/double-side curve. Based on this result, for some cases of two end tangents being nearly parallelwe present a newmethod for the problemof cubic polynomial interpolation. We also point out a flaw in Section 9.5.2 of the book and give the correction result. © 2008 Elsevier Ltd. All rights reserved. Section 9.5.2 (from Page 441 to Page 453) in [1] (hereafter referred to as the “book”) studies the following Cubic Interpolation Problem (CIP): Construct a cubic polynomial interpolating given points P0 (the start point) and P3 (the end point), tangenting the given unit vectors T0 to P0 and T3 to P3, and passing through a given point P̄. For the case T0, T3 and P3 − P0 are not coplanar and T0 and T3 are not nearly parallel, one can find a neat solution of the problem in [1]. But there is a flaw in the book, that is, the method in the book may have no solution when T0, T3 and P3−P0 are coplanar. In addition, in the case that T0 and T3 are nearly parallel we present a new method to solve the CIP. The following are the main steps of constructing the cubic polynomial in the case that T0, T3 and P3 − P0 are coplanar in the book (on page 446): 1. assign a parameter, ū, to the given point P̄ of CIP by accumulating chord lengths. In our case, ū = ‖P̄− P0‖ ‖P̄− P0‖ + ‖P3 − P̄‖ ; (1) 2. assign a tangent, Tp̄, at P̄ by (9.29) and (9.31) in the book (please find the detail on the pages of 384–386 in the book); 3. set up the following equations P̄ = sP0 + 3stP1 + 3stP2 + tP3, (2) Tp̄ × (P1 − P 2 0) = 0, (3) where s = 1− ū, t = ū, and P0 = sP0+2stP1+ tP2,P 2 1 = s P1+2stP2+ tP3. P1 and P 2 0 are obtained from the deCasteljau algorithm and lie on the line defined by P̄ and Tp̄. For the above CIP, the following result holds. ∗ Corresponding author. E-mail address: [email protected] (X. Shi). 0898-1221/$ – see front matter© 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.camwa.2008.02.032 M. Hu et al. / Computers and Mathematics with Applications 56 (2008) 1358–1363 1359 Lemma 0.1. If P3−P0, T0 and T3 are not coplanar, the necessary and sufficient condition that CIP has a solution (i.e., there exists a parameter t, 0 < t < 1, such that (2) holds) is that 0 < M < 1, (4) where M = 〈T0×T3,P̄−P0〉 〈T0×T3,P3−P0〉 . We will prove this lemma later. If we denote T = T0×T3 ‖T0×T3‖ and P3 − P0 = c0T0 + c3T3 + cT, P̄− P0 = c̄0T0 + c̄3, T3 + c̄T then M = 〈T0 × T3, P̄− P0〉 〈T0 × T3,P3 − P0〉 = c̄ c . (5) Eq. (5) shows that the geometric meaning of (4) is that the projections of P̄− P0 and P3 − P0 to T have the same direction but the projection of P̄− P0 is shorter. The key of CIP is to find P1 and P2 such that (2) holds. But we found that in the case T0 ‖ T3, (2) does not always have a solution if the parameter ū is priorly assigned. The following is a counter-example that (2) has no solution. Assuming P1 = P0 + αT0 and P2 = P3 + β T3 and substituting them into (2) yield 3αstT0 + 3stβT3 = P̄− (s3 + 3st)P0 − (t3 + 3st)P3. (6) (6) is the same as (9.107) in [1]. Let P0 = (0, 0),P3 = (3, 0)T, T0 = T3 = (0, 1)T and P̄ = (1, ȳ)T, then (6) is equivalent to { 1− 3(t3 + 3st2) = 0, 3αs2t + 3st2β = ȳ, (7) where s = 1− t and t = ū. Since the first equation of (7) does not rely on α and β, it should be an identity if (6) has a solution. However, this is impossible if the parameter ū of P̄ is priorly assigned. In the current example, from (1), ū = √ 1+ȳ2 √ 1+ȳ2+ √ 4+ȳ2 . Let ȳ = √ 3. Then t = ū = 2 2+ √ 7 , s = √ 7 2+ √ 7 . Denoting f (ȳ) = 1−3(t3+3st2), we have f ( √ 3) = 1− 3(2 3 +3·22 √ 7) (2+ √ 7)3 = 26−17 √ 7 (2+ √ 7)3 6= 0. Therefore, (6) has no solution. From (6), it holds that P̄− P0 = 3αstT0 + 3βstT3 + (t3 + 3st)(P3 − P0). (8) The above equation shows that P̄ has to be on the same plane, sayΠ , determined by P0,P3, T0 and T3. In addition, we rewrite (3) as follows ((9.108) in [1]) s(s− 2t)(Tp̄ × T0)α+ t(2s− t)(Tp̄ × T3)β = 2st(Tp̄ × (P0 − P3)). (9) Proof of Lemma 0.1. According to (8), it holds that t3 + 3st2 = M. (10) Therefore, CIP has a solution is equivalent to (10) has a root in the interval (0, 1). Denoting f (t) = t3+3st2−M = 3t2−2t3−M, we have f (t) = 6t(1− t) > 0 for 0 < t < 1. Therefore, the necessary and sufficient condition of f (t) = 0 having a solution in the interval (0, 1) is that f (0) = −M < 0 and f (1) = 1 − M > 0, i.e, (4) holds. If the condition (4) holds, we can obtain a parameter ū, the only root of f (t) = 0 in interval (0, 1), for P̄. Lemma 0.1 is proved. ♣ In the following, we discuss case by case the solution of CIP in the case T0, T3 and P3 − P0 being coplanar. Case 1: T0 and T3 are not parallel. This case is discussed in the book (last two lines on page 446) Case 2: T0 ‖ T3. In this case, there exist four constants a, b, ā, b̄ such that P̄−P0 = āT0+ b̄T⊥0 and P3−P0 = aT0+ bT ⊥ 0 , where T⊥0 is a unit vector which is perpendicular to both T0 and the normal vector of Π . On the other hand, T3 = εT0 (ε = ±1) because T3 ‖ T0 and both of them are the unit vectors. Thus, (8) can be converted to { b(t3 + 3st2) = b̄, 3αs2t + 3εst2β = ā− a(t3 + 3st2). (11) We can reasonably assume that 0 < t = ū < 1. By this assumption, the first equation of (11) shows that b and b̄ have the same sign. Thus, P̄− P0 ‖ T0 if and only if P3 − P0 ‖ T0. 1360 M. Hu et al. / Computers and Mathematics with Applications 56 (2008) 1358–1363 Case 2.1: P3−P0 is not parallel to T0. In this case, we can not priorly assign a parameter for P̄. We have to obtain a parameter t, 0 < t < 1, for P̄ from the first equation of (11). Since P3 −P0 is not parallel to T0, we know that both b and b̄ are not zeros. Similar to Lemma 0.1, in this case the CIP has a solution if and only if