Computer Aided Geometric Design 20 (2003) 1–10 www.elsevier.com/locate/cagd
A shape preserving representation with an evaluation algorithm of linear complexity J. Delgado, J.M. Peña 1 Departamento de Mathemática Aplicada, Universidad de Zaragoza, Zaragoza, Spain Received 4 October 2002; received in revised form 15 November 2002; accepted 18 November 2002
Abstract We consider a blending basis for which we obtain an algorithm for the evaluation of polynomial curves with linear time complexity and we prove that it is a normalized totally positive basis. Therefore, it simultaneously satisfies efficiency and shape preservation. We also provide the corner cutting algorithm for obtaining the Bézier polygon from the control polygon with respect this basis. Related bases with additional properties are also considered. 2002 Elsevier Science B.V. All rights reserved. Keywords: Shape preservation; Corner cutting; Evaluation algorithms; Total positivity
1. Introduction (t)) (t ∈ I ) a Let U be a vector space of real functions defined on a real interval I and (u0 (t), . . . , un basis of U . If a sequence P0 , . . . , Pn of points in Rk is given then we define a curve γ (t) = ni=0 Pi ui (t), t ∈ I . The points P0 , . . . , Pn are called control points and the corresponding polygon P0 . . . Pn is called the control polygon of γ . In Computer Aided Geometric Design the functions u0 , . . . , un are usually nonnegative and ni=0 ui (t) = 1 ∀t ∈ [a, b] (i.e., the system (u0 , . . . , un ) is normalized) and in this case we say that (u0 , . . . , un ) is a blending system. The convex hull property is an important property for curve design: for any control polygon, the curve always lies in the convex hull of the control polygon. The convex hull property holds if and only if (u0 , . . . , un ) is a blending system.
E-mail address:
[email protected] (J. Delgado). 1 Partially supported by BFM2000-1253 Research Grant, Spain.
0167-8396/02/$ – see front matter 2002 Elsevier Science B.V. All rights reserved. doi:10.1016/S0167-8396(02)00190-5
2
J. Delgado, J.M. Peña / Computer Aided Geometric Design 20 (2003) 1–10
These geometric properties correspond to some properties concerning the collocation matrices of the system of functions. The collocation matrix of (u0 (t), . . . , un (t)) at t0 < · · · < tm in I is given by u0 , . . . , un (1) := (uj (ti ) i=0,...,m;j =0,...,n . M t0 , . . . , tm Clearly, (u0 , . . . , un ) is blending if and only if all its collocation matrices are stochastic (that is, nonnegative and such that the sum of the elements of each row is 1). In interactive design we want that the shape of a parametrically defined curve mimics the shape of its control polygon; thus we can predict or manipulate the shape of the curve by choosing or changing the control polygon suitably. In order to assure this property, we have to require more properties to the collocation matrices (1). A matrix is totally positive (TP) if all its minors are nonnegative (see (Ando, 1987)) and a system of functions is TP when all its collocation matrices (1) are TP. In case of a normalized totally positive (NTP) basis one knows that the curve imitates the shape of its control polygon, due to the variation diminishing properties of TP matrices (see (Peña, 1999)). In fact, shape preserving representations are associated with NTP bases. On the other hand, in order to be able to guide the curve and to put together several pieces of curves, it is desirable for the designer to have a precise control over what happens at the ends of the curve. This leads to the endpoint interpolation property: the first control point always coincides with the start point of the curve and the last control point always coincides with the final point of the curve. Although the Bernstein–Bézier form is the usual representation of a polynomial curve (see (Hoschek and Lasser, 1993; Farin, 1996)) and presents optimal shape preserving properties (see (Carnicer and Peña, 1993)), other polynomial representations can also be useful. The computational cost of the de Casteljau algorithm for a polynomial curve of degree n is quadratic (that is, of O(n2 ) elementary operations), but there are other evaluation algorithms useful in design with linear computational cost (that is, of O(n) elementary operations). This can lead to an important reduction of the computational cost when they are generalized for the evaluation of parametric surfaces. A first example of evaluation algorithm with linear complexity uses the Bernstein–Bézier representation and a simple predivision (see (Schumaker and Volk, 1986)). However, this algorithm cannot be expressed as a corner cutting algorithm. It is well known that corner cutting algorithms present good stability properties because each step is formed by convex combinations. In contrast, the evaluation algorithm for the Wang-Ball basis (see (Wang, 1987; Shi-Min et al., 1996; Dejdumrong and Phien, 2000)) is a corner cutting algorithm and it also has linear complexity. However, this basis is not NTP (see (Delgado and Peña, 2002)) for degree greater than 3. In this section we introduce another family of systems whose evaluation algorithms have linear complexity and form NTP bases. m (t)), m 2, be the system of polynomials on [0, 1] defined by: Definition 1. Let (c0m (t), . . . , cm
c0m (t) = (1 − t)m , cim (t) = t (1 − t)m−i , cim (t) = t i (1 − t),
m − 1, 1i 2 m+1 + 1 i m − 1, 2
m (t) = t m . cm
In addition, if m is even, m
m
cmm (t) = 1 − t 2 +1 − (1 − t) 2 +1 , 2
J. Delgado, J.M. Peña / Computer Aided Geometric Design 20 (2003) 1–10
3
and, if m is odd, cmm−1 (t) = t (1 − t)
m+1 2
2
cmm+1 (t) = 2
1 m+1 1−t 2 2
1 m+1 m+1 1 − t 2 − (1 − t) 2 , 2 m+1 m+1 − (1 − t) 2 + t 2 (1 − t). +
Let us observe that, for m = 2, the previous system coincides with the Bernstein basis on [0, 1]. For m > 2, it is given by the product of the initial m factors of (1 − t
t)·
1−t 0 · 0 0 0
1−t 0
t 0 0 1 0 1/2 0 0 0 0
t 1−t
1 − t 0 · 0 t 0
0 0 0 0 1/2 0 1 0 0 1−t
t 0 1/2 1/2 0 1−t
1−t 0 0 0 0 0 · 0 0 0 t 0
t 0 0 0 0 0
0 0 t 0 1 0 0 0 0
1−t 0 · 0 0
0 0 1 1 0 0
0 0 0 0 0 0 0 0 1 0 0 1−t
t 0 0 0
0 0 1 0 1 0 0 1−t 0 0 0 ···. 0 0 t
0 0 0 t
(2)
It is not obvious to check that this system is formed by linearly independent functions. It will m m m be a consequence of the main results of Section 2. Given the curve C (u) = i=0 Vi ci (u) = m (t))(V0 , . . . , Vm )T , we derive from the expression of the system given by (2) the following (c0m (t), . . . , cm evaluation algorithm, which consists of successive convex combinations (see Fig. 1).
Fig. 1. Evaluation algorithm for m = 5.
4
J. Delgado, J.M. Peña / Computer Aided Geometric Design 20 (2003) 1–10
Algorithm. (1) If m is odd, let (1 − u)V0 + uV1 , Vi+1 , i = 12 V m−1 + 12 V m+1 , V 2 2 V, i (1 − u)Vm−1 + uVm ,
i = 0, 1 i m−3 , 2 i = m−1 , 2 m+1 i m − 2, 2 i = m − 1,
otherwise, let (1 − u)V0 + uV1 , V i+1 , i = V Vi , (1 − u)Vm−1 + uVm ,
i = 0, 1 i m2 − 1, m i m − 2, 2 i = m − 1,
i i = 0, 1, . . . , m − 1). replace the original curve by a new curve C m−1 , with control points (V i (i = 0, 1, . . . , m) and go back to step 1. (2) m ← m − 1. If m 3 then let Vi = V (3) We obtain the value C m (u) as follows 1 , 0 + uV Q = (1 − u)V 2 , 0 + uV R = (1 − u)V m C (u) = (1 − u)Q + uR. m . In step 2 we have m −1 = V Remark 2. If we perform step 1 with m 3 even, then we get V 2 2 i (i = 0, 1, . . . , m), so that V m −1 = V m . Performing again step 1, now with m − 1 odd, we obtain Vi = V 2 2 m −1 = 1 V m −1 + 1 V m = V m −1 . Taking into account these facts, one can check that, when m is odd, the V 2 2 2 2 2 2 number of sums is 2m and the number of multiplications is 4m, and if m is even the number of sums is 2m − 1 and the number of multiplications is 4m − 2. In Section 2, we prove that this system is an NTP basis by constructing the corner cutting algorithm for obtaining the Bézier polygon from its control polygon with respect to this basis. This corner cutting algorithm allows us a stable conversion from our representation to the Bernstein–Bézier form. In Section 3, we also use the corner cutting algorithm of Section 2 to obtain corner cutting algorithms from the control polygons with respect to new NTP bases to the Bézier polygon. These new bases present better shape preserving properties than the NTP basis studied in previous sections, but their corresponding evaluation algorithms have greater computational cost.
2. Total positivity of the basis An elementary corner cutting is a transformation which maps any polygon A0 . . . An into another polygon B0 . . . Bn defined by one of the following ways: either Bj = Aj for j = i and Bi = (1 − λ)Ai + λAi+1 or Bj = Aj for j = i and Bi = (1 − λ)Ai + λAi−1 for some 0 < λ < 1. Any composition of elementary corner cuttings forms a corner cutting algorithm. We now state a result which can be proved with the proof of Lemma 2 of (Goodman and Said, 1991). This proof shows that an elementary corner
J. Delgado, J.M. Peña / Computer Aided Geometric Design 20 (2003) 1–10
5
cutting is associated to a bidiagonal stochastic TP nonsingular matrix relating both control polygons and that a corner cutting algorithm relates two control polygons by a product of such matrices, which is again an stochastic TP nonsingular matrix. This matrix is also a matrix relating both systems of functions and hence, if one system is a basis, the other system is also a basis. Lemma 3. Suppose that we have an NTP basis (b0 , . . . , bn ) of a vector space of functions U and another U . If there is a corner cutting algorithm which obtains B0 . . . Bn from system of functions (c 0 , . . . , cn ) in A0 . . . An whenever ni=0 Ai ci = ni=0 Bi bi then (c0 , . . . , cn ) is also an NTP basis. m (t)) introduced in Section 1 is an NTP In this section, we shall prove that the system (c0m (t), . . . , cm basis. This is a very convenient property from a design point of view because it permits to derive properties of the curve from properties of the control polygon (see (Goodman and Said, 1991; Carnicer and Peña, 1993; Peña, 1999)). For the proof, we shall distinguish the case m even and odd. In addition, we shall obtain for both cases the corner cutting algorithm from the corresponding control polygon to the Bézier control polygon. From now on, we shall use the following notations: m m , qn := n n i n n n−i n t (1 − t) , gi (t) := t n−i (1 − t)i , i = 0, . . . , n, fi (t) := i n−i n−i n + 1 + i k+i n i = 0, . . . , n. t (1 − t)n+1−k , Si (t) := k + i k=1
Let us start with an auxiliary result, whose proof is straightforward taking into account that Sin (t) = (1 − t)Sin (t) + tSin (t). Lemma 4. For each i ∈ {0, 1, . . . , n − 1} one has n+i+1 n+i+1 n (t) + Si+1 (t) + tgi+1 (t). Sin (t) = (1 − t)fi+1 m (t)) is NTP when m is odd and the proof The following result proves that the system (c0m (t), . . . , cm provides the corner cutting algorithm from the corresponding control polygon to the Bézier control polygon. 2n+1 (t)), t ∈ [0, 1], of Definition 1 is an NTP basis for all n 1. Theorem 5. The system (c02n+1 (t), . . . , c2n+1
Proof. Suppose p(t) =
2n+1
Ai ci2n+1 (t).
(3)
i=0
By Lemma 3, it is sufficient to obtain the corner cutting algorithm from A0 . . . A2n+1 to the Bézier polygon of the curve p(t). Let us prove that it corresponds to the following corner cutting algorithm: (1) B0 = A0 B2n+1 = A2n+1
6
J. Delgado, J.M. Peña / Computer Aided Geometric Design 20 (2003) 1–10
(2) A0i = Ai (i = 1, . . . , 2n) (3) For i = 1, . . . , n i + 12 · (n + 1) i−1 1 n+1 An + · Ai−1 Ain = n+1+i 2 n + 1 + i n+1 i + 12 · (n + 1) i−1 n+1 1 i i−1 A + An+1 An+1 = · 2 n+1+i n n+1+i For j = n − 1, . . . , i i 2n + 1 − j Ai−1 Ai + Aij = j 2n + 1 − j + i 2n + 1 − j + i j +1 2n + 1 − j i Ai2n−j + Ai−1 Ai2n+1−j = 2n + 1 − j + i 2n + 1 − j + i 2n+1−j Next j Bi = Aii B2n+1−i = Ai2n+1−i Next i It is sufficient to see that p(t) =
2n+1
Bi bi2n+1 (t),
(4)
i=0 2n+1 (t)) is the corresponding Bernstein basis. where (b02n+1 (t), . . . , b2n+1 Let us define 1 1 (5) M i := Ain + Ain+1 , i = 0, 1, . . . , n, 2 2 which coincide for all i. By Definition 1, (3) and (5), using that A0i = Ai for i = 1, . . . , 2n, and denoting A00 := A0 and 0 A2n+1 := A2n+1 , we have
p(t) =
n−1
A0k ck2n+1 (t) + A0n t (1 − t)n+1 + M 0 1 − t n+1 − (1 − t)n+1
k=0
+ A0n+1 t n+1 (1 − t) +
2n+1
A0k ck2n+1 (t).
(6)
k=n+2
Observe that
n+1 − (1 − t)n+1 − t n+1 1 − (1 − t)n+1 − t n+1 = (1 − t) + t n = qkn+1 t k (1 − t)n+1−k = S0n (t),
(7)
k=1
which, using the function cn2n (t) of Definition 1, can be written as cn2n (t) = S0n (t).
(8)
J. Delgado, J.M. Peña / Computer Aided Geometric Design 20 (2003) 1–10
7
Taking into account (7), (6) can be written as p(t) =
n−1
A0k ck2n+1 (t) + A0n f0n+1 (t)t + M 0 S0n (t)
k=0
+ A0n+1 g0n+1 (t)(1 − t) +
2n+1
A0k ck2n+1 (t).
(9)
k=n+2
Then, using Lemma 4, it can be proved by induction that, for i = 1, . . . , n, p(t) =
n−i
A0k ck2n+1 (t) +
k=0
+
n
n+1+i Ak−n+i fk−n+i (t) + M i Sin (t) k
k=n−i+1 n+i
k=n+1
n+1+i An+1+i−k gn+1+i−k (t) + k
2n+1
A0k ck2n+1 (t).
(10)
k=n+i+1
Then, by (10) for i = n, we obtain (4) and the result follows.
✷
Fig. 2 illustrates the corner cutting algorithm for m = 5 of the previous proof. When m is even, it can be proved in a similar but simpler way as in Theorem 5 (because it is not necessary to use the auxiliary constants M i ) that the corner cutting algorithm from the control polygon 2n 2n A0 · · · A2n of the curve p(t) = i=0 Ai ci (t) to the corresponding Bézier polygon is provided by (1) B0 = A0 Bn = An B2n = A2n (2) A0i = Ai (i = 1, . . . , n − 1, n + 1, . . . , 2n − 1) Akn = An (k = 1, . . . , n) (3) For i = 1, . . . , n − 1 For j = n − 1, . . . , i i 2n − j Ai−1 Ai + Aij = j 2n − j + i 2n − j + i j +1 2n − j i Ai Ai−1 Ai2n−j = + 2n − j + i 2n−j −1 2n − j + i 2n−j Next j Bi = Aii B2n−i = Ai2n−i Next i In conclusion we can state the following result corresponding to the even case. 2n (t)), t ∈ [0, 1], of Definition 1 is an NTP basis for all n 1. Theorem 6. The system (c02n (t), . . . , c2n
m Let us observe that the basis also satisfies the endpoint interpolation property: if γ (t) = m i=0 Ci ci (t) m m then γ (0) = C0 and γ (1) = Cm . It is a consequence of c0 (0) = 1, ci (0) = 0 for all i = 1, . . . , m,
8
J. Delgado, J.M. Peña / Computer Aided Geometric Design 20 (2003) 1–10
Fig. 2. Corner cutting algorithm for m = 5. m cm (1) = 1 and cim (1) = 0 for all i = 0, . . . , m − 1. However, we can observe in Fig. 2 that this basis does not satisfy the boundary tangent property. Let us recall that a basis (u0 , . . . , um ) satisfying the endpoint interpolation property also satisfies the boundary tangent property if the segments P0 P1 , Pm−1 Pm are, P u respectively, tangent to the curve γ (t) = m i=0 i i (t) at the endpoints P0 , Pm . In the next section we introduce new families of bases satisfying this property.
3. New families of NTP bases The Bernstein basis of degree m can be obtained by calculating the product of the initial m factors of 1 − t t 0 0 1−t t 0 · (1 − t t ) · 0 1−t t 0 0 1−t t 0 0 1−t t 1−t t 0 0 0 1−t t 0 0 0 · (11) ···. 0 0 1−t t 0 0 0 0 1−t t Let h be a positive integer. Preserving the initial h rows and last h rows of each matrix in (11) and changing the remaining rows by the corresponding rows of the corresponding matrix in (2), we obtain the following families of systems (observe that for h = 1 we obtain the family of Definition 1): m m (t), . . . , cm,h (t)), m 2 and h 1, be the system of polynomials on [0, 1] that Definition 7. Let (c0,h coincides with the Bernstein basis for 2 m 2h and that for m 2h + 1 is defined by: m i m (t) = t (1 − t)m−i , 0 i h − 1, ci,h i m−1+h−i h m m m−i − 1, t (1 − t) , h i ci,h (t) = 2 h−1
J. Delgado, J.M. Peña / Computer Aided Geometric Design 20 (2003) 1–10
9
i +h−1 i m+1 + 1 i m − h, t (1 − t)h , 2 i m i m t (1 − t)h , m + 1 − h i m. ci,h (t) = i m ci,h (t) =
In addition, if m is even, n+h−1 j 2h h m h j −h+1 c m ,h (t) = + t (1 − t)h t (1 − t) 2 h h − 1 j =2h +
n+h−1 j =2h
and, if m is odd,
j t j −h+1 (1 − t)h , j −h+1
m−1
1 +h h m+1 t (1 − t) 2 + cm−1 c m−1 ,h (t) = m−1 (t), 2 2 2 ,h h−1 m−1 + h m+1 1 m−1 m 2 t 2 (1 − t)h . c m+1 ,h (t) = c m−1 ,h (t) + m+1 2 2 2 2 m
2
In Section 1, we obtained the evaluation algorithm associated with the basis (cim )0im from the expression of the basis provided by using the initial factors of (2). In a similar way, we can derive m )0im . Besides, these systems the evaluation algorithms associated with each of the new systems (ci,h are NTP bases. In fact, given a parametric polynomial curve p(t) with control points Ai with respect of to the basis (cim )0im (p(t) = ni=0 Ai cim (t)), the basis corresponding to the control points Ah−1 i m the corner cutting algorithms of the proofs of Theorems 5 and 6 is precisely the basis (ci,h )0im h−1 m (p(t) = m ci,h (t)). This can be proved analogously to Theorem 5 (and Theorem 6), by replacing i=0 Ai n 2n (t) = cn,h (t) (which can be proved by induction on n h) formula (8) (which is equivalent to (7)) by Sh−1 and formula (10) by p(t) =
n−i+h−1
2n+1 Ah−1 k ck,h (t) +
k=0
+
n+i−h+1 k=n+1
n
n+1+i Ak−n+i fk−n+i (t) + M i Sin (t) k
k=n−i+h n+1+i An+1+i−k gn+1+i−k (t) + k
2n+1
2n+1 Ah−1 k ck,h (t)
k=n+i−h+2
are NTP bases and the corner cutting algorithm of for i = h, . . . , n. In conclusion, the systems the proofs of Theorem 6 and Theorem 5 can be interpreted as a corner cutting algorithm from the control . . . Ah−1 with respect polygon A0 . . . Am with respect to the basis (cim )0im to the control polygon Ah−1 m 0 h−1 m h−1 to the basis (ci,h )0im , followed by a corner cutting algorithm from A0 . . . Am to the Bézier polygon m )0im presents better shape preserving properties than (see Fig. 2). This also shows that the basis (ci,h m a basis (ci,h¯ )0im when h > h¯ (see (Carnicer and Peña, 1993)), although the computational complexity of the evaluation algorithm increases as h increases. The greatest complexity appears with the Bernstein basis (h m/2), which has optimal shape preserving properties (see (Carnicer and Peña, 1993)). m )0im (ci,h
10
J. Delgado, J.M. Peña / Computer Aided Geometric Design 20 (2003) 1–10
We can easily check that the families of bases of Definition 7 for h 2 satisfy the boundary tangent property. Among them, the family of bases with the evaluation algorithm of lowest computational cost satisfying the boundary tangent property corresponds to h = 2, but the computational cost of this evaluation algorithm is approximately (when m is large enough) twice the computational cost of the evaluation algorithm of the family of systems of Definition 1.
Acknowledgement We thank to an anonymous referee for motivating us to look for a shape preserving basis with an evaluation algorithm with linear complexity and satisfying the boundary tangent property. The basis m )0im of Section 3 satisfies all these properties. (ci,2 References Ando, T., 1987. Totally positive matrices. Linear Algebra Appl. 90, 165–219. Carnicer, J.M., Peña, J.M., 1993. Shape preserving representations and optimality of the Bernstein basis. Adv. Comput. Math. 1, 173–196. Dejdumrong, N., Phien, H.N., 2000. Efficient algorithms for Bézier curves. Computer Aided Geometric Design 17, 247–250. Delgado, J., Peña, J.M., 2002. On generalized Ball bases. Preprint. Farin, G., 1996. Curves and Surfaces for Computer Aided Geometric Design, 4th Edition. Academic Press, San Diego. Goodman, T.N.T., Said, H.B., 1991. Shape preserving properties of the generalised Ball basis. Computer Aided Geometric Design 8, 115–121. Hoschek, J., Lasser, D., 1993. Fundamentals of Computer Aided Geometric Design. AK Peters, Wellesley. Peña, J.M., 1999. Shape Preserving Representations in Computer Aided Geometric Design. Nova Science Publishers, Inc. Schumaker, L., Volk, W., 1986. Efficient evaluation of multivariate polynomials. Computer Aided Geometric Design 3, 149– 154. Shi-Min, H., Guo-Zhao, H., Tong-Guang, J., 1996. Properties of two types of generalized Ball curves. Computer-Aided Design 28, 125–133. Wang, G.-J., 1987. Ball curve of high degree and its geometric properties. Appl. Math.: A Journal of Chinese Universities 2, 126–140.
