13.472J/1.128J/2.158J/16.940J COMPUTATIONAL GEOMETRY Lectures 4 and 5 N. M. Patrikalakis Massachusetts Institute of Technology Cambridge, MA 02139-4307, USA Copyright c 2003 Massachusetts Institute of Technology Contents 4 Introduction to Spline Curves 2 4.1 Introduction to parametric spline curves . . . . . . . . . . . . . . . . . . . . . . 2 4.2 Elastic deformation of a beam in bending . . . . . . . . . . . . . . . . . . . . . 2 4.3 Parametric polynomial curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 4.3.1 Ferguson representation . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 4.3.2 Hermite-Coons curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 4.3.3 Matrix forms and change of basis . . . . . . . . . . . . . . . . . . . . . . 7 4.3.4 B ezier curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 4.3.5 B ezier surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.4 Composite curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.4.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.4.2 Lagrange basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.4.3 Continuity conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.4.4 B ezier composite curves (splines) . . . . . . . . . . . . . . . . . . . . . . 20 4.4.5 Uniform Cubic B-splines . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Bibliography 28 Reading in the Textbook Chapter 1, pp.6 - pp.33 1 Lecture 4 Introduction to Spline Curves 4.1 Introduction to parametric spline curves Parametric formulation x = x(u); y = y(u); z = z(u) or R = R(u) (vector notation) Usually applications need a nite range for u (e.g. 0 u 1). For free-form shape creation, representation, and manipulation the parametric representa- tion is preferrable, see Table 1.2 in textbook. Furthermore, we will use polynomials for the following reasons: Cubic polynomials are good approximations of physical splines. (Historical note: Shape of a long exible beam constrained to pass through a set of points ! Draftsman’s Splines). Parametric polynomial cubic spline curves are the \smoothest" curves passing through a set of points; (i.e. they minimize the bending strain energy of the beam / RL0 2ds). 4.2 Elastic deformation of a beam in bending Within the Euler beam theory: R y dA dAN.A. y Fiber Length ds d Center of curvature a0a1a0a1a0a1a0 a0a1a0a1a0a1a0 a2a1a2a1a2a1a2 a2a1a2a1a2a1a2 Figure 4.1: Di erential segment of an Euler beam. 2 Plane sections of the beam normal to the neutral axis (N.A.) remain plane and normal to the neutral axis after deformation. Elongation strain of a ber at distance y above the N.A.: = ( y + R)d Rd Rd = yR The radius of curvature can be related to the di erential angle and arc length: Rd = ds ) = 1R = d ds The stress experienced by a ber at distance y from the N.A. can be found from Hooke’s law relating stress and strain. = E = E yR The bending moment about the N.A. can be calculated by integrating the stress due to bending times the moment arm over the surface of the cross-section: M = Z A ydA = ER Z A y2dA = EI 1R = EI d ds = EI : For small de ections, the displacement Y = Y (x) of the N.A., the following approxima- tions can be made: dYdx ; ds dx; d ds d 2Y dx2 Therefore, MEI d 2Y dx2 (4.1) Consider a beam with simple supports and no distributed loads, as shown in Figure 4.2. Figure 4.2: Simply supported beam. Looking at a section of the beam between two supports (Figure 4.3), the bending moment can be expressed as a linear equation: M = A0 + A1x (4.2) where A0 and A1 are the bending moment and shearing force at x = 0 in Figure 4.3. 3 x A1 M(x)A1A0 Figure 4.3: Section of simply supported beam between pins. The shear force and bending moment at the ends of the section are illustrated. Upon substituting Equation 4.2 into Equation 4.1 and integrating, we get: = dYdx = 1EI [A0x + A1 x 2 2 ] + A2 Y (x) = C0 + C1x + C2x2 + C3x3 (cubic) Therefore, in order to replicate the shape of physical splines, the CAD community developed shape representation methods based on cubic polynomials. Generically, polynomials have the following additional advantages: easy to store as sequences of coe cients; e cient to compute and trace e ciently; easy to di erentiate, integrate, and adapt to matrix and vector algebra; and easy to piece together to construct composite curves with a certain order of continuity, a feature important in increasing complexity of a curve or surface. 4.3 Parametric polynomial curves 4.3.1 Ferguson representation In 1963, Ferguson at Boeing developed a polynomial representation of space curves: R(u) = a0 + a1u + a2u2 + a3u3 (4.3) where 0 u 1 by convention. Note there are 12 coe cients, ai, de ning the curve R(u). This representation is also known as the power basis or monomial form. The coe cients, ai, are di cult to interpret geometrically, so we can express ai in terms of R(0), R(1), _R(0), and _R(1) (where _R denotes derivative with respect to u): R(0) = a0 R(1) = a0 + a1 + a2 + a3 _R(0) = a1 _R(1) = a1 + 2a2 + 3a3 4 Solving the above equations yields expression for the coe cients, ai’s, in terms of the geometric end conditions of the curve. a0 = R(0) a1 = _R(0) a2 = 3[R(1) R(0)] 2 _R(0) _R(1) a3 = 2[R(0) R(1)] + _R(0) + _R(1) 4.3.2 Hermite-Coons curves By substituting the coe cients into the Ferguson representation, we can rewrite Equation 4.3 as: R(u) = R(0) 0(u) + R(1) 1(u) + _R(0) 2(u) + _R(1) 3(u) (4.4) where 0(u) = 1 3u2 + 2u3 1(u) = 3u2 2u3 2(u) = u 2u2 + u3 3(u) = u3 u2 The new basis functions, i(u), are known as Hermite polynomials or blending functions, see Figure 4.4. They were rst used for 3D curve representation in a computer environment in the 60’s by the late Steven Coons, an MIT professor, participant in the famous ARPA project MAC. 0.0 1.0 0.0 1.0 0 1 2 3 Figure 4.4: Plot of Hermite basis functions. Note that these basis functions, i(u), satisfy the following boundary conditions, see Fig- ure 4.4: 0(0) = 1; 0(1) = 00(0) = 00(1) = 0; 1(1) = 1; 1(0) = 01(0) = 01(1) = 0; 02(0) = 1; 2(0) = 02(1) = 02(1) = 0; 03(1) = 1; 3(0) = 03(0) = 03(1) = 0: These boundary conditions also allow the computation of the cubic Hermite polynomials, i(u), by setting up and solving a system of 16 linear equations in the coe cients of these polynomials. 5 Geometric Interpretation of Hermite-Coons curves. If we make the following substitutions: _R(0) = 0t(0) _R(1) = 1t(1) where t(u) is the unit tangent of the curve, the following observations can be made from Figure 4.5, relating the coe cients 0 and 1 to the shape of the curve: 0 "; 1 = const ) bias towards t(0) 1 "; 0 = const ) bias towards t(1) 0; 1 ") increase fullness 0; 1 both LARGE ) cusp forms ( _R(u ) = 0) or self-intersection occurs R(1)t(0) t(1)R(0) simultaneously 0, 1 increasing 1 constant 0 increasing t(1)t(0) R(1)R(0) Figure 4.5: Fullness of a Hermite-Coons curve. Due to the importance of choosing 0 and 1 appropriately, the Hermite-Coons represen- tation can be di cult for designers to use e ciently, but is much easier to understand than the Ferguson or monomial form. 6 4.3.3 Matrix forms and change of basis Monomial or Power or Ferguson form R(u) = h u3 u2 u 1 i 2 66 64 a3 a2 a1 a0 3 77 75 = U FM Hermite-Coons R(u) = h u3 u2 u 1 i 2 66 64 2 2 1 1 3 3 2 1 0 0 1 0 1 0 0 0 3 77 75 2 66 64 R(0) R(1) _R(0) _R(1) 3 77 75 = U MH FH Conversion between the two representations is simply a matter of matrix manipulation: U FM = U MH FH Hence, FM = MH FH FH = M 1H FM 4.3.4 B ezier curves B ezier developed a reformulation of Ferguson curves in terms of Bernstein polynomials for the UNISURF System at Renault in France in 1970. This formulation is expressed mathematically as follows: R(u) = nX i=0 RiBi;n(u) 0 u 1 Ri and Bi;n(u) represent polygon vertices and the Bernstein polynomial basis functions respec- tively. The de nition of a Bernstein polynomial is: Bi;n(u) = n i ! ui(1 u)n i i = 0;1;2;:::n where n i ! = n!i!(n i)! The polygon joining R0;R1;:::;Rn is called the control polygon. Examples of B ezier curves: n=2: Quadratic B ezier Curves (Parabola), see Figures 4.6 and Figure 4.7 R(u) = R0 (1 u)2 + R1 2u(1 u) + R2 u2 = h u2 u 1 i 2 64 1 2 1 2 2 0 1 0 0 3 75 2 64 R0R 1 R2 3 75 7 0.0 0.0 1.0 1.0 i=0 u Bi;2(u) i=2 i=1 Figure 4.6: Plot of the quadratic Bernstein basis functions. u = 1 u = 0 R 2 R1 R0 Figure 4.7: Illustration of end conditions for a quadratic B ezier curve. n=3: Cubic B ezier Curves, see Figures 4.8 and 4.9 R(u) = R0 (1 u)3 + R1 3u(1 u)2 + R2 3u2(1 u) + R3 u3 = h u3 u2 u 1 i 2 66 64 1 3 3 1 3 6 3 0 3 3 0 0 1 0 0 0 3 77 75 2 66 64 R0 R1 R2 R3 3 77 75 Interpolation of R0, Rn and tangency of curve with polygon at u = 0;1: The B ezier curve approximates and smooths control polygon. Properties of Bernstein polynomials 1. Positivity Bi;n(u) 0 in 0 u 1 8 0.0 0.0 1.0 1.0 u i=3i=2i=1i=0 Bi;3(u) Figure 4.8: Plot of the cubic Bernstein basis functions. 2. Partition of unity Pni=0 Bi;n(u) = [1 u + u]n = 1 (by the binomial theorem) 3. Recursion Bi;n(u) = (1 u)Bi;n 1(u) + uBi 1;n 1(u) with Bi;n(u) = 0 for i < 0; i > n and B0;0(u) = 1 4. Linear Precision Property u = nX i=0 i nBi;n(u) Conversion of explicit curves to parametric B ezier curves. Parametrization of straight lines y = y(x) ! x = u;y = y(u) 5. Degree elevation: The basis functions of degree n can be expressed in terms of those of degree n + 1 as: Bi;n(u) = 1 in + 1 Bi;n+1(u) + i + 1n + 1Bi+1;n+1(u); i = 0;1; ;n 6. Symmetry Bi;n(u) = Bn i;n(1 u) 7. Derivative B0i;n(u) = n[Bi 1;n 1(u) Bi;n 1(u)] where B 1;n 1(u) = Bn;n 1(u) = 0 8. Basis conversion (for each of the x;y;z; coordinates) Ferguson or monomial form: f(u) = 3X i=0 fiui = h u3 u2 u 1 i 2 66 64 f3 f2 f1 f0 3 77 75 = U FM 9 R0 X R3 R2 R1 R0 R1 R3 R2 R1 R0 R3 R2 R3 R2 R1 R0 X Figure 4.9: Examples of cubic B ezier curves approximating their control polygons. Hermite form: f(u) = h u3 u2 u 1 i MH 2 66 64 f(0) f(1) _f(0) _f(1) 3 77 75 = U MH FH B ezier form: f(u) = h u3 u2 u 1 i MB 2 66 64 F0 F1 F2 F3 3 77 75 = U MB FB Conversion between forms: U FM = U MH FH = U MB FB ) FM = MH FH = MB FB or ) FH = M 1H FM = M 1H MB FB etc. 10 Properties of B ezier curve: 1. Geometric conditions at ends R(0) = R0; R(1) = Rn _R(0) = n(R1 R0); _R(1) = n(Rn Rn 1) R(0) = n(n 1)[(R2 R1) (R1 R0)] R(1) = n(n 1)[(Rn Rn 1) (Rn 1 Rn 2)] = 1 = j _R Rj j _Rj3 2. Curve symmetry R i = Rn i () R (u) = R(1 u), see Figure 4.10 R1 = R 2 R2 = R 1 R3 = R 0 u = 0 u = 1 R0 = R 3 Figure 4.10: B ezier curve symmetry 3. Bounding box min(xi) Pni=0 xiBi;n(u) max(xi), etc, see Figure 4.11. A tighter bounding box can usually be obtained by axis reorientation. Bounding boxes are useful in interference and intersection problems in CAD. 4. Geometry invariance properties (a) Translation: Recalling that Pni=0 Bi;n(u) = 1, we can translate the whole curve by R0 and trans- late each control vertex back by R0 and show that our resulting curve is the same as the original, see Figure 4.12 R(u) = R0 + nX i=0 (Ri R0)Bi;n(u) = nX i=0 RiBi;n(u) Absolute position of origin is unimportant ) \invariance under translation ". Relative position of control vertices is important. 11 min(xi) max(xi) min(yi) max(yi) y x Figure 4.11: B ezier curve bounding box R0 R1 R2 R3 O Figure 4.12: Geometry of B ezier curves is invariant under translation. (b) Rotation (using a rotation matrix A): Let A be a matrix of rotation about an arbitrary point. We can then demonstrate that the B ezier curve is invariant under rotation: R(u) = [x(u);y(u);z(u)] Ri = [xi;yi;zi] [x0(u);y0(u);z0(u)] = [x(u);y(u);z(u)] A = nX i=0 [xi;yi;zi]ABi;n(u) = nX i=0 [x0i;y0i;z0i]Bi;n(u) Form of B ezier curve is invariant under rotations Curve is easily computable Therefore, B ezier curves are invariant under translations and rotations and they are easily computable under such transformations. 5. End points geometric property: 12 The rst and last control points are the endpoints of the curve. In other words, b0 = R(0) and bn = R(1). The curve is tangent to the control polygon at the endpoints. This can be easily observed by taking the rst derivative of a B ezier curve _R(u) = dR(u) du = n n 1X i=0 (bi+1 bi)Bi;n 1(u); 0 u 1: In particular we have _R(0) = n(b1 b0) and _R(1) = n(bn bn 1). Equation (4.5) can be simpli ed by setting bi = bi+1 bi: _R(u) = n n 1X i=0 biBi;n 1(u); 0 u 1 The rst derivative of a B ezier curve, which is called hodograph, is an another B ezier curve whose degree is lower than the original curve by one and has control points bi, i = 0; ;n 1. Hodograph is useful in the study of intersection and other interrogation problems such as singularities and in ection points. 6. Convex hull property R(u) = nX i=0 RiBi;n(u); 0 u 1 Recall: 0 Bi;n(u) 1 and Pni=0 Bi;n(u) = 1. n=1 The convex hull is the line segment connecting the two vertices, see Figure 4.13 Convex Hull R1R0 = Figure 4.13: Convex hull of a linear B ezier curve. R(u) = (1 u)R0 + uR1; 0 u 1 = R0 + (R1 R0)(u); 0 u 1 n=2 R(u) = R0(1 u)2 + 2u(1 u)R1 + u2R2 0 u 1 = R0 + 2u(1 u)(R1 R0) + u2(R2 R0) = R0 + (R1 R0) + (R2 R0) Here we have made the substitution = 2u(1 u) and = u2 to de ne an oblique coordinate system ( , ). With the restriction 0 u 1 ) 0 ; 1, and + = 1 (1 u)2 1, any point on the curve must also be contained in the triangle de ned by the polygon vertices, as shown in Figure 4.14. 13 + b < 1 R0 R1 R2 + = 1 = 0 = 0 = 1 = 1 Figure 4.14: Convex hull of a quadratic B ezier curve. R0 R1 R2 R3 + + = 1 Figure 4.15: Convex hull of a 3D space cubic B ezier curve is a tetrahedron. n=3 R(u) = R0 + (R1 R0)+ (R2 R0)+ (R3 R0) 0 ; ; 1 for 0 u 1 ; ; and are the oblique coordinates of a point on the curve. With 0 ; ; 1, and + + = 1 (1 u)3 1 the curve is contained in the tetrahedron R0R1R2R3, see Figure 4.15 or the quadrilateral in Figure 4.16. In general, the convex hull is the intersection of all the convex sets containing all vertices or the intersection of the half spaces generated by taking 3 vertices at a time to construct a plane and having all other vertices on one side. The convex hull can also be conceptualized at the shape of a rubber band/sheet stretched taut over the polygon vertices. This is illustrated for the 2D case in Figure 4.16. The convex hull property is useful in Intersection problems: { Subdivide curves into two curve segments at u = 0:5. { Compare the convex hulls of the subdivided curves (see Figure 4.17). { Continue subdivision of intersecting subcurves until the size of their convex hulls are less than a given tolerance. 14 R0 R1 R2 R3 Figure 4.16: Convex hull of a planar cubic B ezier curve Detection of absence of interference. Providing the approximate position of curve through simple bounds. Figure 4.17: Comparison of convex hulls of B ezier curves as a means to detect intersection 7. Variation diminishing property A polygon can be created by the segments connecting the ordered vertices of a B ezier curve. In 2D: The number of intersections of a straight line with a planar B ezier curve is no greater than the number of intersections of the line with the control polygon. A line intersecting the convex hull of a planar B ezier curve may intersect the curve, be tangent to the curve, or not intersect the curve at all. It may not, however, intersect the curve more times than it intersects the control polygon. This property is illustrated in Figure 4.18. In 3D: The same relation holds true for a plane. Results : (a rough interpretation) A B ezier curve oscillates less than its polygon. The polygon’s segments exaggerate the oscillation of the curve. This property is important in intersection algorithms and in detecting the \fairness" of B ezier curves. 15 ImpossiblePossible Figure 4.18: Variation diminishing property of a cubic B ezier curve Algorithms for B ezier curves Evaluation and subdivision algorithm: A B ezier curve can be evaluated at a speci c pa- rameter value t0 and the curve can be split at that value using the de Casteljau algorithm, where bki (t0) = (1 t0)bk 1i 1 + t0bk 1i ; k = 1;2;::: ;n; i = k;::: ;n (4.5) is applied recursively to obtain the new control points. The algorithm is illustrated in Figure 4.19, and has the following properties: { The values b0i are the original control points of the curve. { The value of the curve at parameter value t0 is bnn. { The curve can be represented as two curves, with control points (b00; b11;::: ; bnn) and (bnn; bn 1n ;::: ; b0n). Continuity algorithm: B ezier curves can represent complex curves by increasing the degree and thus the number of control points. Alternatively, complex curves can be represented using composite curves, which can be formed by joining several B ezier curves end to end. If this method is adopted, the continuity between consecutive curves must be addressed. One set of continuity conditions are the geometric continuity conditions, designated by the letter G with an integer exponent. Position continuity, or G0 continuity, requires the endpoints of the two curves to coincide, Ra(1) = Rb(0): The superscripts denote the rst and second curves. Tangent continuity, or G1 continuity, requires G0 continuity and in addition the tangents of the curves to be in the same direction, _Ra(1) = 1t _Rb(0) = 2t 16 where t is the common unit tangent vector and 1, 2 are the magnitude of _ra(1) and _rb(0). G1 continuity is important in minimizing stress concentrations in physical solids and preventing ow separation in uids. Curvature continuity, or G2 continuity, requires G1 continuity and in addition the center of curvature to move continuously past the connection point , Rb(0) = 2 1 2 Ra(1) + _Ra(1): where is an arbitrary constant. G2 continuity is important for aesthetic reasons and also for preventing uid ow separation. More stringent continuity conditions are the parametric continuity conditions, where C k continuity requires the kth derivative (and all lower derivatives) of each curve to be equal at the joining point. In other words, dkRa(1) dtk = dkRb(0) dtk : The C1 and C2 continuity conditions for consecutive segments of a composite degree n B ezier curve can be stated as hi+1 (bni bni 1) = hi (bni+1 bni) ; and (4.6) bni 1 + hi+1h i (bni 1 bni 2) = bni+1 + hih i+1 (bni+1 bni+2) (4.7) where, for the ith B ezier curve segment parameter t runs over the interval [ti, ti+1], hi = ti+1 ti (see Figure 4.20 for the connection of cubic B ezier curve segments). Degree elevation: The degree elevation algorithm permits us to increase the degree and control points of a B ezier curve from n to n+1 without changing the shape of the curve. The new control points bi of the degree n + 1 curve are given by bi = i n + 1bi 1 + 1 in + 1 bi; i = 0;::: ;n + 1 (4.8) where b 1 = bn+1 = 0 ls 4.3.5 B ezier surfaces A tensor product surface is formed by moving a curve through space while allowing deformations in that curve. This can be thought of as allowing each control point bi to sweep a curve in space. If this surface is represented using Bernstein polynomials, a B ezier surface (patch) is formed, with the following formula: R(u;v) = mX i=0 nX j=0 bijBi;m(u)Bj;n(v); 0 u;v 1: 17 Here, the set of straight lines drawn between consecutive control points bij is referred to as the control net. It is easy to see that boundary isoparametric curves (u = 0, u = 1, v = 0 and v = 1) have the same control points as the corresponding boundary points on the net. An example of a bi-quadratic B ezier surface with its control net can be seen in Figure 4.21. Since a B ezier surface is a direct extension of univariate B ezier curve to its bivariate form, it inherits many of the properties of the B ezier curve described before such as: Geometry invariance property End points geometric property Convex hull property However, no variation diminishing property is known for B ezier surface patch. 4.4 Composite curves 4.4.1 Motivation Increase complexity without increasing degree of polynomial. Interpolation/ approximation of large arrays of points. Avoid problems of high degree polynomials such as oscillation, see Section 4.4.2. 4.4.2 Lagrange basis The Lagrange form is another basis for curve representation. It is presented here as motivation for establishing a composite polynomial curve representation of lower order. Given: (ui; Ri) i = 0;1;2; ;N, we create an interpolatory curve in terms of: R(u) = NX i=0 RiLi;N(u) 0 u 1 where Li;N(uj) = ij and Li;N(u) = Nj=0;j6=i u uj ui uj : where ij = 1 if i = j and ij = 0 if i 6= j ( ij is called the Kronecker delta); ui are parameter values (eq. accumulated length of the polygon R0R1:::Ri) and Ri are data points; Li;N is a polynomial of degree N. As with any polynomial of degree N, it has up to N 1 real extrema. As can be seen in Figure 4.22, Lagrange interpolation may lead to highly oscillatory approximants. The Lagrange basis is not frequently used because of oscillation and of the need to use high degree basis functions for interpolation. This leads to developing methods involving low degree composite curves (splines). 18 4.4.3 Continuity conditions See Figure 4.23 for an illustration. Position: R(1)(1) = R(2)(0) Tangent: _R(1)(1) = 1t _R(2)(0) = 2t where t is a common unit tangent vector at the interface and 1 and 2 are real constants. Curvature: _R = _st where _s = j _Rj _t = ts _s = _s n R = st + _s2 n Using the above equations, taking _R R and using t n = b implies b = _R R j _Rj3 Therefore, if ; n are continuous ! b = t n is also continuous ! b is continuous. Hence, _R(1)(1) R(1)(1) j _R(1)(1)j3 = _R(2)(0) R(2)(0) j _R(2)(0)j3 t h R(2)(0) 2 1 2 R(1)(1) i = 0 R(2)(0) = 2 1 2 R(1)(1) + _R(1)(1) where is another real constant. Order of geometric continuity: G0{ position continuity, G1{ unit tangent continuity, and G2{ center of curvature continuity. Order of parametric continuity (more restrictive than geometric continuity): C0{ position continuity, C1{ rst derivative continuity, and C2{ second derivative continuity. 19 4.4.4 B ezier composite curves (splines) Cardinal or interpolatory splines are used for tting data. B ezier composite curves usually used for preliminary design (without tting) R(u) = 3X i=0 ui(1 u)3 iRi 0 u 1 Conditions for position and tangent continuity, see Figure 4.24 { Position (G0) R(1)3 = R(2)0 { Tangent (G1) t = 3 1 R(1)3 R(1)2 = 3 2 R(2)1 R(2)0 Hence, R(1)3 = R(2)0 and R(1)2 , R(2)1 are collinear. { Conditions for curvature continuity: R(2)(0) = 2 R(1)(1) + _R(1)(1) = 2 1 Using R(1)(1) = 6 R(1)3 2R(1)2 + R(1)1 R(2)2 R(1)3 = 2 R(1)1 R(1)2 + 2 + 2 + 2 R(1)3 R(1)2 These conditions imply coplanarity of points R(1)1 , R(1)2 , R(1)3 , R(2)1 , R(2)2 . The design process may involve the following steps: Start with one segment. For each new segment ; ;R3 may be chosen freely for shape creation. However, such construction may lead to not enough freedom for G2 continuity. enough freedom for G1 continuity. If more degrees of freedom are needed, an increase of the degree of the curve is necessary. 20 4.4.5 Uniform Cubic B-splines Motivation: Changing of one polygon vertex of a B ezier curve or one data point of a cardinal or interpolatory spline a ects the curve GLOBALLY. Avoid increasing degree of B ezier curve to construct complex composite shapes. Develop a LOCAL curve scheme. De ne a uniform, cubic B-spline curve span by: Ri(u) = ViN0;4(u) + Vi+1N1;4(u) + Vi+2N2;4(u) + Vi+3N3;4(u) Ri(u) represents the ith span of the composite curve and Ni;4(u) represents the ith blending function of order 4 (or degree 3). Each span is de ned in terms of a local parameter: 0 u 1. Nk;4(u) are degree 3 polynomials de ned from continuity conditions. Let us assume: Nj;4(u) = aj + bju + cju2 + dju3 j = 0; ;3 Consider two adjacent spans, as shown in Figure 4.25 with control points Vi shown in Figure 4.26. For C2 continuity, the following three vector conditions must hold: Ri(1) = Ri+1(0) ) h N0;4(1) N1;4(1) N2;4(1) N3;4(1) i 2 66 64 Vi Vi+1 Vi+2 Vi+3 3 77 75 = h N0;4(0) N1;4(0) N2;4(0) N3;4(0) i 2 66 64 Vi+1 Vi+2 Vi+3 Vi+4 3 77 75 This is valid for every Vi ) N0;4(1) = 0 (4.9) N1;4(1) = N0;4(0) (4.10) N2;4(1) = N1;4(0) (4.11) N3;4(1) = N2;4(0) (4.12) N3;4(0) = 0 (4.13) R0i(1) = R0i+1(0) ) N00;4(1) = 0 (4.14) N01;4(1) = N00;4(0) (4.15) N02;4(1) = N01;4(0) (4.16) N03;4(1) = N02;4(0) (4.17) N03;4(0) = 0 (4.18) 21 R00i (1) = R00i+1(0) ) N000;4(1) = 0 (4.19) N001;4(1) = N000;4(0) (4.20) N002;4(1) = N001;4(0) (4.21) N003;4(1) = N002;4(0) (4.22) N003;4(0) = 0 (4.23) Also if all four vertices coincide ) the ith span reduces to a point ) Ri(u) = Vi 3X i=0 Ni;4(u) = 1 (4.24) Solving the above 16 equations 4.9{4.24 yield aj;bj;cj;dj; j = 0;1;2;3. h N0;4(u) N1;4(u) N2;4(u) N3;4(u) i = 13! h 1 u u2 u3 i MB = 13! h 1 u u2 u3 i 2 66 64 1 4 1 0 3 0 3 0 3 6 3 0 1 3 3 1 3 77 75 Properties: Position, derivatives (see Figure 4.26): Ri(0) = 23Vi+1 + 13 Vi + Vi+2 2 Ri(1) = 23Vi+2 + 13 Vi+1 + Vi+3 2 R0i(0) = 12 Vi+2 Vi R0i(1) = 12 Vi+3 Vi+1 R00i (0) = Vi+2 Vi+1 Vi+1 Vi R00i (1) = Vi+3 Vi+2 Vi+2 Vi+1 Local control { one vertex a ects only 4 spans. Convex hull property { local. Variation diminishing property. Replication of polygon control of B ezier curves at end points can be achieved by setting V 1 = 2V0 V1 (imaginary vertex) Then R 1(0) = 23V0 + 16(V 1 + V1) = V0 R0 1(0) = 12(V1 V 1) = V1 V0 22 Interpolation can also be achieved, via the so called cardinal or interpolatory splines, which involve interpolation of N data points Ri (i = 0;1;:::N 1) 2 3Vi+1 + 1 6(Vi + Vi+2) = Ri or Vi + 4Vi+1 + Vi+2 = 6Ri for i = 0; 2; N 1 We have N + 2 unknowns so we need two boundary conditions. To replicate the natural spline boundary conditions, we set R00(0) = R00(1) = 0, which imply zero curvature at the ends. Natural splines minimize the integral of curvature square and are the smoothest curves that pass through a set of points. We then have R000(0) = V2 2V1 +V0 = 0: Also we have N 1 spans, N data points and N +2 control points, V0;V1;:::VN+1 the condition will be R00N 1(0) = VN+1 2VN + VN 1 = 0: Then we can derive the banded matrix to solve the above linear system. 2 66 66 66 66 66 4 1 2 1 1 4 1 1 4 1 1 4 1 ... 1 4 1 1 2 1 3 77 77 77 77 77 5 2 66 66 66 66 66 4 V0 V1 V2 V3 ... VN VN+1 3 77 77 77 77 77 5 = 6 2 66 66 66 66 66 4 0 R0 R1 R2 ... RN+1 0 3 77 77 77 77 77 5 23 b03 =b02 =b01 =b00 b10 b20 b30 t 1-t t 1-t t 1-t b13 =b12 =b11 b21 b31b2 2 b32 b33 b00 b10 b20 b30 b11 b21 b31 b22 b32 1?t t 1?t 1?t 1?t 1?tt t t t b3=r(t)3 1?t t Figure 4.19: The de Casteljau algorithm. 24 bn i-3 bn i+3 bn i-2 bn i+2 bn i-1 bn i+1 bn ih i : hi+1 hi : hi+1 hi : hi+1 Figure 4.20: Continuity conditions. b20 b21 b22 b10 b12b00 b01 b02 Figure 4.21: A B ezier Surface with Control Net. 25 Exact Lagrange Approximation 0 5 10 15 20 25 30?4 ?3 ?2 ?1 0 1 2 3 4 5 6 x y y = x^(1/3) Figure 4.22: Approximation of the function y = x 13 using a cubic Lagrange curve. 0 u2 1 0 u1 1 R(2)(u2) R(1)(u1) Figure 4.23: Composite curve consisting of two smaller curve segments. R(1)1 R(1)2 R(1)3 = R(2)0 R(1)0 R(2) 3 R(2)2 R(2)1 Figure 4.24: Continuity between B ezier curves restricts placement of control vertices. 26 0 u 1 0 u 1 Ri+1( u) Ri(u) Figure 4.25: Two adjacent curve spans with position, rst, and second derivative continuity. Vi Vi+1 Vi+2 Ri+1(u)Ri(u) Vi+4 Vi+3 Figure 4.26: Two spans of a uniform cubic B-spline curve. 27 Bibliography [1] G. Farin. Curves and Surfaces for Computer Aided Geometric Design: A Practical Guide. Academic Press, Boston, MA, 3rd edition, 1993. [2] I. D. Faux and M. J. Pratt. Computational Geometry for Design and Manufacture. Ellis Horwood, Chichester, England, 1981. [3] J. Hoschek and D. Lasser. Fundamentals of Computer Aided Geometric Design. A. K. Pe- ters, Wellesley, MA, 1993. Translated by L. L. Schumaker. [4] L. A. Piegl and W. Tiller. The NURBS Book. Springer, New York, 1995. [5] F. Yamaguchi. Curves and Surfaces in Computer Aided Geometric Design. Springer-Verlag, NY, 1988. 28