- Email: [email protected]
Extensions of the general polar value based control point specification method in constructing tensor product B-spline surfaces
Computers & Graphics 24 (2000) 493}507
Technical Section
Extensions of the general polar value based control point speci"cation method in constructing tensor product B-spline surfaces Stephen Wang-Cheung Lam* Department of Computing, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong
Abstract Blossoming has been utilized as a powerful tool for understanding Bezier and B-spline curves. Recently, general polar values (GPV) have been employed as an aid to specify control points for polynomial curves. Based on this technique, it is possible for us to use a set of polar values to control a polynomial curve. In this paper, we "rst construct linear systems of tensor product B-spline surfaces (TPBS) under the GPV-based framework. We call this extension as the direct extension scheme. Subsequently, based on similar principles, we investigate the issue of extending the linear systems based on degree elevation and partial derivatives evaluation. For the direct extension and degree elevation schemes, high #exibility is provided though a large set of corresponding free parameters. For the partial derivative scheme, partial derivatives on a surface can be used as parameters. In all schemes, a restriction on parameter selection is imposed by the invertibility of the related matrices. Finally, based on the subtle relationships among the basis functions, a fast computation algorithm is developed for evaluating those related basis-function matrices. 2000 Elsevier Science Ltd. All rights reserved. Keywords: General polar value; Control point; Tensor product surface; Partial derivative; Degree elevation; Basis-function matrices
1. Introduction It is well known that the polar form of a degree n polynomial G(u) is a symmetric, multia$ne polynomial g(u ,2, u ) that satis"es the identity L g( t,..., t )"G(t). GHI L In [1], this principle is applied to the formulation of Bezier and B-spline curves and surfaces that are frequently used in computer-aided geometric design. The word `blossominga is used to name the process of replacement of a polynomial of high degree with its polar form. The main bene"t of that replacement is the dis-
* Tel.: #851-2766-7301; fax: #851-2774-0842. E-mail address: [email protected] (S. WangCheung Lam).
covery of a way to label the Bezier points that control a curve segment or surface patch and the de Boor points that control a B-spline curve, with symmetric, multivariate labels. This labeling gives us great geometric insight into the basic algorithms in this area, such as the de Casteljau algorithm and the de Boor algorithm. A good introduction to this topic can be found in [2,3] and an in-depth treatment of the theory behind blossoming is given in [4]. Blossoming has been applied in developing geometric algorithms and proo"ng existing algorithms with new insights. For example, in [5], the theory of blossoming has been applied to the insertion of a new knot into the knot vector of a B-spline curve. There, a formula that expresses the de Boor points of the re"ned curve in terms of a set of polar values is derived. As a corollary, a new proof of Boehm's knot insertion algorithm has also been obtained. Furthermore, the relationship between multiple knot insertion and the de Boor algorithm is also revealed through blossoming. To extend blossoming to
0097-8493/00/$ - see front matter 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 7 - 8 4 9 3 ( 0 0 ) 0 0 0 5 2 - 2
494
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
B-spline, in [6,7], Seidel attempts to apply blossoming directly to the recurrence relation of B-splines. There, he gives a new derivation for expressing the de Boor points of a B-spline curve in terms of blossoming. Based on that, the Curry}Schoenberg theorem and the de Boor algorithm were re-derived. Although all major theorems in that article are already known, the proofs are shorter than the previous ones and they also give further insights into the basic theory of B-splines construction. Recently, blossoming has been applied to the formulation of triangular B-splines surfaces [8,9] and geometric continuous spline curves [10], computation of product of splines [11] and composition of Bezier simplexes [12,13]. Besides, extensions of the technique to multirational framework can be found in [14,15]. In [16], Duchaineau shows that general polar values (GPV) may also be used to control a polynomial curve when a related matrix is invertible. In that case, those GPV act as blossom control points and the invertible matrix provides a useful translation from these blossom control points to well known ones such as those of Bezier. In [16], the special case in which the polar values can be evaluated through pairwise a$ne combination is also characterized. That characterization provides users of this approach a way of generating GPV. As stated in [16], this GPV approach allows the arguments of the blossom control points to be chosen in a manner akin to choosing the knot vector of a B-spline segment. Furthermore, the main advantage of this approach over the usual B-spline control points is the increased number of free parameters for specifying the arguments of the blossom control points. As stated in [16], one of the important future research works is on the extension of the technique to the area of surfaces. We agree with him and hence in this paper, we extend the GPV approach to surfaces. Besides, we also investigate other possible ways of applying GPV in specifying control points. To ful"ll our goal, we "rst construct the linear system of a tensor product B-spline surface (TPBS) under the GPV-based framework. Subsequently, we extend the system based on degree elevation and partial derivatives calculation. Finally, after a close examination of the calculation of the related basis function matrices, we discover that it is possible to greatly reduce the computation time for evaluating those matrices through using a smarter computing algorithm. The remainder of this paper is structured as follows. In Section 2, we present a brief review of the blossoming of Bezier and B-spline curves. In Section 3, we describe the fundamental concepts and techniques used in general polar values (GPV) based control point speci"cation. In Section 4, we "rst apply GPV in controlling B-spline curve. Then, we extend it to tensor product B-spline surfaces (TPBS). In particular, we focus on controlling an individual region on a surface. We call this extension as direct extension scheme. Moreover, we also investigate
other usage of the related basis function matrix. In Section 5, we modify the GPV system based on degree elevation. Similar to the original scheme [16], the direct extension and degree elevation schemes provide high #exibility by allowing a large set of free parameters. In Section 6, we extend our system through using partial derivative formulation. Furthermore, we study the issue of applying higher-order derivative calculation in controlling a TPBS. In all schemes, a restriction of parameter selection is imposed by the invertibility of the related matrices. In Section 7, a fast computation algorithm is developed for calculating the basis function matrices used in our system. In Section 8, we perform experiments to illustrate the schemes we proposed in controlling a TPBS. In Section 9, we discuss the pros and cons of our approach. Finally, in Section 10, we summarize our contributions and suggest some topics for future research works.
2. Review on blossoming This section aims at brie#y reviewing some of the results in blossoming theory of Bezier and B-spline curves [5]. 2.1. Basis De5nition. A map f : 1P1:(u ,2, u )Pf (u ,2, u ) is L L n-a$ne, if it is a$ne in each argument. In other words, the following formula f (u ,2 a t ,2, u )" a f (u ,2, t ,2, u ) (2.1) G G L G G L G G holds for all a , t 31 with a "1. Furthermore, f is G G G G symmetric, if it keeps itsvalue under any permutation of its arguments. Blossoming Principle: A one-to-one correspondence exists between polynomial functions F : 1P1B of degree n and symmetric n-a$ne maps f : 1P1B. Further, given a function of either type, a unique function of the other type exists that satis"es the identity F(u)"f (u,2,u), where f is called the blossom of F and F is the diagonal of f. 2.2. Bezier and B-spline control points If F is taken as a Bezier curve over [s t], its control points b ,2, b are given by the formula L b "f ( s,..., s, t,..., t ). (2.2) G GHI GHI L\G G Obviously, for spline curves the situation becomes more complicated. The following formula represents
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
a non-periodic B-spline curve K F(t)" d ) N (t) G GI G
where (2.3)
of order k (i.e. degree n"k!1) over the non-decreasing knot sequence ¹"(t "t "2"t , t ,2t , t "2"t ) I\ I K K> K>I that does not contain any intermediate knots t with G t )t )t of multiplicity greater than k!1. In Eq. I G K (2.3), d ,.., d are the de Boor points (B-spline control K points) and N are the normalized B-splines of order GI k over T. For t (t , we denote by F the polynomial of the H H> H B-spline curve F in the segment (t , t ). F is of order H H> H k with blossom f . For every i)j)i#k!1 with H t (t the de Boor point d is related to the blossom H H> G f by the formula H d "f (t ,2, t ) G H G> G>I\
(2.4)
Furthermore, since every knot t appears with multipliG city at most k, such a j always exists. After a brief review of the blossoming theory, we will introduce the GPV-based control point speci"cation scheme in the next section. Readers who are interested in other fundamental theorems and results obtained in blossoming theory can refer to [1,4,5].
3. General polar values In [16], Duchaineau shows that general polar values (GPV) may also be used to control polynomial curves when a related matrix is invertible. In this section, we brie#y review this methodology. Any nth-degree polynomial curve F(x) can be written as L F(x)" Q B (x) G G G
(3.1)
where Q ,2, Q are control points and B (x),2, B (x) L L are any known basis for the nth-degree polynomials. One example is the Bezier basis. The polar form f of F can be written as L f (x ,2, x )" Q b (x ,2, x ) L G G L G
(3.2)
where b ,2, b are the polar forms of the basis functions. L Obviously, the set of polar values P can be obtained by the application of a linear transformation A to the Q control points, that is P"AQ
495
(3.3)
P f (t ,2, t ) L P f (t ,2, t ) L P" " $ $ P L
and
A"
(3.4)
f (t ,2, t ) L LL
b (t ,2, t ) 2 b (t ,2, t ) L L L b (t ,2, t ) 2 b (t ,2, t ) L L L . $ b (t ,2, t ) L LL
2
b (t ,2, t ) L L LL
(3.5)
It is interesting to know whether a set of n#1 polar values with general argument suite is capable of controling the nth-degree polynomial curves. In [16], the word complete is used to describe the polar values [ f (t ,2, t ),2, f (t ,2, t )] when a polynomial can L L LL be produced by setting P appropriately. This occurs precisely when the matrix A is invertible, since this allows Q to be obtained from P using Q"A\P. That implies there is a one-to-one correspondence between P and Q. The following provides the de"nition of completeness. De5nition 2.1 (Duchainean [16]). A set of n#1 polar values [f (t ,2, t ),2, f (t ,2, t )] is complete if and L L LL only if the linear transformation to these polar values form control points in any known basis of the nth-degree polynomials has a nonzero determinant. If the polar values have been proved to be complete, they will be useful as control points and shall therefore be called blossom control points [16].
4. GPV based B-spline system From [6], F (u)" H d ) NN (u) represents the H GH\N G G H polynomial in the non-empty interval I :"[t , t ) in H H H> a B-spline curve F(u), where t (t . The corresponding H H> polar form of F (u) is expressed as follows: H H f (u ,2, u )" d ) gN (u ,2, u ). (4.1) H N G G H N GH\N Based on this formula, we construct the linear system for B-spline curves. Following the notation in [16], we de"ne the GPV linear transformation of a B-spline curve as follows: P"AQ.
(4.2)
The entries of P is f (;) where and f is the polar form H H of the polynomial starting from the knot t and U is an H argument suite. Q represents the control points. An entry
496
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
in A is the polar form of a B-spline basis function. We denote it as gN (;) where p is the degree of the spline. GH In general, we have
f (; ) H f (; ) H $
gN (; ) H\N H gN (; ) H\N H $
gN (; ) H\N> H gN (; ) H\N> H $
gN (; ) H\N H N>
gN (; ) H\N> H N>
... ... ...
gN (; ) H H gN (; ) H H $
... gN (; ) H H N>
(4.3)
d H
For example, given a normalized cubic B-spline curve with knot vector ¹"(0, 0, 0, 0, 1, 2, 4, 5, 6, 6, 6, 6) and control point vector Q"(d , d , d , d , d , d , d , d ) [5], when j"5 and j"6, according to (4.3) we have, respectively, f (; ) g (; ) f (; ) g (; ) " f (; ) g (; ) f (; ) g (; )
;
and
d
d
d
d
f (; ) g (; ) f (; ) g (; ) " f (; ) g (; ) f (; ) g (; )
;
De5nition 4.1. A set of n#1 polar values + f (; ),2, f (; ),, where f is the polar form of the H H > H polynomial between t and t and n is the number of H H> control points, is complete if the linear transformation to these polar values form the corresponding set of control points of a normalized B-spline curve has a nonzero determinant. 4.1. Other usage of the transform
d H\N> . $
(4.6)
Hence, similar to [16], here completeness is further de"ned as follows:
d H\N
;
If A is invertible, we can obtain the control points Q of a B-spline curve segment through the following equation: Q"A\P
f (; ) H N>
"
where ; is a suite of arguments. Notice that ; in Eqs. G G (4.4) and (4.5) may be di!erent.
d
d
d
d
g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; )
(4.4)
Similar to [16], the matrix A given above turns out to be useful for more than just testing for the completeness of the set of polar values. With Bezier basis, matrix A is shown to be useful for Bezier curve composition based on derivatives [16]. Here, we suggest that matrix A can also be used in de"ning new GPV-based blossom control points in case of knot insertion. In particular, after knot insertion, new control points Q* are obtained from the old ones Q through a linear transformation (p. 262 of [17]): QH"BQ
(4.7)
Based on (4.6), the new GPV vector can therefore be expressed as P*"(A*BA\)P.
(4.8)
This equation is now in the same form of Eq. (4.7) and, similar to the ordinary control points, we can then obtain the new set of blossom control points from the old ones. 4.2. Tensor product surface
g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; )
(4.5)
We "rst review the de"nition and some properties of tensor product surfaces. Most of the materials in this review are taken from [1]. Let F : 1;1PO be a polynomial function; that is, each coordinate of the point F(u, v) is given by a polynomial in the two real valuables u and v. In a (m; n) degree tensor product surface, the coordinates of F(u, v) are required to have degree no higher than m in the variable u and degree no higher than n in the variable v separately. In a triangular patch surface of degree n, the coordinates of F(u, v) are required to have total degree no higher than n in the variables u and v jointly. In this paper, we will not consider triangular patch surfaces.
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
The blossom f of a tensor product surface F of degree (m; n) is an (m#n) a$ne function f : 1 ;1PO that is symmetric in its "rst m arguments and its last n arguments. f also satis"es the identity F(u, v)" f (u,2, u;v,2, v). Also, every tensor product surface has a unique blossom. If [p, q];[r, s] is a reference rectangle for the domain plane 1;1, the (m#1)(n#1) Bezier points of the surface F are the blossom value f ( p,2, p, q,2, q, r,2, r, s,2, s ) (4.9) GHI GHI GHI GHI K\G L\H G H for i in [0, m] and j in [0, n]. The de Casteljau algorithm for a tensor product surface calculates the blossom value f (u ,2, u ; v ,2, v ) by performing m#n stages of lin K L ear interpolations, with each blossom argument providing the ratio for one of these stages. Moreover, these stages can be performed in any order. As stated in [18], a tensor product surface F(u, v) can be evaluated as follows: K L F(u, v)" d NK(u)NL (v). VW V W V W The above formula can be rewritten as
(4.10)
L K (4.11) F(u, v)" d NL (v) NK(u). VW W V V W Based on (4.1), we derive the polar form of the above formula as the following equation: f (u ,2, u ; v ,2, v ) GH K L H G d gL (v ,2, v ) gK (u ,2, u ). " VW W H L V G K VG\K WH\L (4.12)
For simplicity, we just concentrate on the case where i"j in the following parts of this paper and (4.12) is therefore rewritten as f (u ,2, u ; v ,2, v ) H K L H H " d gL (v ,2, v ) gK (u ,2, u ). VW W H L V H K VH\K WH\L (4.13)
To construct the linear system for a tensor product B-spline surface, similar to the case of B-spline curve, we write the following: P"XQ>
(4.14)
where P, X, Q and > are given by
P"
f (; , < ) H f (; , < ) H $
f (; , < ) H f (; , < ) H $
2 2 2
f (; , < ) H L> f (; , < ) H L> $
f (; , < ) f (; , < ) 2 f (; ,< ) H K> H K> H K> L>
,
X"
gK (; ) H\K H gK (; ) H\K H $
497
gK (; ) H\K> H gK (; ) H\K> H $
gK (; ) H H gK (; ) H H $
2 2 2
gK (; ) gK (; ) 2 gK (; ) H\K H K> H\K> H K> H H K>
,
d d 2 d H\K H\L H\K H\L> H\K H d d 2 d H\K> H , Q" H\K> H\L H\K> H\L> $ $ 2 $ d H H\L
d H H\L>
d H H
2
gL (< ) gL (< ) 2 gL (< ) H\L H H\L H H\L H L> gL (< ) gL (< ) 2 gL (< ) H\L> H L> . >" H\L> H H\L> H $ $ 2 $ gL (< ) H H
gL (< ) H H
gL (< ) H H L>
2
In our formulation, P and Q may not be square matrices while X and > are constructed as square matrices. In this formulation, the number of free parameters (i.e. the number of di!erent arguments in the suites) is m(m#1)#n(n#1). With some matrix arithmetic, we have the following formula Q"X\P>\.
(4.15)
This shows that if X and > are invertible, Q can be calculated from P. Suppose j"3, n"2 and m"3, P, X, Q and > are given by
f (; , < ) f (; , < ) P" f (; , < ) f (; , < )
f (; , < ) f (; , < ) f (; , < ) f (; , < )
d d Q" d d
d
d
,
d d d d
g (; ) g (; ) X" g (; ) g (; ) and
g (< ) >" g (< ) g (< )
d d
f (; , < ) f (; , < ) , f (; , < ) f (; , < )
g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; )
g (< ) g (< ) g (< )
g (< ) g (< ) . g (< )
498
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
In order to distinguish this technique with the later ones, we herein name it as direct extension scheme.
In this section, we investigate the issue of using degree elevation formulation for constructing a GPV based control point speci"cation system. Let F(u) be a polynomial curve of degree two and f (u, v) its polar form. If we consider F(u) to be a degenerate cubic polynomial curve, the blossom of F(u) is a tria$ne function, say g(u, v, w). Then, we have F(u)"f (u, u)"g(u, u, u). In [1], Ramshaw has stated the following formula: g(u, v, w)"[f (u, v)#f (v, w)#f (w, u)]/3.
(5.1)
We can understand it through revealing that righthand side of the above equation is a symmetric, tria$ne function whose diagonal agrees with F. As stated in [1], for arbitrary degree n and an arbitrary parameter space P, we have the following equation. 1 L> g(u ,2, u )" f (u ,2u( ,2, u ). (5.2) L> L n#1 G The hat over the argument u indicates that the ith G argument is omitted. Similar to the previous section, we can express a B-spline curve based on (5.2) as 1 L P"AM Q where AM " A ? n ? where P, A and Q are given by G
gL\ (;K ) gL\ (;K ) 2 gL\(;K ) H\L> H H\L> H H H gL\ (;K ) gL\ (;K ) 2 gL\(;K ) H\L> H H H A " H\L> H ? $ $ 2 $ 2
gL\(;K ) H H L
d H\L> d Q" H\L> $ d
;gK\(u ,.u( ., u ) V H J K
H
where ; denotes u ,2, u and ;K in A denotes G L G ? u ,2u ,2, u where u is omitted. Similarly, to ? L ?
(5.5)
To construct a degree-elevation-based GPV system, using distributive law, we derive the following: 1 K P"XM Q> where XM " X (5.6) ? m ? and where P, X , Q and > are given by G f (; , < ) f (; , < ) 2 f (; , < ) H H H L> f (; , < ) f (; , < ) 2 f (; , < ) H H L> , P" H $ $ 2 $
f (; , < ) H K
gK\ (;K ) H\K> H gK\ (;K ) X " H\K> H ? $
d H H\L
f (; , < ) 2 f (; , < ) H K H K L> gK\ (;K ) H\K> H gK\ (;K ) H\K> H $
gK\ (;K ) gK\ (;K ) H\K> H K H\K> H K
g (; ) H L
and
f (u ,2, u ;v ,2, v ) H K L 1 K H H " d gL (v ,.., v ) VW W H L m J VH\K> WH\L
d H\K> H\L d Q" H\K> H\L $
gL\ (;K ) H\L> H L
(5.4)
Based on Eqs. (4.13) and (5.2), if we raise the degree of a TPBS in the u direction from m!1 to m, we have
(5.3)
g (; ) H g (; ) P" H , $
gL\ (;K ) H\L> H L
Q"AM \P.
5.1. Using degree elevation in GPV based TPBS system
5. Degree-elevation-based GPV system
calculate the control points, we have
d
H\K> H\L>
d
H\K> H\L> $ d H H\L>
2 2 2
gK\(;K ) H H gK\(;K ) H H , $
2 gK\(;K ) H H K
2 d H\K> H 2 d H\K> H , 2 $ 2
d H H
gL (< ) gL (< ) 2 gL (< ) H\L H H\L H H\L H L> gL (< ) gL (< ) 2 gL (< ) H\L> H L> >" H\L> H H\L> H $ $ 2 $ gL (< ) H H
gL (< ) H H
2
gL (< ) H H L>
where ; denotes u ,2, u and ;K in X denotes G K G a u ,2, u ,2, u where u is omitted. In this case, the ? K ? number of free parameters are m#n(n#1). To evaluate the control points, we have Q"XM \P>\.
(5.7)
A similar set of equations can be derived if we choose to elevate the degree of the other direction. Now, if we want to raise the degree in both directions, we can
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
construct the following formula:
Based on Eq. (4.13), we obtain the following two formulas:
f (u ,2, u ; v ,2, v ) H K L 1 L 1 K H H " d gL\ VW W H n m P R VH\K> WH\L>
;(v ,.,v( ,., v ) gK\(u ,.u( ., u ) P L V H J K
(5.8)
.
To construct a linear system based on the above, we write the following: 1 K P"XM Q>M with XM " X ? m ?
1 L and >M " > @ n @ (5.9)
and where P and > are given by @
f (; , < ) H f (; , < ) P" H $ f (; , < ) H K
and
f (; , < ) H f (; , < ) H $
2 2
f (; , < ) H L f (; , < ) H L $
f (; , < ) 2 f (; , < ) H K H K L
gL\ ( H gL\ ( " H\L> H @ $ gL\(
2
(5.10)
gL\ ( H H\L> H L> gL\ ( H H\L> H L> $ 2 $ gL\(
2
gL\( (5.11)
< denotes v ,2, v and denotes H K H @ v ,2, v ,2, v where v is omitted. Q is then calculated @ K @ by the following formula: Q"XM \P>M \
(5.12)
f ( p,2, p; 1, q,2, q) H H H " d gL (1, q,2, q) gK ( p,2p) V H VW W H VH\K WH\L f ( p,2, p; 0, q,2, q) H
In this section, we will explore the construction of partial-derivatives-based GPV system. An extension of this system to higher-order derivatives is presented in Section 6.1. Let F(u; v) be a degree n;m tensor product surface. From [19], we obtain the following partial evaluation formula: *F (p; q) *v
For convenience, we denote f (p,2, p; 1, q,2, q)! H f (p,2, p; 0, q,2, q) as *f (p, q) and [gL (1, q,2, q)! H H W H gL (0, q,2, q)] as *gL (q). Hence, Eq. (6.4) can be simpliW H W H "ed as
H H *f (p, q)" d ) *gL (q) NK (p). (6.5) H VW W H V H VH\K WH\L To construct a GPV linear system for a TPBS using partial derivatives, we write the following: P "mX Q >
(6.6)
where
P "
*f (p , q ) H *f (p , q ) H $
*f (p , q ) H *f (p , q ) H $
(6.1)
*f (p , q ) H L> *f (p , q ) H L> $
2 2 $
,
*f (p , q ) *f (p , q ) 2 *f (p ,q ) H K> H K> H K> L> NK (p ) H\K H NK (p ) H\K H $
NK (p ) H\K> H NK (p ) H\K> H $
NK (p ) H H NK (p ) H H $
2 2 $
,
NK (p ) NK (p ) 2 NK (p ) H\K H K> H\K> H K> H H K> d H\K H\L
d H\K H\L>
2
d H H\L>
2
d
H\K H d d 2 d H\K> H , Q " H\K> H\L H\K> H\L> $ $ $ $ d
H H\L
> "
"m f ( p,..., p, 1, q,..., q ) GHI GHI L K !f ( p,..., p, 0, q,..., q ) GHI GHI L K
(6.2)
H H " d gL (0, q,2, q) gK ( p,2, p) (6.3) VW W H V H VH\K WH\L For simplicity, we express gK (p,2, p) as its non-polar V H form NK (p). The di!erence of them is therefore V H f ( p,2, p;1, q,2, q)!f (p,2, p; 0, q,2, q)" H H H H d [gL (1, q,2, q)!gL (0, q,2, q)] NK (p). VW W H W H V H VH\K WH\L (6.4)
X "
6. Partial-derivatives-based GPV system
499
d H H
*gL (q ) *gL (q ) 2 *gL (q ) H\L H H\L H H\L H L> *gL (q ) *gL (q ) 2 *gL (q ) H\L> H H\L> H H\L >H L> . $ $ $ $ *gL (q ) H H
*gL (q ) H H
2
*gL (q ) H H L>
500
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
The number of free parameters in this case is m#n#2. Note that the parameters p and q should be G G chosen such that the partial derivatives to be speci"ed are within the 2D polynomial segment under consideration. Similar to the previous case, we calculate the control points by 1 Q " ) X \P > \ m
(6.7)
This shows that Q can be evaluated from P if X and > are invertible. Based on partial derivatives evaluated for the u direction, we also obtain f (1, p,2, p; q,2, q) H H H " d gL (q,2, q) gK (1, p,2, p) VW W H V H VH\K WH\L f (0, p,2, p; q,2, q) H
H H d gL (q,2, q) gK (0, p,2, p) " VW W H V H VH\K WH\L and
P "
X "
(6.9)
*fK (p , q ) 2 *fK (p ,q ) H K> *fK (p , q ) 2 *fK (p ,q ) K> $ $ $
,
NL (q ) NL (q ) H\L> H 2 H H NL (q ) 2 NL (q ) H\L> H H H $ : $
,
NL (q ) NL (q ) NL (q ) H\L H L> H\L> H L> 2 H H L>
Q "
d H\K H\L d H\K H\L> $
d H\K> H\L d H\K> H\L> $
...
d H\K H
d H\K> H
...
... $
> "
P "3X Q >
(6.13)
*fK (p , q ) *fK (p , q ) P " *fK (p , q ) *fK (p , q )
*fK (p , q ) *fK (p , q ) *fK (p , q ) *fK (p , q )
d d Q " d d
d d d d
d d d d
d H H\L d H H\L> , $ d
H H
*gK (p ) H H
2
*gK (p ) H H K>
*fK (p , q ) *fK (p , q ) *fK (p , q ) *fK (p , q )
*fK (p , q ) *fK (p , q ) , *fK (p , q ) *fK (p , q )
d d , d d
N (q ) N (q ) X " N (q ) N (q )
N (q ) N (q ) N (q ) N (q )
*g (p ) *g (p ) > " *g (p ) *g (p )
*g (p ) *g (p ) *g (p ) *g (p )
N (q ) N (q ) N (q ) N (q ) *g (p ) *g (p ) *g (p ) *g (p )
N (q ) N (q ) , N (q ) N (q )
*g (p ) *g (p ) . *g (p ) *g (p )
We then use the following formula to evaluate the control points: 1 Q " X \P > \ 3
(6.14)
6.1. Higher-order derivatives
*gK (p ) *gK (p ) 2 *gK (p ) H\K H H\K H H\K H K> *gK (p ) *gK (p ) 2 *gK (p ) H\K> H H\K> H H\K> H K> . $ $ $ $ *gK (p ) H H
This shows that Q can be evaluated from P if X and > are invertible. For example, assume n"m"3, j"6, to construct a linear system for tensor product B-spline surface based on partial derivatives in the u direction, we can write the following:
(6.8)
*fK (p , q ) *fK (p , q ) *fK (p ,q ) H L> L> 2 K> L> NL (q ) H\L H NL (q ) H\L H $
(6.12)
(6.11)
where
*fK (p , q ) H *fK (p , q ) H $
1 Q " X \P > \ n
where
H H *fK (p, q)" *gK (p) d NL (q) . (6.10) H V H VW W H VH\K WH\L Consequently, we obtain the following equation: P "nX Q >
In this case, we also have
From [1], the formula for calculating the second derivative of a polynomial curve F at u is stated as follows. F (u)"n ( f (u#1, u#1, u,2, u)!2f (u#1, u, u,2, u) #f (u, u, u,2, u))
(6.15)
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
where n "n(n!1) and in general nIM " n(n!1)2(n!k#1). Similarly, we rewrite Eq. (6.1) as
*F H (p;q)"m f (p,2, p, 1,1, q,2, q ) H GHI GFHFI *v L K !2f (p,2, p, 1,0, q,2, q ) H GHI GFHFI L K #f (p,2, p, 0,0, q,2, q ) . (6.16) H GHI GFHFI L K Each function f in the right-hand side can be expressed H as the following equations:
f (p,2, p; 1,1, q,2, q) H H H d gL (1,1, q,2, q) gK (p,2, p) " VW W H V H VH\K WH\L
where P "
*f(p , q ) H *f(p , q ) H $
*f(p , q ) H *f(p , q ) H $
*f(p ,q ) H K>
*f(p ,q ) H K>
f (p,2, p; 1,0, q,2, q) H H H " d gL (1,0, q,2, q) gK (p,2, p), VW W H V H VH\K WH\L
f (p,2, p; 0,0, q,2, q) H H H " d gL (0,0, q,2, q) gK (p,2, p). VW W H V H VH\K WH\L
(6.19)
Therefore,
$
2 *f(p ,q ) H K> L>
NK (p ) H\K H NK (p ) H\K H $
NK (p ) H\K> H NK (p ) H\K> H $
NK (p ) H\K H K>
NK (p ) H\K> H K>
d H\K H\L
d H\K H\L>
d Q " H\K> H\L $
d H\K> H\L> $
NK (p ) H H NK (p ) H H $
2 2 $
2 NK (p ) H H K>
2
d H\K H
2 d
H\K> H , $
$
d H H\L
d H H\L>
*g L (q ) H\L H *g L (q ) H\L> H $
*g L (q ) H\L H *g L (q ) H\L> H $
2
*g L(q ) H H
*g L(q ) H H
2
and (6.18)
*f(p , q ) H L> *f(p , q ) H L> $
2 2
X "
(6.17)
501
2
d
H H
*g L(q ) H H L>
(6.23)
From [1], the general formula for the kth derivative of a curve is
;gK (p,2, p). (6.20) VH Again, for convenience, we denote f (p,2, p; 1,1, q,2, q)!2f (p,2, p; 1,0, q,2, q)#f (p, H H H 2, p; 0,0, q,2, q) as *f(p, q) and W gL (1,1, q,2, q)! H W H 2gL (1,0, q,2, q)#gL (0,0, q,2, q)X as *gL (q). Hence, W H W H W H Eq. (6.20) can be simpli"ed:
H H *f(p, q)" d ) *gL (q) NK (p). (6.21) H VW W H V H VH\K WH\L Similar to Eq. (6.6), a GPV linear system for a TPBS can be constructed based on second partial derivatives as follows: P "m X Q >
*g L (q ) H\L H L> 2 *g L (q ) H\L >H L> . $ $
1 Q " ) X \P > \. m
!2gL (1,0, q,2, q)#gL (0,0, q,2, q)] W H W H
,
> "
Hence, for calculating the control points, we have
f (p,2, p; 1,1, q,2, q)!2f (p,2, p; 1,0, q,2, q) H H #f (p,2, p;0,0, q,2q) H H H " d [gL (1,1, q,2, q) VW W H VH\K WH\L
,
(6.22)
k FI(u)"nIM (!)I\H j WHWI ;f ( u#1,2, u#1, u,2, u ). (6.24) GFFHFFI GHI L\H H Thus, it is possible for us to construct linear systems for higher-order derivatives in a similar way.
7. Fast computation of basis function matrices In the previous sections, we have derived various formulations for constructing GPV-based linear systems. Usually, in each of these schemes, two matrices consisting of polar forms of B-spline basis function are involved. We call these matrices as basis function matrices. In the
502
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
following, we will present a fast computation algorithm for these basis function matrices. As stated in [6], the basis function gN (u ,2, u ) is computed recursively as GH P follows: g ( )"d G H G H and
(7.1)
u !t G gP\(u ,2, u ) gP (u ,2, u )" P G H P P\ t !t G H G>P G t !u P gP\ (u ,2, u ) # G>P> G> H P\ t !t G>P> G> (7.2) for 1)r)p. Eq. (7.1) states that g ( ) equals 1 when G H i"j and it equals zero otherwise while Eq. (7.2) depicts the recursive calculation process. We need to calculate each entry (i.e. gP (;)) in the basis G H function matrix. An example has been given in Eqs. (4.4) and (4.5). Let us take a look at the row Wg (; ), g (; ), g (; ), g (; )X of the basis function matrix in Eq. (4.5). In order to evaluate the entries of this row according to Eqs. (7.1) and (7.2), we perform recursive calculation for each basis function. This evaluation process is illustrated in Figs. 1a}d. Using this strict forward computation strategy, to evaluate all the gN in a row, the GH number of linear combinations required is evaluated by the following formula:
j"
p#2 p
!(p#1) ;(p#1)"
;(p#1)"[p(p#1)]/2
p#1
Fig. 1. This "gure compares the strict forward computation of a row in the basis function matrix (a)}(d) and a more e$cient algorithm (e).
p!1
(7.3)
where
p#1 p!1
is the number of nodes in the computation structure, except that of the bottom row and (p#1) is the number of entries in a row. Similar form of expression for the number of nodes in a triangular computation structure can be found in [13]. The main idea of our computation algorithm is the following. Instead of performing calculation separately for each entry of in a row, we consider the computation of a single basis function, which subsumes all the required calculations. Explicitly, for a row starting with the entry gN , we perform the recursive evaluation process of gN. GH GH Fig. 1e depicts our computation structure. We consider the whole triangular computation structure as four triangular regions as shown in Fig. 2. Obviously, the calculations involved in the upper triangular region are wasted since we do not really need to compute gN, at the top of GH the triangle. Furthermore, since the computations in the
Fig. 2. The computation structure is divided into three parts and only the central part is utilized in the calculation.
right and left triangular regions involved only zero operands, they are in fact not productive to the computation of the row entries. Fig. 1e highlights the productive computations by bold lines. If the linear combinations in these non-productive regions are eliminated, we have the following number of linear combinations:
u"
2p#1
p#1
!3
2p!1
p!1
(7.4)
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
where
2p#1 2p!1
Table 1 Percentage gain across di!erent values of p
and 3
p#1
p Percentage gain (%)
p!1
represents the total number of linear combinations in the whole triangular computation structure and the nonproductive triangular regions. With some arithmetic, the gain of our algorithm compared with the strict forward computation algorithm is formulated as the following equation: h"j!u"p(p#p#2)/2.
(7.5)
Gain in terms of percentage is calculated as follows: o"h/j"(p#p#2)/(p#1).
503
(7.6)
Table 1 shows the percentage gain across di!erent p values. Our overall algorithm is summarized using the following piece of C-like pseudo code, which is constructed by utilizing the above observations. /* Each U[i] is an array of yoating point numbers representing the arguments in a polar form.*/ Matrix}evaluation(int p, int curve}seg}index, array of U) + int row; for(row "1; row("p#1; row##) Row}evaluation(p, curve}seg}index, U[row]); , /* This function calculates a row of a basis function matrix. */ Row}evaluation(int p, int j, #oat u ,u ,2,u ) + Calculate all the multiplying factors (i.e. the values used for multiplying the basis functions) represented in bold lines in Fig. 1e; /* This requires p and j. */ Starting from g (u , u ,..., u ) which is equal to N N N one, calculate the sub-basis function values in the reverse direction according to Fig. 1e; , The above pseudo code only provides the kernel of the basis function matrix calculation. As there are minor changes in calculating di!erent basis function matrices when di!erent conditions are speci"ed (e.g. degree elevation or partial derivatives), appropriate modi"cations in the above code are required. Besides the above algorithm, we also observe that we can save computation in calculating *gL (q)" GH gL (1, q,2, q)!gL (0, q,2, q) in the partial derivativeGH GH based GPV system. As an example, the following equations show the evaluation of g (1, q, q). By the symmetric
2 88.9
3 87.5
4 88.0
5 88.9
6 89.8
property of polar form arguments, we have 1!t g (q, q) g (1, q, q)"g (q, q,1)" t !t t !1 # g (q, q) t !t where
(7.7)
t !q q!t g (q)# g (q, q)" g (q), t !t t !t q!t t !q g (q)# g (q, q)" g (q) t !t t !t t !q q!t g ( )# g ( ), g (q)" t !t t !t q!t t !q g ( )# g (q)" g ( ), t !t t !t t !q q!t g ( )# g (q)" g ( ). t !t t !t From the above formulas, we observe that when we calculate g (0, q, q), most of the computations have al ready done previously in "nding out the value of g (1, q, q). In fact, we only need to do computation for the last recursive step. After showing our fast basis matrix computation algorithm, in the next section, we will present the experimental results based on the various GPV-based schemes derived in the previous sections.
8. Experiments In this section, we will present some examples which the aforementioned GPV-based schemes are employed as control point speci"cation mechanisms. In all cases, we need to calculate a set of general polar values. Strictly speaking, we need to "nd out the piecewise polynomial (p.p.) forms of the underlying B-spline curves. However, it is quite a computationally complicated task. Rather, just to ful"ll the purpose of illustrating the techniques in practice, we employ a simpler way in performing our experiments. That is, we set an arbitrary value to the "rst entry of the matrix P. We then "ll up the other entries by multiplying it with an appropriate factor obtained from the given polar forms. For example, if the "rst entry
504
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
f (u , u , u ; v , v , v ) is set to a value C, then an entry f (u , au , u ; v , v , bv ) will set to the value a;b;C. The shape of the resulting control points will only di!er in scale with the ones, which are computed by actually calculating the exact values. In the following, we will give the parameters used for constructing the examples. In all the cases, the degree of the surface is set to 3 in both directions and the knot vector is set to (0, 0, 0, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 12, 12, 12) with 15 control points along each direction. The control points we want to specify lay on the central part, as illustrated in Fig. 3. In the following cases, j is set to 8. For the direct extension and the degree elevation approaches, we set ; "(3, 6, 5), ; "(3.9, 6, 5), ; "(3, 6, 7), ; "(3, 12, 5), < "(1, 2, 3), < "(1.3, 2, 3), < "(1, 2, 4.2), < "(1, 4, 3). For the direct extension scheme, based on the formulation shown in Eq. (4.14), we obtain
0.50 0.65 0.70 1.00
P"
0.65 0.85 0.91 1.30 0.70 0.91 0.98 1.40
,
Fig. 3. This "gure shows our 15;15 control point mesh structure. The central part is highlighted with circles and is speci"ed through general polar values schemes. In case of degree elevation approach, only the white circle control points are speci"ed while in the other two cases, all the circled control points are speci"ed.
1.00 1.30 1.40 2.00
X"
0.00 1.33 !0.33
0.00
0.00 1.03 !0.03
0.00
0.00 0.00
1.67 !0.67
!3.00 5.33 !1.33
!4.00
!3.70
15.00
13.9
!20.00
!18.6
10.00
9.4
>"
0.00
>"
,
5.33
!10.00 !8.00 6.00
5.00
Note that, in this case, we set 0.5 to the "rst entry in P. The value 0.65 in the next entry is calculated by multiplying it with 1.3 which can be observed from the di!erence between ; and ; . For the degree elevation scheme, based on Eq. (5.6), we have
0.50 0.65 0.70 1.00
P" 0.65 0.85 0.91 1.30 , 0.70 0.91 0.98 1.40
XM "
0.50 0.83 !0.33
0.35 0.83 !0.18 ,
!0.50 2.17 !0.67
!3.70
15.00
13.90
!1.60 !1.33 6.60
5.33
!20.00 !18.60 !10.00 !8.00
!1.60 !1.33 6.60
!4.00
10.00
9.4
6.00
5.00
.
In this case, we assume we raise the degree from 2 to 3 in u direction of the TPBS. For the partial-derivative-based scheme, we set p "5.50, p "5.78, p "6.05, p "6.325, q "5.00, q "5.25, q "5.50, q "5.75. Based on Eq. (6.6), we have
0.20 0.22 0.24 0.26
P "
0.23 0.26 0.28 0.31 0.27 0.29 0.32 0.35
,
0.30 0.34 0.37 0.40
0.02 0.48 0.48 0.02
X "
0.00 0.30 0.62 0.08 0.00 0.14 0.66 0.19
,
0.00 0.07 0.54 0.39
> "
0.00
0.01
0.04
0.09
0.17
0.22
0.21
0.14
0.00 !0.13 !0.21 !0.22
!0.17 !0.09 !0.04 !0.01
.
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
505
A dome-shaped control point mesh is constructed for the purpose of comparison (Fig. 4a). The corresponding TPBS is shown in Fig. 5a while the corresponding control point values (height) of the central part (i.e. circled parts of Fig. 3) are depicted in Fig. 6a. The results of the three GPV-based control point speci"cation schemes with parameter settings given above are shown in Figs. 6b}d. To translate the calculated control points to the top for the sake of better visualization, an o!set of 550 has been added to the calculated values. The control meshes with new control points at the central part and the corresponding TPBS are shown in Figs. 4 and 5.
9. Discussion In this paper, we have extended the GPV-based control point speci"cation method in [16] to tensor product B-spline surfaces and proposed two more alternatives in the speci"cation technique, which are based on degree elevation and partial derivative evaluation, respectively. In the previous sections, we focus on the evaluation of an individual region of a TPBS. In fact, we can potentially extend our formulation to control the whole surface. Here, we give the formulation in the curve level. For example, given a normalized cubic B-spline with knot vector ¹"(0, 0, 0, 0, 1, 2, 4, 5, 6, 6, 6, 6) (this con"guration of the curve is taken from [5]), based on (4.1) we write the following formula: P"AQ
(9.1)
where P"[ f (; ), f (; ), f (; ), f (; ), f (; ), f (; ), f (; ) f (; )]2 Q"[d , d , d , d , d , d , d , d ]2,
g (; ) g (; ) g (; ) g (; ) g (; )
A"
g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; ) g (; ) g (; )
The empty slots in A are "lled with zeros and U is G a suite of arguments. Other con"gurations set by changing the corresponding entries of P and A of the above formulations are of course possible. The control points Q can be obtained by Q"A\P if A is invertible. We can
Fig. 4. The shapes of various control point meshes produced using di!erent kinds of GPV-based scheme: (a) Mesh for comparison, (b) direct extension scheme, (c) degree elevation and (d) partial-derivatives-based schemes.
then follow the technique mentioned in Section 4.2 to extend this formulation to a TPBS. As shown in the experiment section, the drawback of this kind of control point speci"cation scheme is the following. Since the ultimate shape and location of the resulting control points is not easy to predict initially, the user may need to try out many di!erent parameters in the process of setting up a desirable control point mesh. Alternately, as illustrated in the previous section, the user may need to add o!sets to the resulting points to make it
g (; ) g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; ) g (; ) g (; ) g (; )
look "t to the surface as a whole. The key to this problem is a technique for restricting the range of the values of the evaluated control points. This may require the overall shape information of the surface and that will be our next research work.
506
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
Fig. 5. The shapes of the corresponding tensor product surfaces generated by the meshes shown in Fig. 4: (a) Surface for comparison, (b) direct extension scheme, (c) degree elevation and (d) partial-derivatives-based schemes.
Fig. 6. The control point values produced by our schemes: (a) Control points for comparison, (b) direct extension, (c) degree elevation and (d) partial-derivatives-based schemes.
10. Conclusions Recently, a GPV-based control point speci"cation scheme is proposed in [16]. In this paper, we "rst construct linear systems of tensor product B-spline surfaces (TPBS) under the GPV-based framework. We name this
extension as the direct extension scheme. Subsequently, based on similar principles, we investigate the issue of extending the linear systems based on degree elevation and partial derivatives evaluation. As inherited from the original scheme [16], a large number of free parameters are provided in the direct extension and degree elevation schemes. In the partial derivative schemes, partial derivatives on a surface can be used as parameters, which is a desirable feature as slope information is important in controlling the shape of the surface. In all cases, we need to choose appropriate parameters by conforming to the completeness criteria. In other words, the related matrices should be invertible. Furthermore, in case of partial derivative scheme, the parameters should be chosen such that the partial derivatives to be evaluated are within the 2D polynomial segment under consideration. Finally, we introduce an e$cient way of computing the basis function matrices. We discover that the calculations of the entries of a row of a basis function matrix can be reduced if we combine their calculations together. Large computation reduction is found and has
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
been summarized in Table 1. Finally, we further observe that partial results can be reused in the case of partial derivative scheme. Besides the research work suggested in the previous section, the idea of using general polar values as control points for surfaces raises other important issues for future research. They are the extensions of these formulations to other B-spline surface types, for example, triangular B-spline surfaces [9,20] and geometric continuous Bspline (i.e. b-spline) tensor product surfaces [21]. Triangular B-spline surfaces can model complex objects with non-rectangular topology while b-spline provides extra shape control parameters based on geometric continuity. Light has already thrown on how to construct those surfaces using blossoming [8,10] and we plan to extend our research work to those types of surfaces.
Acknowledgements The author thanks reviewer's valuable comments and Prof. Stephen Mann of University of Waterloo for his help in providing the technical report [22] through the Internet.
References [1] Ramshaw L. Bezier and B-splines as multia$ne maps. In: Earnshaw RA, editor. Theoretical foundations of computer graphics and CAD. Berlin: Springer, 1988. p. 757}76. [2] DeRose T, Goldman R. A tutorial introduction to blossoming. In: Hagen H, Roller D, editors. Geometric modeling methods and applications. Berlin: Springer, 1991. p. 267}86. [3] Seidel HP. An Introduction to polar forms. IEEE Computer Graphics and Applications 1993;13(1):38}46. [4] Ramshaw L. Blossoms are polar forms. Computer Aided Geometric Design 1989;6:323}58. [5] Seidel HP. Knot Insertion from a blossoming point of view. Computer Aided Geometric Design 1988;5:81}6. [6] Seidel HP. A new multia$ne approach to B-splines. Computer Aided Geometric Design 1989;6:23}32.
507
[7] Seidel HP. Computing B-spline control points. In: Straber W, Seidel HP, editors. Theory and practice of geometric modeling. Berlin: Springer, 1986. p. 17}32. [8] Dahmen W, Micchelli CA, Seidel HP. Blossoming begets B-spline bases built better by B-patches. Mathematics of Computation 1992;59(199):97}115. [9] Fong P, Seidel HP. An implementation of triangular B-spline surfaces over arbitrary triangulations. Computer Aided Geometric Design 1993;10:267}75. [10] Seidel HP. Polar forms for geometrically continuous spline curves of arbitrary degree. ACM Transactions on Graphics 1993;12(1):1}34. [11] Lee ETY. Computing a chain of blossoms, with application to products of splines. Computer Aided Geometric Design 1994;11(6):597}620. [12] DeRose TD, Goldman RN, Hagen H, Mann S. Functional composition algorithms via blossoming. ACM Transactions on Graphics 1993;12(2):113}35. [13] Liu W, Mann S. An optimal algorithm for expanding the composition of polynomials. ACM Transactions on Graphics 1997;16(2):155}78. [14] Goldman R. Blossoming with cancellation. Computer Aided Geometric Design 1999;16(7):671}89. [15] Goldman R. The rational Bernstein bases and the multirational blossoms. Computer Aided Geometric Design 1999;16(8):701}38. [16] Duchaineau MA. Using general polar values as control points for polynomial curves. Computer Aided Geometric Design 1994;11:411}23. [17] Eck M, Handenfeld J. Knot removal for B-spline curves. Computer Aided Geometric Design 1995;12:259}82. [18] Farin G. Curves and surfaces for computer aided geometric design: a practical guide, 4th ed. NewYork: Academic Press, 1998. [19] Mann S, DeRose T. Computing values and derivatives of bezier and B-spline tensor products. Computer Aided Geometric Design 1995;12:107}10. [20] Seidel HP. Triangular B-splines. In: Strasser W, Wahl F, editors. Graphics and robotics. Berlin: Springer, 1995. p. 149}62. [21] Barsky BA. In: Kunii TL, editor. Computer graphics and geometric modeling using beta-splines. In: Berlin: Springer, 1988. [22] Mann S, DeRose T, Winkenback G. Computing values and derivatives of bezier and B-spline tensor products. Research Report CS-93-31, Computer Science Department, University of Waterloo, 1993.
Technical Section
Extensions of the general polar value based control point speci"cation method in constructing tensor product B-spline surfaces Stephen Wang-Cheung Lam* Department of Computing, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong
Abstract Blossoming has been utilized as a powerful tool for understanding Bezier and B-spline curves. Recently, general polar values (GPV) have been employed as an aid to specify control points for polynomial curves. Based on this technique, it is possible for us to use a set of polar values to control a polynomial curve. In this paper, we "rst construct linear systems of tensor product B-spline surfaces (TPBS) under the GPV-based framework. We call this extension as the direct extension scheme. Subsequently, based on similar principles, we investigate the issue of extending the linear systems based on degree elevation and partial derivatives evaluation. For the direct extension and degree elevation schemes, high #exibility is provided though a large set of corresponding free parameters. For the partial derivative scheme, partial derivatives on a surface can be used as parameters. In all schemes, a restriction on parameter selection is imposed by the invertibility of the related matrices. Finally, based on the subtle relationships among the basis functions, a fast computation algorithm is developed for evaluating those related basis-function matrices. 2000 Elsevier Science Ltd. All rights reserved. Keywords: General polar value; Control point; Tensor product surface; Partial derivative; Degree elevation; Basis-function matrices
1. Introduction It is well known that the polar form of a degree n polynomial G(u) is a symmetric, multia$ne polynomial g(u ,2, u ) that satis"es the identity L g( t,..., t )"G(t). GHI L In [1], this principle is applied to the formulation of Bezier and B-spline curves and surfaces that are frequently used in computer-aided geometric design. The word `blossominga is used to name the process of replacement of a polynomial of high degree with its polar form. The main bene"t of that replacement is the dis-
* Tel.: #851-2766-7301; fax: #851-2774-0842. E-mail address: [email protected] (S. WangCheung Lam).
covery of a way to label the Bezier points that control a curve segment or surface patch and the de Boor points that control a B-spline curve, with symmetric, multivariate labels. This labeling gives us great geometric insight into the basic algorithms in this area, such as the de Casteljau algorithm and the de Boor algorithm. A good introduction to this topic can be found in [2,3] and an in-depth treatment of the theory behind blossoming is given in [4]. Blossoming has been applied in developing geometric algorithms and proo"ng existing algorithms with new insights. For example, in [5], the theory of blossoming has been applied to the insertion of a new knot into the knot vector of a B-spline curve. There, a formula that expresses the de Boor points of the re"ned curve in terms of a set of polar values is derived. As a corollary, a new proof of Boehm's knot insertion algorithm has also been obtained. Furthermore, the relationship between multiple knot insertion and the de Boor algorithm is also revealed through blossoming. To extend blossoming to
0097-8493/00/$ - see front matter 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 7 - 8 4 9 3 ( 0 0 ) 0 0 0 5 2 - 2
494
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
B-spline, in [6,7], Seidel attempts to apply blossoming directly to the recurrence relation of B-splines. There, he gives a new derivation for expressing the de Boor points of a B-spline curve in terms of blossoming. Based on that, the Curry}Schoenberg theorem and the de Boor algorithm were re-derived. Although all major theorems in that article are already known, the proofs are shorter than the previous ones and they also give further insights into the basic theory of B-splines construction. Recently, blossoming has been applied to the formulation of triangular B-splines surfaces [8,9] and geometric continuous spline curves [10], computation of product of splines [11] and composition of Bezier simplexes [12,13]. Besides, extensions of the technique to multirational framework can be found in [14,15]. In [16], Duchaineau shows that general polar values (GPV) may also be used to control a polynomial curve when a related matrix is invertible. In that case, those GPV act as blossom control points and the invertible matrix provides a useful translation from these blossom control points to well known ones such as those of Bezier. In [16], the special case in which the polar values can be evaluated through pairwise a$ne combination is also characterized. That characterization provides users of this approach a way of generating GPV. As stated in [16], this GPV approach allows the arguments of the blossom control points to be chosen in a manner akin to choosing the knot vector of a B-spline segment. Furthermore, the main advantage of this approach over the usual B-spline control points is the increased number of free parameters for specifying the arguments of the blossom control points. As stated in [16], one of the important future research works is on the extension of the technique to the area of surfaces. We agree with him and hence in this paper, we extend the GPV approach to surfaces. Besides, we also investigate other possible ways of applying GPV in specifying control points. To ful"ll our goal, we "rst construct the linear system of a tensor product B-spline surface (TPBS) under the GPV-based framework. Subsequently, we extend the system based on degree elevation and partial derivatives calculation. Finally, after a close examination of the calculation of the related basis function matrices, we discover that it is possible to greatly reduce the computation time for evaluating those matrices through using a smarter computing algorithm. The remainder of this paper is structured as follows. In Section 2, we present a brief review of the blossoming of Bezier and B-spline curves. In Section 3, we describe the fundamental concepts and techniques used in general polar values (GPV) based control point speci"cation. In Section 4, we "rst apply GPV in controlling B-spline curve. Then, we extend it to tensor product B-spline surfaces (TPBS). In particular, we focus on controlling an individual region on a surface. We call this extension as direct extension scheme. Moreover, we also investigate
other usage of the related basis function matrix. In Section 5, we modify the GPV system based on degree elevation. Similar to the original scheme [16], the direct extension and degree elevation schemes provide high #exibility by allowing a large set of free parameters. In Section 6, we extend our system through using partial derivative formulation. Furthermore, we study the issue of applying higher-order derivative calculation in controlling a TPBS. In all schemes, a restriction of parameter selection is imposed by the invertibility of the related matrices. In Section 7, a fast computation algorithm is developed for calculating the basis function matrices used in our system. In Section 8, we perform experiments to illustrate the schemes we proposed in controlling a TPBS. In Section 9, we discuss the pros and cons of our approach. Finally, in Section 10, we summarize our contributions and suggest some topics for future research works.
2. Review on blossoming This section aims at brie#y reviewing some of the results in blossoming theory of Bezier and B-spline curves [5]. 2.1. Basis De5nition. A map f : 1P1:(u ,2, u )Pf (u ,2, u ) is L L n-a$ne, if it is a$ne in each argument. In other words, the following formula f (u ,2 a t ,2, u )" a f (u ,2, t ,2, u ) (2.1) G G L G G L G G holds for all a , t 31 with a "1. Furthermore, f is G G G G symmetric, if it keeps itsvalue under any permutation of its arguments. Blossoming Principle: A one-to-one correspondence exists between polynomial functions F : 1P1B of degree n and symmetric n-a$ne maps f : 1P1B. Further, given a function of either type, a unique function of the other type exists that satis"es the identity F(u)"f (u,2,u), where f is called the blossom of F and F is the diagonal of f. 2.2. Bezier and B-spline control points If F is taken as a Bezier curve over [s t], its control points b ,2, b are given by the formula L b "f ( s,..., s, t,..., t ). (2.2) G GHI GHI L\G G Obviously, for spline curves the situation becomes more complicated. The following formula represents
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
a non-periodic B-spline curve K F(t)" d ) N (t) G GI G
where (2.3)
of order k (i.e. degree n"k!1) over the non-decreasing knot sequence ¹"(t "t "2"t , t ,2t , t "2"t ) I\ I K K> K>I that does not contain any intermediate knots t with G t )t )t of multiplicity greater than k!1. In Eq. I G K (2.3), d ,.., d are the de Boor points (B-spline control K points) and N are the normalized B-splines of order GI k over T. For t (t , we denote by F the polynomial of the H H> H B-spline curve F in the segment (t , t ). F is of order H H> H k with blossom f . For every i)j)i#k!1 with H t (t the de Boor point d is related to the blossom H H> G f by the formula H d "f (t ,2, t ) G H G> G>I\
(2.4)
Furthermore, since every knot t appears with multipliG city at most k, such a j always exists. After a brief review of the blossoming theory, we will introduce the GPV-based control point speci"cation scheme in the next section. Readers who are interested in other fundamental theorems and results obtained in blossoming theory can refer to [1,4,5].
3. General polar values In [16], Duchaineau shows that general polar values (GPV) may also be used to control polynomial curves when a related matrix is invertible. In this section, we brie#y review this methodology. Any nth-degree polynomial curve F(x) can be written as L F(x)" Q B (x) G G G
(3.1)
where Q ,2, Q are control points and B (x),2, B (x) L L are any known basis for the nth-degree polynomials. One example is the Bezier basis. The polar form f of F can be written as L f (x ,2, x )" Q b (x ,2, x ) L G G L G
(3.2)
where b ,2, b are the polar forms of the basis functions. L Obviously, the set of polar values P can be obtained by the application of a linear transformation A to the Q control points, that is P"AQ
495
(3.3)
P f (t ,2, t ) L P f (t ,2, t ) L P" " $ $ P L
and
A"
(3.4)
f (t ,2, t ) L LL
b (t ,2, t ) 2 b (t ,2, t ) L L L b (t ,2, t ) 2 b (t ,2, t ) L L L . $ b (t ,2, t ) L LL
2
b (t ,2, t ) L L LL
(3.5)
It is interesting to know whether a set of n#1 polar values with general argument suite is capable of controling the nth-degree polynomial curves. In [16], the word complete is used to describe the polar values [ f (t ,2, t ),2, f (t ,2, t )] when a polynomial can L L LL be produced by setting P appropriately. This occurs precisely when the matrix A is invertible, since this allows Q to be obtained from P using Q"A\P. That implies there is a one-to-one correspondence between P and Q. The following provides the de"nition of completeness. De5nition 2.1 (Duchainean [16]). A set of n#1 polar values [f (t ,2, t ),2, f (t ,2, t )] is complete if and L L LL only if the linear transformation to these polar values form control points in any known basis of the nth-degree polynomials has a nonzero determinant. If the polar values have been proved to be complete, they will be useful as control points and shall therefore be called blossom control points [16].
4. GPV based B-spline system From [6], F (u)" H d ) NN (u) represents the H GH\N G G H polynomial in the non-empty interval I :"[t , t ) in H H H> a B-spline curve F(u), where t (t . The corresponding H H> polar form of F (u) is expressed as follows: H H f (u ,2, u )" d ) gN (u ,2, u ). (4.1) H N G G H N GH\N Based on this formula, we construct the linear system for B-spline curves. Following the notation in [16], we de"ne the GPV linear transformation of a B-spline curve as follows: P"AQ.
(4.2)
The entries of P is f (;) where and f is the polar form H H of the polynomial starting from the knot t and U is an H argument suite. Q represents the control points. An entry
496
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
in A is the polar form of a B-spline basis function. We denote it as gN (;) where p is the degree of the spline. GH In general, we have
f (; ) H f (; ) H $
gN (; ) H\N H gN (; ) H\N H $
gN (; ) H\N> H gN (; ) H\N> H $
gN (; ) H\N H N>
gN (; ) H\N> H N>
... ... ...
gN (; ) H H gN (; ) H H $
... gN (; ) H H N>
(4.3)
d H
For example, given a normalized cubic B-spline curve with knot vector ¹"(0, 0, 0, 0, 1, 2, 4, 5, 6, 6, 6, 6) and control point vector Q"(d , d , d , d , d , d , d , d ) [5], when j"5 and j"6, according to (4.3) we have, respectively, f (; ) g (; ) f (; ) g (; ) " f (; ) g (; ) f (; ) g (; )
;
and
d
d
d
d
f (; ) g (; ) f (; ) g (; ) " f (; ) g (; ) f (; ) g (; )
;
De5nition 4.1. A set of n#1 polar values + f (; ),2, f (; ),, where f is the polar form of the H H > H polynomial between t and t and n is the number of H H> control points, is complete if the linear transformation to these polar values form the corresponding set of control points of a normalized B-spline curve has a nonzero determinant. 4.1. Other usage of the transform
d H\N> . $
(4.6)
Hence, similar to [16], here completeness is further de"ned as follows:
d H\N
;
If A is invertible, we can obtain the control points Q of a B-spline curve segment through the following equation: Q"A\P
f (; ) H N>
"
where ; is a suite of arguments. Notice that ; in Eqs. G G (4.4) and (4.5) may be di!erent.
d
d
d
d
g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; )
(4.4)
Similar to [16], the matrix A given above turns out to be useful for more than just testing for the completeness of the set of polar values. With Bezier basis, matrix A is shown to be useful for Bezier curve composition based on derivatives [16]. Here, we suggest that matrix A can also be used in de"ning new GPV-based blossom control points in case of knot insertion. In particular, after knot insertion, new control points Q* are obtained from the old ones Q through a linear transformation (p. 262 of [17]): QH"BQ
(4.7)
Based on (4.6), the new GPV vector can therefore be expressed as P*"(A*BA\)P.
(4.8)
This equation is now in the same form of Eq. (4.7) and, similar to the ordinary control points, we can then obtain the new set of blossom control points from the old ones. 4.2. Tensor product surface
g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; )
(4.5)
We "rst review the de"nition and some properties of tensor product surfaces. Most of the materials in this review are taken from [1]. Let F : 1;1PO be a polynomial function; that is, each coordinate of the point F(u, v) is given by a polynomial in the two real valuables u and v. In a (m; n) degree tensor product surface, the coordinates of F(u, v) are required to have degree no higher than m in the variable u and degree no higher than n in the variable v separately. In a triangular patch surface of degree n, the coordinates of F(u, v) are required to have total degree no higher than n in the variables u and v jointly. In this paper, we will not consider triangular patch surfaces.
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
The blossom f of a tensor product surface F of degree (m; n) is an (m#n) a$ne function f : 1 ;1PO that is symmetric in its "rst m arguments and its last n arguments. f also satis"es the identity F(u, v)" f (u,2, u;v,2, v). Also, every tensor product surface has a unique blossom. If [p, q];[r, s] is a reference rectangle for the domain plane 1;1, the (m#1)(n#1) Bezier points of the surface F are the blossom value f ( p,2, p, q,2, q, r,2, r, s,2, s ) (4.9) GHI GHI GHI GHI K\G L\H G H for i in [0, m] and j in [0, n]. The de Casteljau algorithm for a tensor product surface calculates the blossom value f (u ,2, u ; v ,2, v ) by performing m#n stages of lin K L ear interpolations, with each blossom argument providing the ratio for one of these stages. Moreover, these stages can be performed in any order. As stated in [18], a tensor product surface F(u, v) can be evaluated as follows: K L F(u, v)" d NK(u)NL (v). VW V W V W The above formula can be rewritten as
(4.10)
L K (4.11) F(u, v)" d NL (v) NK(u). VW W V V W Based on (4.1), we derive the polar form of the above formula as the following equation: f (u ,2, u ; v ,2, v ) GH K L H G d gL (v ,2, v ) gK (u ,2, u ). " VW W H L V G K VG\K WH\L (4.12)
For simplicity, we just concentrate on the case where i"j in the following parts of this paper and (4.12) is therefore rewritten as f (u ,2, u ; v ,2, v ) H K L H H " d gL (v ,2, v ) gK (u ,2, u ). VW W H L V H K VH\K WH\L (4.13)
To construct the linear system for a tensor product B-spline surface, similar to the case of B-spline curve, we write the following: P"XQ>
(4.14)
where P, X, Q and > are given by
P"
f (; , < ) H f (; , < ) H $
f (; , < ) H f (; , < ) H $
2 2 2
f (; , < ) H L> f (; , < ) H L> $
f (; , < ) f (; , < ) 2 f (; ,< ) H K> H K> H K> L>
,
X"
gK (; ) H\K H gK (; ) H\K H $
497
gK (; ) H\K> H gK (; ) H\K> H $
gK (; ) H H gK (; ) H H $
2 2 2
gK (; ) gK (; ) 2 gK (; ) H\K H K> H\K> H K> H H K>
,
d d 2 d H\K H\L H\K H\L> H\K H d d 2 d H\K> H , Q" H\K> H\L H\K> H\L> $ $ 2 $ d H H\L
d H H\L>
d H H
2
gL (< ) gL (< ) 2 gL (< ) H\L H H\L H H\L H L> gL (< ) gL (< ) 2 gL (< ) H\L> H L> . >" H\L> H H\L> H $ $ 2 $ gL (< ) H H
gL (< ) H H
gL (< ) H H L>
2
In our formulation, P and Q may not be square matrices while X and > are constructed as square matrices. In this formulation, the number of free parameters (i.e. the number of di!erent arguments in the suites) is m(m#1)#n(n#1). With some matrix arithmetic, we have the following formula Q"X\P>\.
(4.15)
This shows that if X and > are invertible, Q can be calculated from P. Suppose j"3, n"2 and m"3, P, X, Q and > are given by
f (; , < ) f (; , < ) P" f (; , < ) f (; , < )
f (; , < ) f (; , < ) f (; , < ) f (; , < )
d d Q" d d
d
d
,
d d d d
g (; ) g (; ) X" g (; ) g (; ) and
g (< ) >" g (< ) g (< )
d d
f (; , < ) f (; , < ) , f (; , < ) f (; , < )
g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; )
g (< ) g (< ) g (< )
g (< ) g (< ) . g (< )
498
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
In order to distinguish this technique with the later ones, we herein name it as direct extension scheme.
In this section, we investigate the issue of using degree elevation formulation for constructing a GPV based control point speci"cation system. Let F(u) be a polynomial curve of degree two and f (u, v) its polar form. If we consider F(u) to be a degenerate cubic polynomial curve, the blossom of F(u) is a tria$ne function, say g(u, v, w). Then, we have F(u)"f (u, u)"g(u, u, u). In [1], Ramshaw has stated the following formula: g(u, v, w)"[f (u, v)#f (v, w)#f (w, u)]/3.
(5.1)
We can understand it through revealing that righthand side of the above equation is a symmetric, tria$ne function whose diagonal agrees with F. As stated in [1], for arbitrary degree n and an arbitrary parameter space P, we have the following equation. 1 L> g(u ,2, u )" f (u ,2u( ,2, u ). (5.2) L> L n#1 G The hat over the argument u indicates that the ith G argument is omitted. Similar to the previous section, we can express a B-spline curve based on (5.2) as 1 L P"AM Q where AM " A ? n ? where P, A and Q are given by G
gL\ (;K ) gL\ (;K ) 2 gL\(;K ) H\L> H H\L> H H H gL\ (;K ) gL\ (;K ) 2 gL\(;K ) H\L> H H H A " H\L> H ? $ $ 2 $ 2
gL\(;K ) H H L
d H\L> d Q" H\L> $ d
;gK\(u ,.u( ., u ) V H J K
H
where ; denotes u ,2, u and ;K in A denotes G L G ? u ,2u ,2, u where u is omitted. Similarly, to ? L ?
(5.5)
To construct a degree-elevation-based GPV system, using distributive law, we derive the following: 1 K P"XM Q> where XM " X (5.6) ? m ? and where P, X , Q and > are given by G f (; , < ) f (; , < ) 2 f (; , < ) H H H L> f (; , < ) f (; , < ) 2 f (; , < ) H H L> , P" H $ $ 2 $
f (; , < ) H K
gK\ (;K ) H\K> H gK\ (;K ) X " H\K> H ? $
d H H\L
f (; , < ) 2 f (; , < ) H K H K L> gK\ (;K ) H\K> H gK\ (;K ) H\K> H $
gK\ (;K ) gK\ (;K ) H\K> H K H\K> H K
g (; ) H L
and
f (u ,2, u ;v ,2, v ) H K L 1 K H H " d gL (v ,.., v ) VW W H L m J VH\K> WH\L
d H\K> H\L d Q" H\K> H\L $
gL\ (;K ) H\L> H L
(5.4)
Based on Eqs. (4.13) and (5.2), if we raise the degree of a TPBS in the u direction from m!1 to m, we have
(5.3)
g (; ) H g (; ) P" H , $
gL\ (;K ) H\L> H L
Q"AM \P.
5.1. Using degree elevation in GPV based TPBS system
5. Degree-elevation-based GPV system
calculate the control points, we have
d
H\K> H\L>
d
H\K> H\L> $ d H H\L>
2 2 2
gK\(;K ) H H gK\(;K ) H H , $
2 gK\(;K ) H H K
2 d H\K> H 2 d H\K> H , 2 $ 2
d H H
gL (< ) gL (< ) 2 gL (< ) H\L H H\L H H\L H L> gL (< ) gL (< ) 2 gL (< ) H\L> H L> >" H\L> H H\L> H $ $ 2 $ gL (< ) H H
gL (< ) H H
2
gL (< ) H H L>
where ; denotes u ,2, u and ;K in X denotes G K G a u ,2, u ,2, u where u is omitted. In this case, the ? K ? number of free parameters are m#n(n#1). To evaluate the control points, we have Q"XM \P>\.
(5.7)
A similar set of equations can be derived if we choose to elevate the degree of the other direction. Now, if we want to raise the degree in both directions, we can
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
construct the following formula:
Based on Eq. (4.13), we obtain the following two formulas:
f (u ,2, u ; v ,2, v ) H K L 1 L 1 K H H " d gL\ VW W H n m P R VH\K> WH\L>
;(v ,.,v( ,., v ) gK\(u ,.u( ., u ) P L V H J K
(5.8)
.
To construct a linear system based on the above, we write the following: 1 K P"XM Q>M with XM " X ? m ?
1 L and >M " > @ n @ (5.9)
and where P and > are given by @
f (; , < ) H f (; , < ) P" H $ f (; , < ) H K
and
f (; , < ) H f (; , < ) H $
2 2
f (; , < ) H L f (; , < ) H L $
f (; , < ) 2 f (; , < ) H K H K L
gL\ (
2
(5.10)
gL\ (
2
gL\(
< denotes v ,2, v and
(5.12)
f ( p,2, p; 1, q,2, q) H H H " d gL (1, q,2, q) gK ( p,2p) V H VW W H VH\K WH\L f ( p,2, p; 0, q,2, q) H
In this section, we will explore the construction of partial-derivatives-based GPV system. An extension of this system to higher-order derivatives is presented in Section 6.1. Let F(u; v) be a degree n;m tensor product surface. From [19], we obtain the following partial evaluation formula: *F (p; q) *v
For convenience, we denote f (p,2, p; 1, q,2, q)! H f (p,2, p; 0, q,2, q) as *f (p, q) and [gL (1, q,2, q)! H H W H gL (0, q,2, q)] as *gL (q). Hence, Eq. (6.4) can be simpliW H W H "ed as
H H *f (p, q)" d ) *gL (q) NK (p). (6.5) H VW W H V H VH\K WH\L To construct a GPV linear system for a TPBS using partial derivatives, we write the following: P "mX Q >
(6.6)
where
P "
*f (p , q ) H *f (p , q ) H $
*f (p , q ) H *f (p , q ) H $
(6.1)
*f (p , q ) H L> *f (p , q ) H L> $
2 2 $
,
*f (p , q ) *f (p , q ) 2 *f (p ,q ) H K> H K> H K> L> NK (p ) H\K H NK (p ) H\K H $
NK (p ) H\K> H NK (p ) H\K> H $
NK (p ) H H NK (p ) H H $
2 2 $
,
NK (p ) NK (p ) 2 NK (p ) H\K H K> H\K> H K> H H K> d H\K H\L
d H\K H\L>
2
d H H\L>
2
d
H\K H d d 2 d H\K> H , Q " H\K> H\L H\K> H\L> $ $ $ $ d
H H\L
> "
"m f ( p,..., p, 1, q,..., q ) GHI GHI L K !f ( p,..., p, 0, q,..., q ) GHI GHI L K
(6.2)
H H " d gL (0, q,2, q) gK ( p,2, p) (6.3) VW W H V H VH\K WH\L For simplicity, we express gK (p,2, p) as its non-polar V H form NK (p). The di!erence of them is therefore V H f ( p,2, p;1, q,2, q)!f (p,2, p; 0, q,2, q)" H H H H d [gL (1, q,2, q)!gL (0, q,2, q)] NK (p). VW W H W H V H VH\K WH\L (6.4)
X "
6. Partial-derivatives-based GPV system
499
d H H
*gL (q ) *gL (q ) 2 *gL (q ) H\L H H\L H H\L H L> *gL (q ) *gL (q ) 2 *gL (q ) H\L> H H\L> H H\L >H L> . $ $ $ $ *gL (q ) H H
*gL (q ) H H
2
*gL (q ) H H L>
500
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
The number of free parameters in this case is m#n#2. Note that the parameters p and q should be G G chosen such that the partial derivatives to be speci"ed are within the 2D polynomial segment under consideration. Similar to the previous case, we calculate the control points by 1 Q " ) X \P > \ m
(6.7)
This shows that Q can be evaluated from P if X and > are invertible. Based on partial derivatives evaluated for the u direction, we also obtain f (1, p,2, p; q,2, q) H H H " d gL (q,2, q) gK (1, p,2, p) VW W H V H VH\K WH\L f (0, p,2, p; q,2, q) H
H H d gL (q,2, q) gK (0, p,2, p) " VW W H V H VH\K WH\L and
P "
X "
(6.9)
*fK (p , q ) 2 *fK (p ,q ) H K> *fK (p , q ) 2 *fK (p ,q ) K> $ $ $
,
NL (q ) NL (q ) H\L> H 2 H H NL (q ) 2 NL (q ) H\L> H H H $ : $
,
NL (q ) NL (q ) NL (q ) H\L H L> H\L> H L> 2 H H L>
Q "
d H\K H\L d H\K H\L> $
d H\K> H\L d H\K> H\L> $
...
d H\K H
d H\K> H
...
... $
> "
P "3X Q >
(6.13)
*fK (p , q ) *fK (p , q ) P " *fK (p , q ) *fK (p , q )
*fK (p , q ) *fK (p , q ) *fK (p , q ) *fK (p , q )
d d Q " d d
d d d d
d d d d
d H H\L d H H\L> , $ d
H H
*gK (p ) H H
2
*gK (p ) H H K>
*fK (p , q ) *fK (p , q ) *fK (p , q ) *fK (p , q )
*fK (p , q ) *fK (p , q ) , *fK (p , q ) *fK (p , q )
d d , d d
N (q ) N (q ) X " N (q ) N (q )
N (q ) N (q ) N (q ) N (q )
*g (p ) *g (p ) > " *g (p ) *g (p )
*g (p ) *g (p ) *g (p ) *g (p )
N (q ) N (q ) N (q ) N (q ) *g (p ) *g (p ) *g (p ) *g (p )
N (q ) N (q ) , N (q ) N (q )
*g (p ) *g (p ) . *g (p ) *g (p )
We then use the following formula to evaluate the control points: 1 Q " X \P > \ 3
(6.14)
6.1. Higher-order derivatives
*gK (p ) *gK (p ) 2 *gK (p ) H\K H H\K H H\K H K> *gK (p ) *gK (p ) 2 *gK (p ) H\K> H H\K> H H\K> H K> . $ $ $ $ *gK (p ) H H
This shows that Q can be evaluated from P if X and > are invertible. For example, assume n"m"3, j"6, to construct a linear system for tensor product B-spline surface based on partial derivatives in the u direction, we can write the following:
(6.8)
*fK (p , q ) *fK (p , q ) *fK (p ,q ) H L> L> 2 K> L> NL (q ) H\L H NL (q ) H\L H $
(6.12)
(6.11)
where
*fK (p , q ) H *fK (p , q ) H $
1 Q " X \P > \ n
where
H H *fK (p, q)" *gK (p) d NL (q) . (6.10) H V H VW W H VH\K WH\L Consequently, we obtain the following equation: P "nX Q >
In this case, we also have
From [1], the formula for calculating the second derivative of a polynomial curve F at u is stated as follows. F (u)"n ( f (u#1, u#1, u,2, u)!2f (u#1, u, u,2, u) #f (u, u, u,2, u))
(6.15)
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
where n "n(n!1) and in general nIM " n(n!1)2(n!k#1). Similarly, we rewrite Eq. (6.1) as
*F H (p;q)"m f (p,2, p, 1,1, q,2, q ) H GHI GFHFI *v L K !2f (p,2, p, 1,0, q,2, q ) H GHI GFHFI L K #f (p,2, p, 0,0, q,2, q ) . (6.16) H GHI GFHFI L K Each function f in the right-hand side can be expressed H as the following equations:
f (p,2, p; 1,1, q,2, q) H H H d gL (1,1, q,2, q) gK (p,2, p) " VW W H V H VH\K WH\L
where P "
*f(p , q ) H *f(p , q ) H $
*f(p , q ) H *f(p , q ) H $
*f(p ,q ) H K>
*f(p ,q ) H K>
f (p,2, p; 1,0, q,2, q) H H H " d gL (1,0, q,2, q) gK (p,2, p), VW W H V H VH\K WH\L
f (p,2, p; 0,0, q,2, q) H H H " d gL (0,0, q,2, q) gK (p,2, p). VW W H V H VH\K WH\L
(6.19)
Therefore,
$
2 *f(p ,q ) H K> L>
NK (p ) H\K H NK (p ) H\K H $
NK (p ) H\K> H NK (p ) H\K> H $
NK (p ) H\K H K>
NK (p ) H\K> H K>
d H\K H\L
d H\K H\L>
d Q " H\K> H\L $
d H\K> H\L> $
NK (p ) H H NK (p ) H H $
2 2 $
2 NK (p ) H H K>
2
d H\K H
2 d
H\K> H , $
$
d H H\L
d H H\L>
*g L (q ) H\L H *g L (q ) H\L> H $
*g L (q ) H\L H *g L (q ) H\L> H $
2
*g L(q ) H H
*g L(q ) H H
2
and (6.18)
*f(p , q ) H L> *f(p , q ) H L> $
2 2
X "
(6.17)
501
2
d
H H
*g L(q ) H H L>
(6.23)
From [1], the general formula for the kth derivative of a curve is
;gK (p,2, p). (6.20) VH Again, for convenience, we denote f (p,2, p; 1,1, q,2, q)!2f (p,2, p; 1,0, q,2, q)#f (p, H H H 2, p; 0,0, q,2, q) as *f(p, q) and W gL (1,1, q,2, q)! H W H 2gL (1,0, q,2, q)#gL (0,0, q,2, q)X as *gL (q). Hence, W H W H W H Eq. (6.20) can be simpli"ed:
H H *f(p, q)" d ) *gL (q) NK (p). (6.21) H VW W H V H VH\K WH\L Similar to Eq. (6.6), a GPV linear system for a TPBS can be constructed based on second partial derivatives as follows: P "m X Q >
*g L (q ) H\L H L> 2 *g L (q ) H\L >H L> . $ $
1 Q " ) X \P > \. m
!2gL (1,0, q,2, q)#gL (0,0, q,2, q)] W H W H
,
> "
Hence, for calculating the control points, we have
f (p,2, p; 1,1, q,2, q)!2f (p,2, p; 1,0, q,2, q) H H #f (p,2, p;0,0, q,2q) H H H " d [gL (1,1, q,2, q) VW W H VH\K WH\L
,
(6.22)
k FI(u)"nIM (!)I\H j WHWI ;f ( u#1,2, u#1, u,2, u ). (6.24) GFFHFFI GHI L\H H Thus, it is possible for us to construct linear systems for higher-order derivatives in a similar way.
7. Fast computation of basis function matrices In the previous sections, we have derived various formulations for constructing GPV-based linear systems. Usually, in each of these schemes, two matrices consisting of polar forms of B-spline basis function are involved. We call these matrices as basis function matrices. In the
502
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
following, we will present a fast computation algorithm for these basis function matrices. As stated in [6], the basis function gN (u ,2, u ) is computed recursively as GH P follows: g ( )"d G H G H and
(7.1)
u !t G gP\(u ,2, u ) gP (u ,2, u )" P G H P P\ t !t G H G>P G t !u P gP\ (u ,2, u ) # G>P> G> H P\ t !t G>P> G> (7.2) for 1)r)p. Eq. (7.1) states that g ( ) equals 1 when G H i"j and it equals zero otherwise while Eq. (7.2) depicts the recursive calculation process. We need to calculate each entry (i.e. gP (;)) in the basis G H function matrix. An example has been given in Eqs. (4.4) and (4.5). Let us take a look at the row Wg (; ), g (; ), g (; ), g (; )X of the basis function matrix in Eq. (4.5). In order to evaluate the entries of this row according to Eqs. (7.1) and (7.2), we perform recursive calculation for each basis function. This evaluation process is illustrated in Figs. 1a}d. Using this strict forward computation strategy, to evaluate all the gN in a row, the GH number of linear combinations required is evaluated by the following formula:
j"
p#2 p
!(p#1) ;(p#1)"
;(p#1)"[p(p#1)]/2
p#1
Fig. 1. This "gure compares the strict forward computation of a row in the basis function matrix (a)}(d) and a more e$cient algorithm (e).
p!1
(7.3)
where
p#1 p!1
is the number of nodes in the computation structure, except that of the bottom row and (p#1) is the number of entries in a row. Similar form of expression for the number of nodes in a triangular computation structure can be found in [13]. The main idea of our computation algorithm is the following. Instead of performing calculation separately for each entry of in a row, we consider the computation of a single basis function, which subsumes all the required calculations. Explicitly, for a row starting with the entry gN , we perform the recursive evaluation process of gN. GH GH Fig. 1e depicts our computation structure. We consider the whole triangular computation structure as four triangular regions as shown in Fig. 2. Obviously, the calculations involved in the upper triangular region are wasted since we do not really need to compute gN, at the top of GH the triangle. Furthermore, since the computations in the
Fig. 2. The computation structure is divided into three parts and only the central part is utilized in the calculation.
right and left triangular regions involved only zero operands, they are in fact not productive to the computation of the row entries. Fig. 1e highlights the productive computations by bold lines. If the linear combinations in these non-productive regions are eliminated, we have the following number of linear combinations:
u"
2p#1
p#1
!3
2p!1
p!1
(7.4)
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
where
2p#1 2p!1
Table 1 Percentage gain across di!erent values of p
and 3
p#1
p Percentage gain (%)
p!1
represents the total number of linear combinations in the whole triangular computation structure and the nonproductive triangular regions. With some arithmetic, the gain of our algorithm compared with the strict forward computation algorithm is formulated as the following equation: h"j!u"p(p#p#2)/2.
(7.5)
Gain in terms of percentage is calculated as follows: o"h/j"(p#p#2)/(p#1).
503
(7.6)
Table 1 shows the percentage gain across di!erent p values. Our overall algorithm is summarized using the following piece of C-like pseudo code, which is constructed by utilizing the above observations. /* Each U[i] is an array of yoating point numbers representing the arguments in a polar form.*/ Matrix}evaluation(int p, int curve}seg}index, array of U) + int row; for(row "1; row("p#1; row##) Row}evaluation(p, curve}seg}index, U[row]); , /* This function calculates a row of a basis function matrix. */ Row}evaluation(int p, int j, #oat u ,u ,2,u ) + Calculate all the multiplying factors (i.e. the values used for multiplying the basis functions) represented in bold lines in Fig. 1e; /* This requires p and j. */ Starting from g (u , u ,..., u ) which is equal to N N N one, calculate the sub-basis function values in the reverse direction according to Fig. 1e; , The above pseudo code only provides the kernel of the basis function matrix calculation. As there are minor changes in calculating di!erent basis function matrices when di!erent conditions are speci"ed (e.g. degree elevation or partial derivatives), appropriate modi"cations in the above code are required. Besides the above algorithm, we also observe that we can save computation in calculating *gL (q)" GH gL (1, q,2, q)!gL (0, q,2, q) in the partial derivativeGH GH based GPV system. As an example, the following equations show the evaluation of g (1, q, q). By the symmetric
2 88.9
3 87.5
4 88.0
5 88.9
6 89.8
property of polar form arguments, we have 1!t g (q, q) g (1, q, q)"g (q, q,1)" t !t t !1 # g (q, q) t !t where
(7.7)
t !q q!t g (q)# g (q, q)" g (q), t !t t !t q!t t !q g (q)# g (q, q)" g (q) t !t t !t t !q q!t g ( )# g ( ), g (q)" t !t t !t q!t t !q g ( )# g (q)" g ( ), t !t t !t t !q q!t g ( )# g (q)" g ( ). t !t t !t From the above formulas, we observe that when we calculate g (0, q, q), most of the computations have al ready done previously in "nding out the value of g (1, q, q). In fact, we only need to do computation for the last recursive step. After showing our fast basis matrix computation algorithm, in the next section, we will present the experimental results based on the various GPV-based schemes derived in the previous sections.
8. Experiments In this section, we will present some examples which the aforementioned GPV-based schemes are employed as control point speci"cation mechanisms. In all cases, we need to calculate a set of general polar values. Strictly speaking, we need to "nd out the piecewise polynomial (p.p.) forms of the underlying B-spline curves. However, it is quite a computationally complicated task. Rather, just to ful"ll the purpose of illustrating the techniques in practice, we employ a simpler way in performing our experiments. That is, we set an arbitrary value to the "rst entry of the matrix P. We then "ll up the other entries by multiplying it with an appropriate factor obtained from the given polar forms. For example, if the "rst entry
504
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
f (u , u , u ; v , v , v ) is set to a value C, then an entry f (u , au , u ; v , v , bv ) will set to the value a;b;C. The shape of the resulting control points will only di!er in scale with the ones, which are computed by actually calculating the exact values. In the following, we will give the parameters used for constructing the examples. In all the cases, the degree of the surface is set to 3 in both directions and the knot vector is set to (0, 0, 0, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 12, 12, 12) with 15 control points along each direction. The control points we want to specify lay on the central part, as illustrated in Fig. 3. In the following cases, j is set to 8. For the direct extension and the degree elevation approaches, we set ; "(3, 6, 5), ; "(3.9, 6, 5), ; "(3, 6, 7), ; "(3, 12, 5), < "(1, 2, 3), < "(1.3, 2, 3), < "(1, 2, 4.2), < "(1, 4, 3). For the direct extension scheme, based on the formulation shown in Eq. (4.14), we obtain
0.50 0.65 0.70 1.00
P"
0.65 0.85 0.91 1.30 0.70 0.91 0.98 1.40
,
Fig. 3. This "gure shows our 15;15 control point mesh structure. The central part is highlighted with circles and is speci"ed through general polar values schemes. In case of degree elevation approach, only the white circle control points are speci"ed while in the other two cases, all the circled control points are speci"ed.
1.00 1.30 1.40 2.00
X"
0.00 1.33 !0.33
0.00
0.00 1.03 !0.03
0.00
0.00 0.00
1.67 !0.67
!3.00 5.33 !1.33
!4.00
!3.70
15.00
13.9
!20.00
!18.6
10.00
9.4
>"
0.00
>"
,
5.33
!10.00 !8.00 6.00
5.00
Note that, in this case, we set 0.5 to the "rst entry in P. The value 0.65 in the next entry is calculated by multiplying it with 1.3 which can be observed from the di!erence between ; and ; . For the degree elevation scheme, based on Eq. (5.6), we have
0.50 0.65 0.70 1.00
P" 0.65 0.85 0.91 1.30 , 0.70 0.91 0.98 1.40
XM "
0.50 0.83 !0.33
0.35 0.83 !0.18 ,
!0.50 2.17 !0.67
!3.70
15.00
13.90
!1.60 !1.33 6.60
5.33
!20.00 !18.60 !10.00 !8.00
!1.60 !1.33 6.60
!4.00
10.00
9.4
6.00
5.00
.
In this case, we assume we raise the degree from 2 to 3 in u direction of the TPBS. For the partial-derivative-based scheme, we set p "5.50, p "5.78, p "6.05, p "6.325, q "5.00, q "5.25, q "5.50, q "5.75. Based on Eq. (6.6), we have
0.20 0.22 0.24 0.26
P "
0.23 0.26 0.28 0.31 0.27 0.29 0.32 0.35
,
0.30 0.34 0.37 0.40
0.02 0.48 0.48 0.02
X "
0.00 0.30 0.62 0.08 0.00 0.14 0.66 0.19
,
0.00 0.07 0.54 0.39
> "
0.00
0.01
0.04
0.09
0.17
0.22
0.21
0.14
0.00 !0.13 !0.21 !0.22
!0.17 !0.09 !0.04 !0.01
.
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
505
A dome-shaped control point mesh is constructed for the purpose of comparison (Fig. 4a). The corresponding TPBS is shown in Fig. 5a while the corresponding control point values (height) of the central part (i.e. circled parts of Fig. 3) are depicted in Fig. 6a. The results of the three GPV-based control point speci"cation schemes with parameter settings given above are shown in Figs. 6b}d. To translate the calculated control points to the top for the sake of better visualization, an o!set of 550 has been added to the calculated values. The control meshes with new control points at the central part and the corresponding TPBS are shown in Figs. 4 and 5.
9. Discussion In this paper, we have extended the GPV-based control point speci"cation method in [16] to tensor product B-spline surfaces and proposed two more alternatives in the speci"cation technique, which are based on degree elevation and partial derivative evaluation, respectively. In the previous sections, we focus on the evaluation of an individual region of a TPBS. In fact, we can potentially extend our formulation to control the whole surface. Here, we give the formulation in the curve level. For example, given a normalized cubic B-spline with knot vector ¹"(0, 0, 0, 0, 1, 2, 4, 5, 6, 6, 6, 6) (this con"guration of the curve is taken from [5]), based on (4.1) we write the following formula: P"AQ
(9.1)
where P"[ f (; ), f (; ), f (; ), f (; ), f (; ), f (; ), f (; ) f (; )]2 Q"[d , d , d , d , d , d , d , d ]2,
g (; ) g (; ) g (; ) g (; ) g (; )
A"
g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; ) g (; ) g (; )
The empty slots in A are "lled with zeros and U is G a suite of arguments. Other con"gurations set by changing the corresponding entries of P and A of the above formulations are of course possible. The control points Q can be obtained by Q"A\P if A is invertible. We can
Fig. 4. The shapes of various control point meshes produced using di!erent kinds of GPV-based scheme: (a) Mesh for comparison, (b) direct extension scheme, (c) degree elevation and (d) partial-derivatives-based schemes.
then follow the technique mentioned in Section 4.2 to extend this formulation to a TPBS. As shown in the experiment section, the drawback of this kind of control point speci"cation scheme is the following. Since the ultimate shape and location of the resulting control points is not easy to predict initially, the user may need to try out many di!erent parameters in the process of setting up a desirable control point mesh. Alternately, as illustrated in the previous section, the user may need to add o!sets to the resulting points to make it
g (; ) g (; ) g (; ) g (; ) g (; )
g (; ) g (; ) g (; ) g (; ) g (; ) g (; ) g (; )
look "t to the surface as a whole. The key to this problem is a technique for restricting the range of the values of the evaluated control points. This may require the overall shape information of the surface and that will be our next research work.
506
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
Fig. 5. The shapes of the corresponding tensor product surfaces generated by the meshes shown in Fig. 4: (a) Surface for comparison, (b) direct extension scheme, (c) degree elevation and (d) partial-derivatives-based schemes.
Fig. 6. The control point values produced by our schemes: (a) Control points for comparison, (b) direct extension, (c) degree elevation and (d) partial-derivatives-based schemes.
10. Conclusions Recently, a GPV-based control point speci"cation scheme is proposed in [16]. In this paper, we "rst construct linear systems of tensor product B-spline surfaces (TPBS) under the GPV-based framework. We name this
extension as the direct extension scheme. Subsequently, based on similar principles, we investigate the issue of extending the linear systems based on degree elevation and partial derivatives evaluation. As inherited from the original scheme [16], a large number of free parameters are provided in the direct extension and degree elevation schemes. In the partial derivative schemes, partial derivatives on a surface can be used as parameters, which is a desirable feature as slope information is important in controlling the shape of the surface. In all cases, we need to choose appropriate parameters by conforming to the completeness criteria. In other words, the related matrices should be invertible. Furthermore, in case of partial derivative scheme, the parameters should be chosen such that the partial derivatives to be evaluated are within the 2D polynomial segment under consideration. Finally, we introduce an e$cient way of computing the basis function matrices. We discover that the calculations of the entries of a row of a basis function matrix can be reduced if we combine their calculations together. Large computation reduction is found and has
S.W.-C. Lam / Computers & Graphics 24 (2000) 493}507
been summarized in Table 1. Finally, we further observe that partial results can be reused in the case of partial derivative scheme. Besides the research work suggested in the previous section, the idea of using general polar values as control points for surfaces raises other important issues for future research. They are the extensions of these formulations to other B-spline surface types, for example, triangular B-spline surfaces [9,20] and geometric continuous Bspline (i.e. b-spline) tensor product surfaces [21]. Triangular B-spline surfaces can model complex objects with non-rectangular topology while b-spline provides extra shape control parameters based on geometric continuity. Light has already thrown on how to construct those surfaces using blossoming [8,10] and we plan to extend our research work to those types of surfaces.
Acknowledgements The author thanks reviewer's valuable comments and Prof. Stephen Mann of University of Waterloo for his help in providing the technical report [22] through the Internet.
References [1] Ramshaw L. Bezier and B-splines as multia$ne maps. In: Earnshaw RA, editor. Theoretical foundations of computer graphics and CAD. Berlin: Springer, 1988. p. 757}76. [2] DeRose T, Goldman R. A tutorial introduction to blossoming. In: Hagen H, Roller D, editors. Geometric modeling methods and applications. Berlin: Springer, 1991. p. 267}86. [3] Seidel HP. An Introduction to polar forms. IEEE Computer Graphics and Applications 1993;13(1):38}46. [4] Ramshaw L. Blossoms are polar forms. Computer Aided Geometric Design 1989;6:323}58. [5] Seidel HP. Knot Insertion from a blossoming point of view. Computer Aided Geometric Design 1988;5:81}6. [6] Seidel HP. A new multia$ne approach to B-splines. Computer Aided Geometric Design 1989;6:23}32.
507
[7] Seidel HP. Computing B-spline control points. In: Straber W, Seidel HP, editors. Theory and practice of geometric modeling. Berlin: Springer, 1986. p. 17}32. [8] Dahmen W, Micchelli CA, Seidel HP. Blossoming begets B-spline bases built better by B-patches. Mathematics of Computation 1992;59(199):97}115. [9] Fong P, Seidel HP. An implementation of triangular B-spline surfaces over arbitrary triangulations. Computer Aided Geometric Design 1993;10:267}75. [10] Seidel HP. Polar forms for geometrically continuous spline curves of arbitrary degree. ACM Transactions on Graphics 1993;12(1):1}34. [11] Lee ETY. Computing a chain of blossoms, with application to products of splines. Computer Aided Geometric Design 1994;11(6):597}620. [12] DeRose TD, Goldman RN, Hagen H, Mann S. Functional composition algorithms via blossoming. ACM Transactions on Graphics 1993;12(2):113}35. [13] Liu W, Mann S. An optimal algorithm for expanding the composition of polynomials. ACM Transactions on Graphics 1997;16(2):155}78. [14] Goldman R. Blossoming with cancellation. Computer Aided Geometric Design 1999;16(7):671}89. [15] Goldman R. The rational Bernstein bases and the multirational blossoms. Computer Aided Geometric Design 1999;16(8):701}38. [16] Duchaineau MA. Using general polar values as control points for polynomial curves. Computer Aided Geometric Design 1994;11:411}23. [17] Eck M, Handenfeld J. Knot removal for B-spline curves. Computer Aided Geometric Design 1995;12:259}82. [18] Farin G. Curves and surfaces for computer aided geometric design: a practical guide, 4th ed. NewYork: Academic Press, 1998. [19] Mann S, DeRose T. Computing values and derivatives of bezier and B-spline tensor products. Computer Aided Geometric Design 1995;12:107}10. [20] Seidel HP. Triangular B-splines. In: Strasser W, Wahl F, editors. Graphics and robotics. Berlin: Springer, 1995. p. 149}62. [21] Barsky BA. In: Kunii TL, editor. Computer graphics and geometric modeling using beta-splines. In: Berlin: Springer, 1988. [22] Mann S, DeRose T, Winkenback G. Computing values and derivatives of bezier and B-spline tensor products. Research Report CS-93-31, Computer Science Department, University of Waterloo, 1993.