Constrained interpolation with rational cubics

Constrained interpolation with rational cubics

Computer Aided Geometric Design 20 (2003) 253–275 www.elsevier.com/locate/cagd Constrained interpolation with rational cubics D.S. Meek a , B.H. Ong ...

484KB Sizes 29 Downloads 182 Views

Computer Aided Geometric Design 20 (2003) 253–275 www.elsevier.com/locate/cagd

Constrained interpolation with rational cubics D.S. Meek a , B.H. Ong b,∗ , D.J. Walton a a Department of Computer Science, University of Manitoba, Winnipeg R3T 2N2, Canada b School of Mathematical Sciences, University Sains Malaysia, 11800 Penang, Malaysia

Received 9 April 2002; received in revised form 27 March 2003; accepted 9 April 2003

Abstract For a given set of ordered points lying on one side of a polyline, the problem of constructing a planar interpolating G2 curve to these planar data which lies on the same side of the polyline as the data is considered. An existing method addresses this problem using rational cubics with multiple lines as constraints and where all the data points lie strictly on one side of these lines. The method presented here extends that method to allow a polyline as a more general constraint where all the data points need not lie on one side of the infinite line through each of its edges. It also allows data points to lie on the constraint polyline. The interpolating curve is shape preserving in the sense that it has the minimal number of inflections consistent with the data.  2003 Elsevier Science B.V. All rights reserved. Keywords: Constrained interpolation; G2 rational cubic splines; Shape preserving

1. Introduction The problem of curve interpolation to given data has been studied with various requirements. One may be concerned with the smoothness of the interpolating curves, the preservation of the shape inherent in the data, the computational complexity, the precision of the curve representation, or the fulfillment of certain constraints. Here we are concerned with constrained interpolation using G2 rational cubic splines where the interpolating curve has to lie on the same side of a polyline as the interpolation points. This type of constrained interpolation could be useful in problems like designing a smooth curve on a polygonal piece of material or generating a smooth robot’s path that avoids corners or polygonal objects for a given planned path which is a sequence of vectors from point to point in the navigation of mobile robots (McKerrow, 1991). * Corresponding author.

E-mail address: [email protected] (B.H. Ong). 0167-8396/03/$ – see front matter  2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0167-8396(03)00044-X

254

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

In (Goodman et al., 1991) interpolation to data points that lie on one side of one or more lines has been considered for generating a G2 rational cubic spline which also lies on the same side of each of these lines. The non-parametric C 1 rational cubic scheme in (Ong and Unsworth, 1992) extends the linear constraints to include quadratic curves as constraints. The schemes in these two papers adjust the weights of the rational cubics so as to satisfy the conditions that a rational cubic curve does not cross a given line. The interpolating curve generated has the minimal number of inflections compatible with the data. The construction of a cubic C 1 interpolating spline curve that preserves the convexity of data points and lies on the same side of a given polygon as the data points is considered in (Zhang and Cheng, 2001). The construction process hinges on the determination of appropriate tangent vectors at the data points which can also incorporate criterion for minimising the maximum curvature of the interpolating curve. A method of constructing a guided G1 rational quadratic spline curve that falls within a closed boundary is described in (Meek et al., 2003) where the boundary is composed of straight line segments and circular arcs. This problem of constrained interpolation has not been restricted to curves and there are a number of interpolating schemes which generate range restricted surfaces interpolating data on a uniform mesh or scattered data, extending the results of the univariate cases (Brodlie et al., 1995; Chan and Ong, 2001; Herrmann et al., 1996; Mulansky and Schmidt, 1994; Schmidt and Heß, 1988; Opfer and Oberle, 1988). For a given set of ordered planar data points lying on one side of a polyline, we would like to construct an interpolating G1 curve which lies on the same side of the polyline as the data. Henceforth this polyline shall be referred to as the constraint polyline or the boundary, and a boundary segment shall mean a line segment (an edge) from the constraint polyline. The method presented extends the method in (Goodman et al., 1991) in two ways. First, it allows a polyline as the constraint while in (Goodman et al., 1991) the constraint was a set of infinite lines. Thus in (Goodman et al., 1991) all data points had to lie strictly on one side of these lines whereas here the data points need not lie on one side of the infinite line through each boundary segment. Second, it allows data points on the constraint polyline. The interpolating curve is, as in (Goodman et al., 1991), shape preserving in the sense that it has the minimal number of inflections consistent with the data. In Section 2, we describe the properties of three families of curves derived from a rational cubic by scaling its weights with positive factors. Section 3 considers finding a curve from these families which passes through a given point or is tangent to a given line. The construction of the constrained interpolating curve is described in Section 4. This involves the determination of a set of boundary avoiding λ-values for a rational cubic curve. These λ-values are positive scaling factors on the weights such that the corresponding curves obtained do not cross the constraint polyline. The cases where the data points either lie on the constraint polyline or where there are more than two collinear consecutive data points are considered in Section 5. We conclude in Section 6 by illustrating the method described with two graphical examples.

2. Families of curves derived from a rational cubic For λα , λβ ∈ (0, ∞), t ∈ (0, 1), let 3 2 2 3  λα , λβ ) = αλα (1 − t) A + 3t (1 − t) B + 3t (1 − t)C + βλβ t D , R(t, αλα (1 − t)3 + 3t (1 − t)2 + 3t 2 (1 − t) + βλβ t 3

(2.1)

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

255

(t, 1, 1) as R(t), we where α, β > 0, the control points A, B, C, D ∈ 2 . Writing the rational cubic R have R(0) = A, R(1) = D, R  (0) = 3(B − A)/α,

R  (1) = 3(D − C)/β,

(2.2) 2β[(C − B) × (D − C)] κ(1) = , 3|D − C|3  λ, 1), R 2 (t, λ) = R(t,  1, λ) where κ(0) and κ(1) are the curvatures of R(t) at 0 and 1. Let R 1 (t, λ) = R(t,  and R 3 (t, λ) = R(t, λ, λ). Three families of curves are derived from R(t) by scaling its weights α or β. They are Fi = {R i (t, λ), t ∈ [0, 1]: λ ∈ (0, ∞)}, i = 1, 2, 3. By the variation diminishing property (Goodman, 1989), when the control polyline ABCD is C-shaped the curves in all the Fi are convex or C-shaped; when ABCD is S-shaped they have an inflection or are S-shaped. 2α[(B − A) × (C − B)] , κ(0) = 3|B − A|3

2.1. Properties of R 1 (t, λ) and R 2 (t, λ) For t ∈ (0, 1), as λ → ∞ R 1 (t, λ) approaches A while R 2 (t, λ) tends to D; and when λ → 0 R 1 (t, λ) and R 2 (t, λ) respectively approach the rational quadratics 3(1 − t)2 B + 3t (1 − t)C + βt 2 D , 3(1 − t)2 + 3t (1 − t) + βt 2 α(1 − t)2 A + 3t (1 − t)B + 3t 2 C . R 2 (t, 0) = α(1 − t)2 + 3t (1 − t) + 3t 2 It is convenient for later reference to have the following simple property in a lemma. R 1 (t, 0) =

(2.3)

Lemma 1. For t ∈ (0, 1) and λ ∈ (0, ∞),   ∂R 1 (t, λ) = h1 (t, λ) A − R 1 (t, 0) , ∂λ   ∂R 2 (t, λ) = h2 (t, λ) D − R 2 (t, 0) , ∂λ where the functions h1 (t, λ) and h2 (t, λ) are positive. Proof. We only consider ∂R 1 (t, λ)/∂λ as the argument for ∂R 2 (t, λ)/∂λ is similar. Denote g1 (t, λ) = λα(1 − t)3 + 3t (1 − t)2 + 3t 2 (1 − t) + βt 3 . A simple computation gives  α(1 − t)3  ∂R 1 (t, λ) = 3t (1 − t)2 + 3t 2 (1 − t) + βt 3 A 2 ∂λ g1 (t, λ)   − 3t (1 − t)2 B + 3t 2 (1 − t)C + βt 3 D   α(1 − t)3 g1 (t, 0) A − R 1 (t, 0) , = 2 g1 (t, λ) so h1 (t, λ) =

α(1 − t)3 g1 (t, 0) > 0. g1 (t, λ)2



256

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

(a)

(b) Fig. 1. Effect of varying λ on R 1 (t, λ) where the control polyline is (a) C-shaped (b) S-shaped.

The behaviour of R 1 (t, λ) will be discussed; and that of R2 (t, λ) is analogous. By Lemma 1, for any fixed t ∈ (0, 1), as λ increases from 0 to ∞, R 1 (t, λ) moves along a line from R 1 (t, 0) towards A. This effect is shown in Fig. 1 where the curve R 1 (t, 1), t ∈ [0, 1], is drawn with a solid dark line while the limiting quadratic curve R 1 (t, 0), t ∈ [0, 1], and the curve R 1 (t, λ), t ∈ [0, 1], with a λ-value greater than 1 are drawn with a lighter shade. A set of curves where each curve does not cross any of the other curves in the set is called a nested set of curves. The behaviour of the partial derivative illustrated in Lemma 1 leads to significant and useful geometric properties. If the control polyline ABCD is C-shaped, let Ω1 be the region bounded by the AB, AD and R1 (t, 0). The segments R 1 (t, λ), λ ∈ (0, ∞) (each is an open line segment joining A to R 1 (t, 0)) for t ∈ (0, 1)

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

257

completely fill the interior of Ω1 . It is shown in Section 3.1 that for any point Q in the interior of Ω1 , there exists a unique (t, λ) such that R 1 (t, λ) = Q. It follows that the family of curves F1 is a nested set that completely fills Ω1 . When ABCD is S-shaped, the family of curves F1 is no longer nested (see Fig. 1(b)). Let U be the point of intersection of the quadratic curve R1 (t, 0) with AD at t = t0 , or U = R 1 (t0 , 0) and let T be the point of R 1 (t, 0) where AT is tangent to R 1 (t, 0) at t = t1 , or T = R 1 (t1 , 0) (see Appendix A). By Lemma 1, every curve in the family F1 intersects AD at t = t0 . Now let Ω1 be the region bounded by line segments AT, DU, BA, and the curve portions R1 (t, 0), t ∈ [0, t0 ] and t ∈ [t1 , 1]. The curves in F1 , with the exception of points A and D, lie in and fill the interior of Ω1 . Note that the open line segment AT = {R(t1 , λ): 0 < λ < ∞} forms part of the boundary of the region Ω1 , thus each curve of the family F1 is tangent to AT at t = t1 . Let Γ1 be the region bounded by line segments AU, AT and the curve R 1 (t, 0), t ∈ [t0 , t1 ]. As λ varies from 0 to ∞, the curve segments R 1 (t, 0), t ∈ (t1 , 1), and R 1 (t, 0), t ∈ (t0 , t1 ), will both sweep through the interior of the region Γ1 and this results in the nonuniqueness of the images in Γ1 . For any Q in the interior of Γ1 there exist two distinct (t, λ), one with t ∈ (t1 , 1) and the other with t ∈ (t0 , t1 ) such that R 1 (t, λ) = Q. However, for any point lying in the interior of Ω1 − Γ1 , there exists a unique (t, λ) such that R 1 (t, λ) = Q. The regions Ω2 and Γ2 are defined analogously for R 2 (t, λ). 2.2. Properties of R 3 (t, λ) For any fixed t ∈ (0, 1), as λ → ∞, R 3 (t, λ) tends to R 3 (t, ∞) =

α(1 − t)3 A + βt 3 D α(1 − t)3 + βt 3

on the line segment AD, while as λ → 0, it tends to R 3 (t, 0) = (1 − t)B + tC on the line segment BC. We note the following result for later reference. Lemma 2. For t ∈ (0, 1) and λ ∈ (0, ∞),   ∂R 3 (t, λ) = h3 (t, λ) R 3 (t, ∞) − R 3 (t, 0) ∂λ where the function h3 (t, λ) =

3t (1 − t)[α(1 − t)3 + βt 3 ] > 0. (λα(1 − t)3 + 3t (1 − t)2 + 3t 2 (1 − t) + λβt 3 )2

If the control polyline ABCD is C-shaped, the family of curves F3 is nested and completely fills the interior of Ω3 , where Ω3 is the region ABCD. In Section 3.1 it will be shown that for any Q in the interior of ABCD, there is a unique (t, λ) with R 3 (t, λ) = Q. If ABCD is S-shaped, let U be the intersection of AD with BC. The family F3 fills the interior of Ω3 which is the union of the regions ABU, CDU and some region Γ3 in the neighbourhood of U . Since U is on BC, there is a unique t0 ∈ [0, 1] such that U = R3 (t0 , 0). Let t1 be the t-value at which R(t) intersects BC. By Lemma 2 every curve R3 (t, λ) in the family will also intersect with BC at t = t1 and R 3 (t1 , ∞) = U . For the curve shown in Fig. 2(b), t0 < t1 , but depending upon the curve, it could happen that t0  t1 . Let V = R 3 (t0 , ∞) and W = R 3 (t1 , 0) (see Appendix A).

258

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

(a)

(b) Fig. 2. Effect of varying λ on R 3 (t, λ) where the control polyline is (a) C-shaped (b) S-shaped.

By Lemma 2, R 3 (t, λ), 0 < λ < ∞, is a line segment from R 3 (t, 0) to R 3 (t, ∞); R 3 (t0 , λ) is a line from U to V , and R 3 (t1 , λ) is a line from W to U as illustrated in Fig. 2(b). If t0 = t1 , the curve R 3 (t, λ) does not pass through U , then let Γ3 be the region bounded by line segments UV, UW and the curve

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

259

Fig. 3. Line segments defining E(x, y) = 0.

E(x, y) = 0 from V to W , excluding the three vertices U , V , W . Here E(x, y) = 0 is the envelope to the one parameter family of straight line segments joining R 3 (t, 0) to R 3 (t, ∞), t ∈ [t0 , t1 ] (see Fig. 3). If t0 = t1 , R 3 (t, λ) passes through U , i.e., R 3 (t0 , λ) = U , then by Lemma 2 the three points U , V and W are identical and thus Γ3 is empty. The family of curves F3 is nested in the regions bounded by the triangles ABU and UDC. In the region Γ3 it is not nested and there are two curves in the family passing through each interior point of Γ3 . Fig. 4 shows the family F3 where R 3 (t, 1) is a S-curve drawn in solid line and the envelope E(x, y) = 0 is drawn as a dotted curve. To get the envelope E(x, y) = 0, first we can find the point of intersection U = R 3 (t0 , 0) between AD and BC and parameter t0 . Then the point V = R 3 (t0 , ∞) can be obtained, and by using the relation R 3 (t0 , 0) = R 3 (t1 , ∞), t1 and W = R 3 (t1 , 0) can be found. The equation of the line through the points R 3 (t, ∞) and R 3 (t, 0) is Ψ (x, y, t) = 0,

(2.4)

where    Ψ (x, y, t) = X − (1 − t)B + tC     × α(1 − t)3 A + βt 3 D − α(1 − t)3 + βt 3 (1 − t)B + tC , X = (x, y) and × is the scalar cross product. Eq. (2.4) and ∂Ψ (x, y, t) = 0 ∂t define implicitly an equation E(x, y) = 0 for the envelope to the lines (2.4). It is still possible to determine whether a given point is in the region Γ3 . To determine if Q(x0 , y0 ) ∈ Γ3 , first determine whether the point Q lies in the triangle UVW. If it does, determine whether there exists any curve R 3 (t, λ), t ∈ [0, 1] which passes through it. This involves solving the quartic equation Ψ (x0 , y0 , t) = 0 for t in the interval with end points t0 and t1 , and then solving for λ the equation R 3 (t, λ) = Q. If there / Γ3 . exists such a curve through Q then Q ∈ Γ3 , otherwise Q ∈

260

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

Fig. 4. Family of curves F3 and the envelope E(x, y) = 0.

Fig. 5. R(t, α, β), R(t, α ∗ , β) and R(t, α ∗ , β ∗ ) where α ∗ > α and β ∗ > β.

2.3. Region for R(t, α ∗ , β ∗ ) where α ∗  α and β ∗  β (t, 1, 1) in (2.1) with weights α and β as R(t, α, β). Let us now denote the rational cubic R(t) = R When the control polyline ABCD is C-shaped, for α ∗  α and β ∗  β, R(t, α ∗ , β ∗ ) is between the curve R(t, α, β) and AD. Indeed by Lemma 1, if α is increased to α ∗ , the curve R(t, α ∗ , β) is between R(t, α, β) and AD. Similarly if β is increased to β ∗ , R(t, α ∗ , β ∗ ) is between R(t, α ∗ , β) and AD (see Fig. 5). Let the region bounded by the curve R(t, α, β) and AD be Φαβ . When ABCD is S-shaped,

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

261

Fig. 6. Region for R(t, α ∗ , β ∗ ) where α ∗  α and β ∗  β.

R(t, α ∗ , β ∗ ) lies in the region bounded by parts of the curve R(t, α, β) and by the tangents to that curve that pass through A and D respectively (see Fig. 6 and Fig. 1(b)). As above denote this region by Φαβ . When α ∗ increases from α to ∞, for fixed t, R(t, α ∗ , β) moves along a line from R(t, α, β) towards A and the curve R(t, α ∗ , β) is tangent to AT where AT is the tangent to the rational quadratic R 1 (t, 0) from A. Thus AT is a boundary of all the S-curves R(t, α ∗ , β), α ∗  α. A similar argument applies to the tangent line from D to R 2 (t, 0) in (2.3). It follows from the above discussion that the regions described above for C-shaped or S-shaped R(t, α, β) have the property Φα ∗ β ∗ ⊂ Φαβ

for α ∗  α, β ∗  β.

Moreover by Lemma 1, the region Φαβ shrinks towards AD as α, β → ∞.

3. Finding curves through given points or tangent to given line segments 3.1. The curve which interpolates a given point For any point Q which lies in the interior of Ωi , i = 1, 2, 3, we shall first consider finding a curve from the corresponding family {R i (t, λ), t ∈ [0, 1]: λ ∈ (0, ∞)} which passes through Q. When ABCD is C-shaped, it will be shown algebraically that the curve is unique. Case 1. The curve R 1 (t, λ) that interpolates a given point Q in the interior of Ω1 . To obtain this, first we find the value t = t0 where the line through AQ intersects the curve R 1 (t, λ), t ∈ [0, 1]. For Q to be a point on the curve R 1 (t, λ) and on the line, we require   3(1 − t)2 (B − A) + 3t (1 − t)(C − A) + βt 2 (D − A) × (Q − A) = 0.

262

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

As Q lies in the interior of Ω1 , this quadratic equation has at least one solution t0 ∈ (0, 1). Then for each solution t0 ∈ (0, 1), λ can be found uniquely from either of the two linear equations of R 1 (t0 , λ) = Q. Observe that the left hand side of this quadratic equation is 3(B − A) × (Q − A) when t = 0 and is 3(D − A) × (Q − A) when t = 1. If the control polyline ABCD C-shaped, these values are of opposite signs. Thus the quadratic equation has a unique solution t0 in (0, 1) and the curve which passes through Q is unique. Case 2. The curve R 2 (t, λ) that interpolates a given point Q in the interior of Ω2 . This case is similar to Case 1 above. Case 3. The curve R 3 (t, λ) that interpolates a given point Q in the interior of Ω3 . A curve from the family {R3 (t, λ), t ∈ (0, 1): λ ∈ (0, ∞)} which passes through Q is obtained by finding the value t = t0 where the line through R 3 (t0 , 0) and R 3 (t0 , ∞) passes through Q, and then getting λ from the relation R 3 (t0 , λ) = Q. For Q to be on the line through R 3 (t, 0) and R 3 (t, ∞), we require (as in (2.4))         Q − (1 − t)B + tC × α(1 − t)3 A + βt 3 D − α(1 − t)3 + βt 3 (1 − t)B + tC = 0, that is f (t) = 0 where f (t) = α(1 − t)4 (A − Q) × (B − Q) + αt (1 − t)3 (A − Q) × (C − Q) + βt 3 (1 − t)(D − Q) × (B − Q) + βt 4 (D − Q) × (C − Q). The solution(s) for this quartic equation in t can be obtained numerically (Press et al., 2002). If ABCD is C-shaped, we show below that the solution in (0, 1) of this quartic equation is unique and thus the curve R 3 (t, λ) which passes through Q is unique. The result breaks into cases depending on the position of Q with respect to the diagonals (see Fig. 7). For convenience let (A − Q) × (B − Q) = a, (A − Q) × (C − Q) = b, (D − Q) × (B − Q) = c and (D − Q) × (C − Q) = d. Notice that if Q is inside ABCD, then a and d are of opposite signs, and it suffices to consider the case where a > 0 and d < 0 as similar arguments hold for the other case. Since f (t) is in Bernstein form, with the center coefficient 0, the signs of a, b, c, d limit the number of times f (t) can be zero. If Q is in the top, left, or right regions, there is one sign change, so there must be exactly one zero of f (t). If Q is in the bottom region, b < 0, c > 0, there are three sign changes and there may be 1 or 3 zeros of f (t). Now consider the case a > 0, b < 0, c > 0, d < 0. The proof that there is actually one zero follows from intersecting two curves. Let y1 (t) and y2 (t) be     y2 (t) = −βt 3 c(1 − t) + dt . y1 (t) = α(1 − t)3 a(1 − t) + bt ,

Fig. 7. Regions cut by diagonals and position of Q relative to the diagonals.

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

263

Fig. 8. Sketch of functions y1 and y2 .

y1 (0) = αa > 0 and y1 (1) = 0, 1 is a triple zero; y2 (0) = 0, 0 is a triple zero and y2 (1) = −βd > 0. The other zero of y1 (t) is t = a/(a − b), which is in (0, 1); the other zero of y2 (t) is t = c/(c − d), which is also in (0, 1). We will now establish that c/(c − d) < a/(a − b). Let 0 < θ1 , θ2 , θ3 < π be the angles as shown in the diagram for point Q. We have θ1 + θ2 + θ3 < 2π , θ1 + θ2 > π , θ2 + θ3 > π . a = A − QB − Q sin θ1 ,

b = A − QC − Q sin(θ1 + θ2 ),

c = −D − QB − Q sin(θ2 + θ3 ),

d = −D − QC − Q sin θ3 .

The difference c bc − ad a − = a − b c − d (a − b)(c − d) and since the denominator is positive, we look at the numerator.   bc − ad = A − QB − QC − QD − Q sin θ1 sin θ3 − sin(θ1 + θ2 ) sin(θ2 + θ3 ) . Using the formulas for the product of sines and the difference of cosines, sin θ1 sin θ3 − sin(θ1 + θ2 ) sin(θ2 + θ3 )   = cos(θ1 − θ3 ) − cos(θ1 + θ3 ) − cos(θ1 − θ3 ) + cos(θ1 + 2θ2 + θ3 ) /2 = − sin(θ1 + θ2 + θ3 ) sin θ2 . The above expression is positive because π < θ1 + θ2 + θ3 < 2π . We can now see that the curves y = y1 (t) and y = y2 (t) cross exactly once for t ∈ (0, 1) as illustrated in Fig. 8 and that proves there is a unique curve R 3 (t, λ) through Q from the family. 3.2. The curve which is tangent to a given line segment Let us now write down some conditions for the multiple zeros of a Bernstein polynomial as the case of a cubic polynomial will enable us to find the curve R 3 (t, λ), if any, which is tangent to a given line segment. Let us consider

n  n j bj t (1 − t)n−j , t ∈ [0, 1], f (t) = j j =0 which is the Bernstein polynomial of degree n where bj ∈ , j = 0, . . . , n.

264

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

It is clear that f has a zero of multiplicity n at 0 (respectively at 1) if and only if b0 = b1 = · · · = bn−1 = 0,

bn = 0 (respectively b1 = b2 = · · · = bn = 0,

Lemma 3. Let

n  n j bj t (1 − t)n−j , f (t) = j j =0

b0 = 0).

t ∈ [0, 1], bj ∈ , j = 0, . . . , n.

Then f has a zero of multiplicity n in (0, 1) if and only if bj bj −1 = < 0, j = 1, . . . , n − 1. bj bj +1 In this case the zero of multiplicity n occurs at ξ ∈ (0, 1) where bj −1 ξ = < 0, j = 1, . . . , n − 1. ξ −1 bj Proof.



n  n  ξ −1 j n j kξ n t (1 − t)n−j , k(ξ − t)n = k ξ(1 − t) + (ξ − 1)t = ξ j j =0

where k ∈ . Therefore f (t) = k(ξ − t)n if and only if

j n ξ −1 , j = 0, . . . , n. bj = kξ ξ ξ < 0 if and only if ξ ∈ (0, 1). Note that ξ −1 Thus f (t) has a zero of multiplicity n at ξ ∈ (0, 1) if and only if bj ξ bj −1 < 0, j = 1, . . . , n − 1. ✷ = = bj bj +1 ξ − 1

Lemma 4. Let g(t) = a(1 − t)3 + 3bt (1 − t)2 + 3ct 2 (1 − t) + dt 3 , t ∈ [0, 1], a, b, c, d ∈ . Then g has a zero of multiplicity 2 in (0, 1) if and only if    (3.1) (ad − bc)2 = 4 ac − b2 bd − c2 and

      ac − b2 (ad − bc) < 0 equivalently bd − c2 (ad − bc) < 0 .

(3.2)

In this case, the zero of multiplicity 2 occurs at ξ ∈ (0, 1) where 2(ac − b2 ) ad − bc ξ = = < 0. ξ −1 ad − bc 2(bd − c2 ) Proof. Suppose that g has a zero of multiplicity 2 in (0, 1). First we shall consider the case where a = 0. Then   g(t) = t 3b(1 − t)2 + 3ct (1 − t) + dt 2

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

265

and by Lemma 3 g(t) has a zero of multiplicity 2 in (0, 1) if and only if 3 c 3b 2 = < 0, 3 d c 2

i.e.,

3c2 = 4bd and bc < 0, cd < 0,

and so (3.1) and (3.2) hold and this zero occurs at ξ where 2b 3c ξ = = < 0. ξ −1 c 2d Similar arguments hold for the case where d = 0. Henceforth we may assume that a = 0 and d = 0. Since g has a zero of multiplicity 2, then by Lemma 3    g(t) = γ (1 − t) + µt p(1 − t)2 + 2qt (1 − t) + rt 2 , where p/q = q/r = λ < 0, for some λ, γ , µ ∈  and the zero of multiplicity 2 occurs at ξ where ξ /(ξ − 1) = λ < 0. Equating the coefficients of g, we have a = γp,

3b = µp + 2γ q,

3c = 2µq + γ r

and

d = µr.

Eliminating γ and µ gives 3b =

2a dp 2aq + = λ2 d + , r p λ

3c =

a 2dq ar + = 2λd + 2 . r p λ

Therefore 2λ2 d = 6b −

a 4a = 3λc − λ λ

which gives a + λc = 2b. λ Similarly

(3.3)

b + λd = 2c. (3.4) λ Note that if b = 0 then by (3.3) c = 0 while if c = 0 then by (3.4) b = 0. Hence at least one of them, b or c, is not zero. Suppose that b = 0 and c = 0. Then eliminating λ from (3.3) and (3.4) gives ad 2 + 4c3 = 0 and thus (3.1) holds. Similarly if b = 0 and c = 0, then (3.1) holds. Suppose b = 0 and c = 0. From the above two equations, eliminating 1/λ gives   λ(ad − bc) = 2 ac − b2 while eliminating λ gives   ad − bc = 2 bd − c2 λ and so we obtain (3.1), i.e.,    (ad − bc)2 = 4 ac − b2 bd − c2 .

266

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

It is clear that from above that ac = b2 if and only if bd = c2 . Since g has a zero of multiplicity 2, by Lemma 3, ac = b2 and bd = c2 cannot hold simultaneously. Hence    4 ac − b2 bd − c2 > 0. As λ < 0,   ac − b2 (ad − bc) < 0. Moreover λ=

ad − bc 2(ac − b2 ) = < 0. ad − bc 2(bd − c2 )

Conversely, if (ad − bc)2 = 4(ac − b2 )(bd − c2 ) and (ac − b2 )(ad − bc) < 0, observe that by Lemma 3 g cannot have a zero of multiplicity 3 in (0, 1). Let p = 2(ac − b2 ), q = ad − bc and r = 2(bd − c2 ), then q 2 = pr. By Lemma 3, the polynomial p(1 − t)2 + 2qt (1 − t) + rt 2 has a zero of multiplicity 2 at ξ ∈ (0, 1) where 2(ac − b2 ) ad − bc ξ = = < 0. ξ −1 ad − bc 2(bd − c2 ) Let γ = a/p and µ = d/r, then by a simple calculation and using (3.1) we obtain    g(t) = γ (1 − t) + µt p(1 − t)2 + 2qt (1 − t) + rt 2 and thus g has a zero of multiplicity 2 at ξ .



Clearly, the polynomial g in Lemma 4 has a zero of multiplicity 2 at 0 (respectively at 1) if and only if a = b = 0, c = 0 (respectively c = d = 0, b = 0). Lemmas 3 and 4 give us the conditions for a non-parametric cubic Bézier polynomial curve to be tangent to the x axis. This result extends directly to the case where a parametric rational cubic Bézier curve is to be tangent to an arbitrary line. First note that a smooth rational function with a non-vanishing denominator has a multiple zero if and only if its numerator has a zero of the same multiplicity at the  1, 1), t ∈ [0, 1], given same point. Now let us consider the parametric rational cubic Bézier R(t) = R(t, in (2.1) and a given line L. Recall that the numerator of R(t) is the cubic polynomial (1 − t)3 αA + 3t (1 − t)2 B + 3t 2 (1 − t)C + t 3 βD, where A, B, C, D ∈ 2 , α, β > 0, its denominator is always positive and that the Bézier representation is affine invariant. Thus by applying Lemmas 3 and 4 and an affine transformation, the necessary and sufficient conditions for the rational cubic R(t) to be tangent to the given line L can be obtained as recorded in the following corollary. (t, 1, 1), t ∈ [0, 1], given Corollary 1. Consider the parametric rational cubic Bézier curve R(t) = R in (2.1) and a given line L. Then R(t) is tangent to L at some ξ ∈ (0, 1) if and only if one of the following two holds. c = bc = βd < 0, (i) αa b (ii) (αβad − bc)2 = 4(αac − b2 )(βbd − c2 ) and (αac − b2 )(αβad − bc) < 0 [equivalently (βbd − c2 )(αβad − bc) < 0],

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

267

where a, b, c, d are respectively the signed distances of the control points A, B, C, D from L, with points on one side of L having positive distances while points on the other side having negative distances. With condition (i) a triple zero occurs at ξ ∈ (0, 1) where αa b c ξ = = = < 0; ξ −1 b c βd with condition (ii) a zero of exact multiplicity 2 occurs at ξ ∈ (0, 1) where 2(αac − b2 ) αβad − bc ξ = = < 0. ξ −1 αβad − bc 2(βbd − c2 ) We note that to avoid getting a cusp instead of a tangent in the above corollary, the curve should be regular at the point concerned. The conditions above bear some resemblance to those derived in (Goodman et al., 1991) but there are differences. In (Goodman et al., 1991), conditions are given for a parametric rational cubic curve not to cross a given line when the two end points of the curve lie strictly on one side of the line (i.e., ad > 0). Here we are interested in the conditions for a parametric rational cubic curve to be tangent to a given line where the two end points of the curve need not lie on the same side of the given line (i.e., ad < 0 is allowed).

4. Construction of the constrained interpolating curve Suppose that we are given an ordered sequence of planar data points {I i : 0  i  n} lying on one side of a constraint polyline that forms a boundary, and the polyline joining these data points consecutively does not intersect the boundary. We would like to construct a G2 piecewise rational cubic interpolating curve that lies on the same side of the boundary as the data. Let us assume first that no three consecutive data points are collinear. The case when there are such data points is discussed in Section 5. 4.1. Initial G2 interpolating curve We shall construct an initial G2 interpolating curve through these ordered data points by using the G2 piecewise rational cubic shape preserving scheme in (Goodman, 1988). It gives a curve with the minimum  1, 1) number of inflections consistent with the data. A rational cubic segment of the form R(t) = R(t, in (2.1) is generated between each pair of consecutive data points. If I i−1 , I i , I i+1 are three data points, the curvature at I i is derived from the circle passing through those three points, and the tangent direction at I i is a linear combination of the two vectors I i − I i−1 and I i+1 − I i . The inner control points B and C are then fixed by some choice on the lengths of |B − I i | and |I i+1 − C| which reproduces circular arcs when the data come from a circle. Finally the weights α and β are defined from (2.2). The reader is referred to (Goodman, 1988) for the formulas and details. 4.2. Determining the boundary avoiding λ-values Some of the curve segments of the initial interpolating curve described in Section 4.1 could have crossed the given boundary since the existence of the boundary has not been taken into account. We shall consider the boundary now.

268

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

Each curve segment of the initial interpolating curve will be treated locally in the order of the given data points. The curve segment R(t) is the same as R 3 (t, 1), and so is a member of the family of curves F3 . We may consider replacing it by another member of that family R 3 (t, λ) whereby R 3 (t, λ) lies on one side of the boundary; a λ-value that gives this property is called a boundary avoiding λ. Since the curves R 3 (t, λ) of F3 fill Ω3 , each of the boundary segments which intersects the interior of Ω3 will be considered and the set of all such line segments is denoted by Λ. If Λ is empty, then it is clear that the set of boundary avoiding λ-values for R3 (t, λ) is (0, ∞), otherwise to determine the boundary avoiding λ-values, one or both of the following operations will be performed on each L ∈ Λ. (a) For each end point of L which lies in the interior of region Ω3 , find the rational cubic R 3 (t, λ) which passes through it (possibly more than one such curve). (b) For each L, find the rational cubic R 3 (t, λ) which touches it, if there is any (and there may be more than one). Recall that as λ decreases from ∞, the curve R 3 (t, λ) moves from the line segment AD to control polyline ABCD. Since the given boundary does not intersect AD, one could start by considering λ = ∞ and decreasing λ until the moment when some member of the family of curves F3 first hits the boundary. Since this family of curves varies continuously with λ, the first hit on the boundary would be through an end point of a boundary segment or a point tangent to a boundary segment. Let the λ of that curve be λ∗ , then [λ∗ , ∞) is a boundary avoiding set of λ. If ABCD is C-shaped, the set of all boundary avoiding λ-values is the interval [λM , ∞), where λM is the largest λ obtained from (a) and (b) for all L ∈ Λ. If ABCD is S-shaped this set may not be so easily described, however there is still some positive q where [q, ∞) is a set of boundary avoiding λ-values. The family of curves F3 is only “partially” nested and its locus contains the region Γ3 around the point U . There may be two rational cubics R 3 (t, λ) which pass through the end point of some L and there may be more than one curve which is tangent to some L. Let the largest λ obtained from all L ∈ Λ be λM , then the range [λM , ∞) is a set of boundary avoiding λ-values. Note that the default curve segment is R 3 (t, 1). After obtaining a set of boundary avoiding λ-values for R 3 (t, λ), if the value 1 belongs to this set then the default segment does not cross the constraint polyline and thus no modification is required, otherwise the smallest boundary avoiding λ larger than 1 is used. 4.3. Restoring G2 continuity When a curve segment, say segment i, denoted by R i (t), is modified by increasing both of its weights, the curvatures at its end points change accordingly and this disrupts the G2 continuity between adjacent curve segments. If both of the weights of the curve segment Ri (t) have been scaled by a factor λi , then to restore G2 continuity the weight β of the preceding segment R i−1 (t) and the weight α of the succeeding segment R i+1 (t) have to be scaled by the same factor λi (see (2.2)). When the adjacent segment is Cshaped, increasing one or both of its weights will draw it closer to the line segment AD and thus this will not cause it to cross the constraint polyline. However if the adjacent segment, say R i+1 (t) is S-shaped, then the increase in its weight α might cause it to cross some boundary segments in the region Γ1 which it initially does not. This problem will be treated in Section 4.4. In Fig. 9, the initial adjacent segment drawn as dotted lines does not cross the constraint polyline but the new segment, which is drawn as a solid line

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

269

Fig. 9. New adjacent curve segment crossing over the polyline which the initial curve segment does not.

and obtained by increasing the weight α with some scaling factor, crosses the constraint polyline in the region Γ1 . 4.4. Algorithm Based upon the above discussion, there are different approaches that could be used to generate a boundary avoiding G2 interpolating curve. We have adopted the approach described in the following algorithm. Algorithm. Given planar data points {I i : 0  i  n} where there are no three consecutive collinear points, and a constraint polyline with boundary segments {Li : 0  i  m} that do not intersect with the polyline joining the data points.  1, 1) in (2.1) (1) For i = 0, 1, . . . , n − 1, construct the rational cubic segment i of the form R(t) = R(t, interpolating I i and I i+1 as outlined in Section 4.1. (2) For i = 0, 1, . . . , n − 1, (a) determine the boundary avoiding set Si of λ-values as described in Section 4.2 for the rational cubic curve segment i by considering each of the boundary segments Lj which intersects the interior of Ω3 and performing if necessary the two operations: (i) find the rational cubic R 3 (t, λ) which passes through the end point(s) of the boundary segment concerned, (ii) find the rational cubic R 3 (t, λ) which touches the boundary segment concerned; (b) determine the scaling factor λi for the weights by considering Si . If 1 ∈ Si , then λi = 1, otherwise λi is the smallest value in Si which is larger than 1. (3) For i = 1, . . . , n − 1, (a) let λ¯ i = max{λi−1 , λi }, (b) if λ¯ i > 1, then scale the weight β of curve segment i − 1 and weight α of segment i by λ¯ i . (4) If in (3) there is some λ¯ i > 1, then repeat (2) and (3) until λ¯ i = 1 for all i.

270

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

In view of the limiting behaviour described in Section 2.3 of the region Φαβ as α, β → ∞ and the fact that there are only a finite number of constraint line segments, the algorithm will terminate in a finite number of iterations. In all the examples that were tried, at most three iterations were required. We note here that the problem mentioned in Section 4.3 could happen in step 3(b), but this problem will be solved by repeating steps (2) and (3) as described in step (4). Each repetition of (2) and (3) pulls each curve segment, if necessary, closer to the corresponding line segment AD so that eventually the curve will not intersect the boundary. From the above discussion, we have following proposition. Proposition 1. Suppose that an ordered finite set of planar data points and a constraint polyline are given. If the polyline joining these data points does not intersect the constraint polyline, then a G2 piecewise rational cubic interpolating curve which lies on the same side of the constraint polyline as the data points can be constructed through these data points and it has the minimal number of inflections consistent with the data.

5. Data points on the boundary and collinear data 5.1. One isolated data point on the boundary An isolated data point in the interior of a boundary segment can be allowed and G2 continuity can still be achieved at that data point. We must use the boundary for the tangent direction at the data point on the boundary to ensure that the curve does not cross the boundary (Fig. 10). We can also allow an isolated data point at the “V” join of two boundary segments if a tangent that stays in the interior can be chosen to keep the curve inside the boundary locally (Fig. 11).

Fig. 10. One isolated data point in the interior of a boundary segment.

Fig. 11. One isolated data point at a “V” join of two boundary segments.

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

271

5.2. More than one successive data point on a boundary segment and collinear data More than one successive data point is allowed on the boundary in special circumstances. The curve segment between two successive points on the boundary will be a straight line, with curvature 0. To preserve G2 continuity, the cubic that joins a straight line segment must have curvature 0 at one end. Since an S-shaped cubic cannot have that property, the cubic must be C-shaped. Thus the cubics joining to each end of a straight line section must both be C-shaped. If the data points from I b to I k are on the boundary (Fig. 12), the cubic joining I b−1 to I b must have its inner control points B and C on the line through I b and I k . This will force the cubic to have zero curvature at I b rather than the value that may have been calculated as in Section 4.1. In Example 2, C was chosen as the midpoint of BIb and the weight β as 1. The cubic joining at I k is treated similarly.

Fig. 12. More than one successive data point on a boundary segment.

Fig. 13. Example 1.

272

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

6. Graphical examples We shall illustrate our discussion with two examples. In both of the examples data points are marked by “•”, the constraint polyline is drawn with dotted lines, the final interpolating curve is drawn in black while the initial curve is drawn in gray. The polynomial equations involved are solved numerically by using the routines in (Press et al., 2002) for the eigenvalue method and for the balancing of matrices. Fig. 13 shows the first example which is a closed curve with 6 data points inside a polygon. Three initial curve segments have crossed the polyline. The final curve which lies inside the polygon is obtained with steps (2) and (3) of the algorithm being repeated once. The curvature plot of the curve is shown in Fig. 14. The second example shown in Fig. 15 is an open G2 curve consisting of 16 data points. Data points E and F lie in the interior of a boundary segment and the interpolating curve segment between them is

Fig. 14. Curvature plot for Example 1.

Fig. 15. Example 2.

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

273

linear. Point G lies at a vertex of the polyline. Linear segments are constructed joining the last three data points H , J and K which are collinear. The final curve is obtained with steps (2) and (3) of the algorithm being applied three times.

Acknowledgements This work was supported by the Natural Sciences and Engineering Research Council of Canada. Special thanks are due to Prof. T.N.T. Goodman for his helpful comments and suggestions which improved this paper.

Appendix A (a) For finding t0 at which the quadratic curve R 1 (t, 0) intersects AD, consider the equation of the line through A and D which is (A − X) × (A − D) = 0 where X = (x, y). Substituting X = R 1 (t0 , 0) into the above equation yields   3(1 − t0 )2 (A − B) + 3t0 (1 − t0 )(A − C) + βt02 (A − D) × (A − D) = 0 from which we obtain (A − B) × (A − D) . t0 = (C − B) × (A − D) (b) To find t1 of R 1 (t1 , 0) = T where AT is tangent to the curve R 1 (t, 0), we require at t = t1 ,  d  R 1 (t, 0) − A × R 1 (t, 0) = 0. dt Let Q(t) = 3(1 − t)2 (B − A) + 3t (1 − t)(C − A) + βt 2 (D − A) and q(t) = 3(1 − t)2 + 3t (1 − t) + βt 2 , then R 1 (t, 0) =

Q(t) + A. q(t)

Thus Q(t) q(t) × Q (t) − q  (t) × Q(t) × =0 q(t) q(t)2 which reduces to Q(t) × Q (t) = 0

274

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

or 3(1 − t)2 (B − A) × (C − A) + 2βt (1 − t)(B − A) × (D − A) + βt 2 (C − A) × (D − A) = 0. Solving this quadratic equation yields t1 in (t0 , 1). (c) We would like to find t0 such that R 3 (t0 , 0) = U where U is the intersection of AD with BC, and t1 where R3 (t1 , ∞) = U . Since R 3 (t0 , 0) = U , there exists v such that (1 − t0 )B − t0 C = (1 − v)A + vD. Then t0 (C − B) = (A − B) + v(D − A) and by taking cross products with A − D, we have t0 =

(A − B) × (A − D) . (C − B) × (A − D)

Similarly v=

(A − B) × (C − B) . (A − D) × (C − B)

As U = R 3 (t1 , ∞), (1 − v)A + vD =

α(1 − t1 )3 A + βt13 D . α(1 − t1 )3 + βt13

Comparing the coefficients of D, v= Thus



βt13 . α(1 − t1 )3 + βt13

1 − t1 t1

3 =

β(1 − v) αv

and hence t1 =

1 1 + p 1/3

where p=

β(1 − v) β(B − D) × (C − B) = . αv α(A − B) × (C − B)

References Brodlie, K.W., Butt, S., Mashwama, P., 1995. Visualization of surface data to preserve positivity and other simple constraints. Computers and Graphics 19, 585–594.

D.S. Meek et al. / Computer Aided Geometric Design 20 (2003) 253–275

275

Chan, E.S., Ong, B.H., 2001. Range restricted scattered data interpolation using convex combination of cubic Bézier triangles. J. Comput. Appl. Math. 136, 135–147. Goodman, T.N.T., 1988. Shape preserving interpolation by parametric rational cubic splines. In: International Series of Numerical Mathematics, Vol. 86. Birkhäuser Verlag, Basel, pp. 149–158. Goodman, T.N.T., 1989. Shape preserving representations. In: Lyche, T., Schumaker, L.L. (Eds.), Mathematical Methods in Computer Aided Geometric Design. Academic Press, Boston, pp. 333–351. Goodman, T.N.T., Ong, B.H., Unsworth, K., 1991. Constrained interpolation using rational cubic splines. In: Farin, G. (Ed.), Nurbs for Curve and Surface Design. SIAM, Philadelphia, pp. 59–74. Herrmann, M., Mulansky, B., Schmidt, J.W., 1996. Scattered data interpolation subject to piecewise quadratic range restrictions. J. Comput. Appl. Math. 13, 209–223. McKerrow, P.J., 1991. Introduction to Robotics. Addison-Wesley, New York. Meek, D.S., Ong, B.H., Walton, D.J., 2003. A constrained guided G1 continuous spline curve. Computer-Aided Design. Mulansky, B., Schmidt, J.W., 1994. Nonnegative interpolation by biquadratic splines on refined rectangular grids. In: Laurent, P.J., Le Mehaute, A., Schumaker, L.L. (Eds.), Wavelets, Images and Surface Fitting. A K Peters, Wellesley, pp. 379–386. Ong, B.H., Unsworth, K.K., 1992. On non-parametric constrained interpolation. In: Lyche, T., Schumaker, L.L. (Eds.), Mathematical Methods in Computer Aided Geometric Design II. Academic Press, pp. 419–430. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 2002. Numerical Recipes in C++, 2nd Edition. Cambridge University Press, Cambridge. Schmidt, J.W., Heß, W., 1988. Positivity of cubic polynomials on intervals and positive spline interpolation. Bit 28, 340–362. Opfer, G., Oberle, H.J., 1988. The derivation of cubic splines with obstacles by methods of optimization and optimal control. Numer. Math. 52, 17–31. Zhang, C., Cheng, F., 2001. Constrained shape-preserving curve interpolation with minimum curvature. Preprint.