Trimming analytic functions using right sided ... - Semantic Scholar
COMPUTER-AIDED DESIGN Computer-Aided Design 33 (2001) 813±824
www.elsevier.com/locate/cad
Trimming analytic functions using right sided Poisson subdivision GeÂraldine Morin, Ron Goldman* Computer Science Department, Rice University, 6100 Main Street, Houston, TX 77005-1892, USA Accepted 15 January 2001
Abstract Just as the Bernstein basis and BeÂzier points give us tools to manipulate polynomials, the Poisson basis and control points provide a new framework for investigating analytic functions. Using the Poisson representation, two different subdivision algorithms for analytic functions can be derived as extensions of the BeÂzier polynomial subdivision algorithm. The ®rst procedure generates an approximation algorithm, the second a trimming algorithm. Three applications of these two subdivision algorithms are developed here: an evaluation algorithm, a method for extending the convergence domain of the Poisson representation of an analytic function, and an intersection algorithm for analytic curves. q 2001 Elsevier Science Ltd. All rights reserved. Keywords: Analytic function; Poisson distribution; Subdivision; Approximation; Trimming; Intersection
1. Introduction The blending functions for BeÂzier curves and surfaces are provided by the binomial distribution. Although the Poisson distribution has long been known in probability theory to be a limiting case of the binomial distribution [1], in geometric modeling the idea of considering limiting cases of the polynomial setting only appeared later on [2,5,12]. In approximation theory, however, the Poisson distribution has already been used in the Szasz±Mirakian operator [8]. Analytic functions provide a richer collection of modeling tools than polynomials. In Computer Aided Geometric Design, the Bernstein basis is acknowledged to be more suitable for representing polynomials than the monomial basis. Analogously, we shall ®nd it more natural to use a Poisson series to describe analytic functions, rather than the usual Taylor expansion. The generalization from polynomials in Bernstein±BeÂzier form to analytic functions in the Poisson representation is quite natural and extensive [7], but convergence issues speci®c to the Poisson representation need to be addressed. Moreover, to provide useful algorithms, we need to consider ®nite procedures. Based on the Poisson representation of an analytic function, we shall present two subdivision algorithms for analytic curves. These algorithms are limiting cases of the BeÂzier subdivision algorithm for polynomial curves [4]. The Poisson representation and these associated subdivision algorithms provide new tools for representing and * Corresponding author. E-mail address: [email protected] (R. Goldman).
manipulating analytic functions intuitively and ef®ciently. The ®rst part of this paper establishes the theoretical foundations and addresses convergence issues. The second part concentrates on several applications of these tools: representing an analytic function both inside and outside the convergence domain of its Taylor series, evaluating an analytic function or its derivative at a point, and, ®nally, intersecting analytic curves. The Poisson representation has already been exploited in [12], where the Poisson representation for analytic functions is introduced along with a subdivision algorithm that generates a rendering procedure for analytic functions inside the domain of convergence of their Taylor series. Here, we shall begin by recalling this work, although our de®nition of a Poisson function differs slightly. We shall then present a new subdivision procedure, together with corresponding analysis algorithms. In the last part of this paper, we shall develop some applications of these new procedures.
2. The Poisson representation 2.1. The Poisson basis The Poisson basis is an in®nite collection of functions that we shall denote by bk tk$0 . If Bnk t denotes the k-th Bernstein polynomial of degree n, then: bk t limn!1 Bnk
0010-4485/01/$ - see front matter q 2001 Elsevier Science Ltd. All rights reserved. PII: S 0010-448 5(01)00099-9
t ; n
k $ 0:
1
814
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
Moreover, the convergence is uniform on any bounded interval of [0, 1) (for a proof see the appendix of [12]). The Poisson basis functions can also be expressed explicitly: for any non negative integer k; bk t e2t tk =k!: These basis functions are positive and form a partition of unity on [0, 1), just as the Bernstein basis functions do on the interval [0, 1]. 2.2. Poisson curves We call F a Poisson function or a parametric Poisson curve on the interval [0, R), if for any x0 [ [0, R) the function F is de®ned by a Poisson series in a neighborhood U of P x0: that is, for any t [ U, F t k$0 Pk bk t 2 x0 converges. 1 The vector of coef®cients p (P0, P1, ¼) is the vector of Poisson control points of F at x0.
3. Poisson subdivision 3.1. BeÂzier subdivision Poisson subdivision is an extension of the well-known de Casteljau subdivision algorithm for BeÂzier curves. In this section, we shall recall some classical results for polynomials, in order to be able to refer to them later on. A BeÂzier curve is a polynomial P parameterized over an arbitrary interval [a, b]. If n is the degree of the curve, the control points (P0, P1, ¼, Pn) of the curve are the coef®cients of P in the basis t2a Bnk ; b 2 a 0#k#n that is: P t
2.3. Representation of analytic functions in the Poisson basis A Poisson curve F is de®ned by its Poisson control points at 0 over the interval [0, R) if and only if the Taylor series of F P at 0 converges on this same interval. Indeed, if F k$0 Pk bk on the interval [0, R), then F is P the product of the function t 7 ! e2t and the Taylor series k$0 Pk tk =k! et F t which converges, by hypothesis, over the P interval [0, R). Conversely, it is easy to see that if F t Pk$0 Ak tk =k!; t k then there P exist (P0, P1, ¼) such that e F t Pk t =k!; so F t k$0 Pk bk t: Moreover, there are simple and ef®cient algorithms to convert between the Taylor coef®cients Ak k$0 and the Poisson coef®cients P Pk k$0 [7,12]. The Taylor series F t k$0 Ak tk =k!; converges uniformly to F on P any subinterval [0, b] of [0, R). Similarly, the Poisson series k$0 Pk bk t converges uniformly to F on any interval [0, b] , [0, R). However, if R 1 and limt!1F(t) 0, then, the Poisson series converges uniformly to F on [0, 1) [14]. We shall illustrate the difference of convergence between the Taylor and the Poisson series in the application section. Any function de®ned by a Taylor series over the interval [0, R) is C 1[0, R) and its derivatives are also de®ned by a Taylor series on [0, R) generated by differentiating the original series term by term. Therefore, a Poisson function, de®ned by a Poisson series on [0, R), is also C 1[0, R) and its derivatives are also Poisson functions de®ned by Poisson series generated by differentiating the original series term by term. Moreover, since b 0k t bk21 t 2 bk t; it is easy to show that the vector of Poisson control points of the m-th derivative, F (m), of a Poisson function F with vector of control point p is the m-th discrete difference of p. Thus, if D denotes the difference matrix, i.e. Di,i 21, Di,i11 1, otherwise Di, j 0, then D mp are the control points of F (m). 1 The authors want to thank Jean Gallier for suggesting this de®nition of a Poisson curve.
n X k0
Pk Bnk
t2a : b2a
The polygonal curve de®ned by the control points is called the control polygon of the BeÂzier curve, and is independent of the interval of parameterization [a , b ]. Therefore, a common choice is to take [a , b ] [0, 1]. Here, however, we shall consider the parameterization over the interval [0, n] for two reasons: ®rst, the Poisson basis is a limiting case of the Bernstein basis of degree n over the interval [0, n]; second, we are going to extend the subdivision algorithm in two distinct ways by considering two different approaches to subdivision, but, over the interval [0, 1], these two approaches merge into one. Geometrically, a step of subdivision corresponds to splitting a BeÂzier curve into the union of two BeÂzier curves, as illustrated in Fig. 1; that is, from the control polygon of P, one round of subdivision generates two new control polygons, Qk 0#k#n and Rk 0#k#n ; corresponding to these two new BeÂzier curves. Iterating the subdivision process produces a sequence of control polygons that converges to the corresponding BeÂzier curve [10]. In the polynomial setting, the de Casteljau algorithm is both an evaluation algorithm and a subdivision algorithm (Fig. 2). This dynamic programming procedure offers a nice way to visualize, and an ef®cient way to compute, subdivision. In the following discussion, we shall consider a BeÂzier curve P of degree n parameterized over the interval [0, n]. The subdivision parameter can be de®ned in two different ways: either we can choose r in (0, 1) and subdivide at nr, in
Fig. 1. A BeÂzier curve of degree 3, with its initial control points (on the left) and its control points after one round of subdivision (on the right). The ®rst four of the new control points are generated by the left subdivision algorithm, the last four by the right subdivision algorithm.
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
815
braic identities for the change of basis are: X k a t a n2j t n a 1 12 Bk2j Bnj ; 0 # k # n: Bk n n n n n j0 4 From these equalities, we obtain an expression for the new control points: nX 2k a Rk a Bn2k 0 # k # n: 5 P ; j n k1j j0 Fig. 2. The de Casteljau algorithm for BeÂzier curves. A point emerging from two arrows is an af®ne combination of the points at the base of these arrows, e.g. Q1 (1 2 .)P0 1 .P1.
which case the coef®cients along the arrows in the de Casteljau algorithm would be 1 2 r and r; or we can choose an absolute parameter a, in which case the coef®cients in the de Casteljau algorithm would be 1 2 a/n and a/n. We shall extend both of these approaches to subdivision to the Poisson setting: the ®rst approach generates an algorithm for approximating a Poisson curve; the second approach yields both an evaluation algorithm and a trimming algorithm. In the next paragraph, we present algebraic identities characterizing both of these aspects of subdivision in the polynomial setting. In the ®rst approach, subdivision occurs at the parameter rn where 0 , r , 1. Left-sided subdivision is equivalent to a scaling of the parameter interval: from the control points (Pk) of the curve over the interval [0, n], we compute the control points (Qk) of the restriction of the curve to the interval [0, rn]. The following algebraic identities, from [6], express the Bernstein basis functions over the interval [0, rn] in terms of the original Bernstein basis functionsÐthat is, the basis functions de®ned over the interval [0, n]: Bnk
rt n
n X jk
Bjk rBnj
t ; n
In the next two sections we shall present generalizations of these two subdivision algorithms to the Poisson setting, and in the remainder of the paper we shall develop several applications of these algorithms for analytic functions. 3.2. Approximation: left-sided Poisson subdivision In this section, we brie¯y recall some of the results presented in [12], which we shall use later on to justify certain convergence properties. In the following discussion, F will denote a Poisson curve de®ned by the control points (P0, P1, ¼) (at 0) over the interval [0, R), that is: X Pk bk t; 0 # t , R: F t k$0
The subdivision scheme that we shall describe shortly, induces an approximation algorithm for any such function F. This scheme is an extension of the corresponding polynomial subdivision algorithmÐthat is, the left side of the BeÂzier subdivision algorithm over the interval [0, n] at the parameter rn, where 0 , r , 1. Eqs. (2) and (3) de®ning left subdivision in the polynomial case generalize to the Poisson setting: X j bk rt Bk rbj t; k $ 0; 6 j$0
and 0 # k # n:
2
F rt
X k$0
From these equalities, follow expressions for the points Qk as functions of the points Pk : Qk r
k X j0
Bkj rPj ;
0 # k # n:
3
In the same way, we can derive algebraic identities modeling right-sided subdivision [6]. In this second approach, the subdivision parameter is a number a in [0, n]. One step of this right-sided subdivision is equivalent to computing the control points of the restriction of the curve to be interval [a, n] from the original control points over the interval [0, n]. The corresponding alge-
Qk rbk t with Qk r
k X j0
Pj Bkj r;
k $ 0: 7
Although recursive subdivision converges only for r in (0, 1), equalities (6) and (7) hold for arbitrary real values r as do Eqs. (2) and (3) in the polynomial setting. Moreover, the new control points (Qk(r))k$0 can be generated by applying an extension of the polynomial de Casteljau algorithm to the original control points (Pk)k$0 (Fig. 3). The following theorem, proved in [12], establishes the uniform convergence, as r goes to zero, of the control polygon Q r Q0 r; Q1 r; ¼ to the analytic curve F. Theorem 1. Let L(F, r) denote the piecewise linear function de®ned by the control polygon Q(r), generated after one
816
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
approximation to the Poisson curve de®ned by F over the interval [a, R). These approximations can be re®ned, as desired, using the approximation algorithm presented in Section 3.2.
Fig. 3. The left subdivision algorithm for Poisson curves.
round of subdivision at r, so that L(F, r)(jr) Qj(r), for any integer j, and L(F, r) is linear on [jr, (j 1 1)r], for all j in N. If the sequence of subdivision parameters (ri) converges to zero, then the sequence of piecewise linear functions L(F, ri) converges uniformly to F over any closed subinterval of [0, R). Thus, left-sided subdivision generates a sequence of control polygons converging to a Poisson curve. This convergence provides a rendering algorithm for Poisson curves on any interval [0, b] , [0, R). 3.3. Trimming: right-sided Poisson subdivision 3.3.1. Extension of the polynomial procedure In this section, we shall extend right-sided polynomial subdivision to analytic functions. Equalities (4) and (5) from the polynomial case, generalize readily to the Poisson setting. Indeed, it follows easily from the binomial theorem that: bk a 1 t
k X j0
bk2j tbj a;
so, for a in (2R, R): F a 1 t
X
k$0
Rk abk t where Rk a
8
X j$0
Pk1j bj a;
3.3.2. Approximation of right-sided subdivision Eq. (9) de®nes the control points (Rk(a))k$0 of the Poisson curve F restricted to the interval [a, R), but no algorithm to ef®ciently compute these control points is provided. In fact, it is not possible to compute these points by extending the de Casteljau algorithm (Fig. 2) directly, as in the case of left-sided subdivision. Since in the Poisson setting we start with an in®nite number of control points, the points corresponding to the points on the right lateral edge of the de Casteljau diagram in the polynomial case recede to in®nity. The in®nite sum de®ning these points in Eq. (9) expresses this same phenomenon. To obtain ®nite sums, for every positive integer k we shall de®ne a computable sequence Rnk an$0 so that the Poisson curves Fan de®ned by the control points Rnk ak$0 converge uniformly to the Poisson curve Fa, that is, to the trimmed curve with control points Rk ak$0 : Since the scaled Bernstein basis converges to the Poisson basis, we de®ne: n X a Rnk a Pk1j Bnj ; k $ 0; n $ 0: 10 n j0 Once again, the de Casteljau algorithm allows us to visualize the construction of the points Rnk ak$0 and de®nes a dynamic programming algorithm to generate these points from the initial control points (Pk)k$0 (Fig. 4). We shall use the following lemma to prove that Fan is well de®ned and that the sequences Rnk an$0 and Fan n$0 converge to Rk a and Fa , respectively. P Lemma 2. Let F t Pk bk t be a Poisson curve over the interval P [0, R). Then for any positive integer k, the series i$0 Pi1k bi t de®nes a Poisson curve over
k $ 0:
9 These equations provide a tool for trimming Poisson curves, since the points (Rk(a))k$0 are the Poisson control points of the curve parameterized by the analytic function Fa F(a 1 .). We shall discuss, in detail, in the applications section, the domain where the Poisson series at 0 of Fa convergesÐthat is, where we can approximate Fa by our left-sided subdivision algorithm. In this section, though, we shall consider only the interval 2R 1 uau; R 2 uau; where the function Fa is de®ned by the control points Rk ak$0 : if t [ 2R 1 uau; R 2 uau; then a 1 t [ (2R, R) and, therefore, Eq. (9) converges. Geometrically, the points (Rk(a))k$0 de®ne a control polygonÐthat is, a piecewise linear
Fig. 4. A de Casteljau-like algorithm to compute the points Rnk a:
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
P
the interval [0, R), and X i$0
Pi1k bi
k X j0
Cjk F j ;
11
where Cjk denotes the binomial coef®cient k!= k 2 j!j!: Proof. This result is P proved by induction on k. First, for k 0: by hypothesis i$0 Pi bi F, which is a Poisson function on the interval [0, R). For k 1, observe that: b 0k bk21 2 bk : Using this identity, we can easily deduce that: X Pi11 bi F 1 F 0 : 12 i$0
From Section 1, both F and F 0 are Poisson functions over the P interval [0, R), and, by addition, so is the function i$0 Pi11 bi : Eq. (12) yields, by induction: X i$0
Pi1k bi
k X j0
Cjk F j :
13
Thus, by induction, for any positive k, the function P i$0 Pi11 bi is a Poisson function on the interval [0, R).A Now consider the points (Rkn(a)) and (Rk(a)) all de®ned for a ®xed parameter value a in [0, R). Theorem 3.
limn!1 Rnk a
Rk a
Proof. Note that, if we call and R0(a) the ®rst points of F, then Rnk an$0 and Rk a are the ®rst points of the function with Poisson coef®cients Pm m$k . From Lemma 2, this function is a Poisson function over the interval [0, R). Thus, it is suf®cient to prove that limn!1 Rn0 a R0 a, for a generic Poisson function, that is, for the function F. For a in [0, R), it follows from Theorem 1 with rn a=n that: a lim L F; a F a: n!1 n P But, by construction L F; a=n aP ni0 Pi Bni a=n Rn0 a and by de®nition F a i$0 Pi bi a R0 a. So limn!1 Rn0 a R0 a.A Fan ;
Before proving the convergence of the sequence we P need to verify that the series Rnk abk t converges, that is, the function Fan t is well de®ned: n X n XX a Rk abk t Pk1j Bnj b t n k k$0 k$0 j0
j0
Bnj
X a P b t n k$0 j1k k
But, from Lemma 2, P Pj1k bk t converges on [0, R). Therefore, the in®nite series Rnk abk t Fan t is well de®ned. The following theorem establishes the convergence of the functions Fan ; de®ned by the sequences Rnk ak$0 ; to the function Fa , de®ned by the sequence Rk ak$0 ; where a [ (2R, R). Note that this proof does not use Theorem 3 in which a had to be positive. Theorem 4. As n goes to in®nity, the sequence Fan converges uniformly to Fa on any closed subinterval of 0; R 2 uau. Proof. Let (P0, P1, ¼) be the vector of control points for the Poisson function F. We apply on this initial set of control points one round of the subdivision algorithm at the parameter a/n as illustrated in Fig. 4. Denote by Qn0 a; Qn1 a; ¼ the sequence of points generated by left subdivision: that is, the points that appear on the left lateral edge of the algorithm. In addition to F and Fan ; we now introduce the Poisson functions Qn and Qna de®ned by the control points Qn0 a; Qn1 a; ¼ and Qnn a; Qnn11 a; ¼: From Section 3.2, on left-sided subdivision, we have: a n t ; 14 Q t F n and Qna t
Fan
a t : n
15
Moreover, from the proof of Lemma 2:
(R0n(a))n$0
n X
817
Qna t
n X k0
Ckm Qn k t:
From Eqs. (14) and (16), we get: n X n! ak k a t : F Qna t k k! n k0 n 2 k!n
16
17
Thus, by Eqs. (15) and (17) and using the change of variable x (a/n)t: Fan x
n X k0
n! ak k F x: n 2 k!nk k!
18
But by Taylor's theorem: Fa x F a 1 x
X ak k F x: k! k$0
19
From Eqs. (18) and (19) it appears already that Fan x converges pointwise to Fa(x), but to establish uniform convergence we must show that k n X n! k a F F a 1 x 2 k! n 2 k!nk k0
goes to zero as n goes to in®nity for every x in
818
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
0; r , 0; R 2 uau. For any integer n0, we have:
This last inequality implies:
n0 n X X n! ak ak k k F x x # F a 1 x 2 F F a 1 x 2 k k! k! k0 n 2 k!n k0
n k X 0 n! a k e F x # ; 12 k k! 3 n 2 k!n k0
n k X 0 n! a k F 1 12 x n 2 k!nk k! k0 X k n n! a k : 1 F x kn 11 n 2 k!nk k!
27
for all x [ [0, r ] and n $ n4.Using Eqs. (21) and (24) with m n0 and Eq. (27) for any n . n4, we can bound each term of the left-hand-side of Eq. (20) by e=3: Thus, Fan converges uniformly to Fa on [0, r ].
0
P
20
Let us ®x e . 0. The series F a 1 x k$0 ak =k!F k x converges uniformly relative to x on any closed subinterval [0, r ] of 0; R 2 uau (see the proof in the Appendix). That is, there exists n1 so that, for all m $ n1 and x [ [0, r ]: m X ak k e F x , : 21 F a 1 x 2 k! 3 k0 and, from the absolute and uniform convergence of this same sum, there exists n2 so that for all n $ m $ n2: n X ak k e uF xu , : k! 3 km
22
Moreover, for any positive integers n and k: 0,
n! , 1; n 2 k!nk
23
so from Eqs. (22) and (23), for all m $ n2: n n X X n! ak k ak k e # F uF xu , : x k k! k! 3 km 1 1 n 2 k!n km 1 1
The following corollary extends Theorem 4 to the derivatives of Fa: Corollary 5. For any integer k $ 0, the functions Fan k converge uniformly to Fa k on [0, r ], where r is an arbitrary number in 0; R 2 uau: Proof. This corollary is a consequence of a general result on convergence for analytic functions, that is, if we have a uniform approximation of an analytic function by analytic functions, then the k-th derivative of the approximants also provides a uniform approximation of the k-th derivative of the original function. See Davis [3], p. 112.A In this section, we constructed the Poisson control points of a trimmed Poisson curve, but, since these points were not computable, they were expressed as the limit of order n approximating control points. We then proved that the Poisson functions de®ned by these approximating families of points converge uniformly to the trimmed functions. This convergence property is true for any function de®ned by its Poisson expansion at zero, as well as for all its derivatives.
24 We know (for a detailed proof see the appendix of [12]): n! # k k 2 1 : 1 2 25 k 2n n 2 k!n P Let n0 max(n1, n2). Since ak =k!F k x converges absolutely and uniformly on [0, r ], there exists a constant M . 0 so that, for all x [ [0, r ]: n0 k X a k 26 F x # M: k! k0 Choose n4 so that n4 $ M
n0 n0 2 1 3 : 2 e
From Eqs. (25) and (26), for any x [ [0, r ] and n $ n4: k n0 a X n! e j F max j1 2 x # : j j 0#j#n0 k! 3 n 2 j!n k0
4. Applications Here, we ®rst compare the approximation of an analytic function F by its Taylor polynomial and its Poisson partial sum with the same number of terms. In the remainder of this section we shall use the two subdivision algorithms introduced in Section 3 to develop practical analysis algorithms for Poisson curves and surfaces. First, we present an evaluation procedure for Poisson curves via an approximation algorithm that comes about as a direct consequence of right-sided subdivision. Then, we shall propose a rendering method for analytic curves, on different intervals. In the case of analytic functions with poles, we shall consider more closely than in Section 3 the domain of convergence of a trimmed Poisson function; the domain where a Poisson curve can be approximated using a subdivision algorithm can actually be extended thanks to the trimming algorithm. Finally, we shall describe an intersection algorithm for analytic functions, similar to the intersection algorithm induced by subdivision for polynomial curves.
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
819
Fig. 5. The left ®gures show the functions e2t sint at the top and sin t at the bottom (solid) together with their Taylor (long dashes) and Poisson (short dashes) approximations of order 10. The right ®gures show the difference between the function and the corresponding Taylor (long dashes) and Poisson (short dashes) approximations.
4.1. Comparison of Taylor and Poisson approximations of analytic functions The convergence of a Poisson series to a function F such that limt!1 F t 0 is uniform on [0, 1) (Section 2.3), whereas the convergence of the Taylor series is uniform only on ®nite intervals. As an illustration, Fig. 5 (top) shows the function e 2t sin t and its approximations by Taylor and Poisson sums of order 10Ðthat is, the ®rst 11 terms. The Taylor approximation is a polynomial and diverges rapidly (in this example around 2.5). Moreover, the fact that the Taylor approximation diverges and the Poisson approximation is bounded also gives an advantage to the Poisson series for approximating bounded functions. The bottom part of Fig. 5 illustrates the approximation of the function sin t by the ®rst 11 terms of the Taylor and Poisson series.
and R0n(a) can be computed by applying the de Casteljau algorithm at a/n (Fig. 4). The derivatives of F can also be approximated using the de Casteljau algorithm. From Section 2.3, the coef®cients of the m-th derivative of F can be derived from the coef®cients of F by taking discrete differences (Fig. 6). Then, using the same procedure for F m as for F, we can ®nd a sequence of values converging to F m a for any a in 0; R: 4.2.2. Approximation on a domain [0, b), where b . 0 The approximation algorithm for a Poisson curve on an interval [0, b], where 0 , b , R, is a direct application of
4.2. Approximation of a Poisson function inside the convergence domain of its Poisson series 4.2.1. Evaluation of a Poisson curve and its derivatives at a given parameter value As in the polynomial setting, the de Casteljau algorithm provides an evaluation procedure for analytic curves. At a given parameter a in the domain of convergence, F a R0 a; but, from Theorem 3, we have: lim
n!1
Rn0 a
R0 a;
Fig. 6. From this dynamic programming algorithm, the Poisson control points P m i of the derivative of order m of an analytic function are computed from the Poisson control points Pi of the function.
820
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
Fig. 7. The function F x 1= 1 1 x (solid) and successive approximations using the left-sided Poisson subdivision algorithm of order 1 (clearly piecewise linear curve), 4 (close to the original curve) and 8 (overlapping with the original curve for values less than 1.5 and then diverging).
the left subdivision algorithm and has already been presented and detailed in [12] where, in particular, a justi®cation that only a ®nite number of control points is necessary to compute the approximation is presented. For all the applications given below the subdivision parameter is r 1/2, and the order of the subdivision is the number of iterations of subdivision at r 1/2. But, one step of subdivision at 1/2 i is equivalent to i steps of subdivision at 1/2 [12] and, if the intermediate results are not needed, is more ef®cient. Therefore, we actually compute i rounds of subdivision at r 1/2 by applying one round of subdivision at 1/2 i. Fig. 7 represents approximations of the function 1/(1 1 x) by several control polygons generated from subdivision. This function has a pole at 21; the convergence result for the approximation algorithm holds only in the domain of convergence of the Taylor seriesÐthat is, in this example, for values less than 1. We shall see in Section 4.2 how to approximate this function for values larger than 1. The approximation algorithm extends to surfaces. Fig. 8 illustrates the use of the subdivision algorithm to generate approximations of the sphere. Note that this algorithm is the only known stationary subdivision scheme providing an exact representation of the sphereÐthat is, such that the sphere is the limit of the subdivision process. To compute the initial control
points (Pjk) of the sphere cosucosv; cosusinv; sinu; u [ [0, 2p], v [ [0, 2p], we ®rst generate the Poisson control points ck and (sk) of the trigonometric functions cos t and sin t from their Taylor coef®cients by the algorithm proposed in [12]. Then, the left-sided subdivision algorithm is applied to Pjk j;k$0 cj ck ; cj sk ; sj j;k$0 in both parameter directions. That is, ®rst, one iteration of subdivision is applied to the points Pjk j$0 for each k, to generate the new control points Q~ jk j$0 . Then, subdividing the family of control points Q~ jk k$0 for each j generates the new control points Qij : As with the curves, the subdivision parameter is 1/2 and the order of subdivision is the number of iterations in each direction. The approximation algorithm presented in this paragraph generates arbitrarily good piecewise linear approximations of a Poisson curve on an interval originating at 0 or of a tensor product Poisson surface on a domain originating at (0, 0). In the next section, we shall explore how the trimming algorithm allows us to relax this domain constraint. 4.2.3. Approximation on a domain not containing zero Using both the approximation and trimming algorithms presented in Sections 3.2 and 3.3, we can construct a more general approximation procedure; that is, we can approximate a Poisson function F on a; b; an arbitrary closed subinterval of 2R; R where R is the radius of convergence at zero. From the control points, Pk ; of F we can ®nd a control polygon Rnk a that de®nes a function Fan approximating the trimmed function Fa : We can then apply the approximation algorithm induced by left subdivision to the control points Rnk ak$0 : From Theorem 4, the convergence of the function Fan is uniform on the interval [a, b], and by Theorem 1 the control polygons generated by the leftsided subdivision algorithm also converge uniformly to the function Fan on the interval [a, b]. Therefore, combining these two algorithms, we can construct a sequence of control polygons converging uniformly to the function F on the interval [a, b]. Fig. 9 illustrates this algorithm. This approximation depends on the two subdivision algorithms: for the trimming algorithm, the relevant parameter is the degree4 n, chosen here to be powers of two; in our example F12 is the trimmed function generated for n 2 4 16. For the approximation algorithm, the order,
Fig. 8. A Poisson subdivision scheme converging to the sphere (cosucosv; cosusinv; sinu), u [ [0, 2p], v [ [0, 2p]: control polyhedra of the sphere of orders 1, 2, 3 and 6.
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
Fig. 9. Approximation of the function F(x) cos x (solid line) over the interval [1,4 4] by 8 the control polygons of order 5 of the functions 1 2 F12 ; F12 ; F12 and F12 (dashed lines).
as in the previous section, is the number of iterations at 1/2. In Fig. 9 the order of the approximation algorithm is kept constant in order to observe the convergence of the trimming algorithm. The algorithm presented above is general enough to approximate entire functionsÐthat is, analytic functions having no polesÐon any interval [a, b]. Nevertheless, for an analytic function having a ®nite domain of convergence (2R, R) around zero, we are limited to intervals [a, b] , (2 R, R). In the next section, we shall see how to approximate outside this domain of convergence. In particular, this new method will solve the divergence problem illustrated in the example in Fig. 7. 4.3. Approximation of a Poisson function outside the initial convergence domain Here, we address the approximation of analytic functions with poles. So far, we are able to approximate functions with poles on any closed subinterval of their domain of convergence around zero. Below, trimming is used to extend the domain where a Poisson function is de®ned by a Poisson series. The rendering process induced by left subdivision then applies on this extended domain.
821
Fig. 10. Approximation of the function F x 1= 1 1 x over the interval [1/2, 2) 2(solid4 line) by the control polygons of order 5 of the functions 20 2 2 F1=2 ; F1=2 ; F1=2 (dashed lines). Note that the singularity appearing in Fig. 7 for x . 1 has been eliminated.
tion; that is, the radius of convergence Ra of the Poisson series de®ning Fa may be greater than R 2 uau. Moreover, a Poisson series is associated with a rendering algorithm on the interval [0, Ra). Thus, the function F can be rendered on [0, R) using the Poisson control points of F at 0 and on [a, a 1 Ra) using the trimming and rendering algorithm presented in Section 4.2.3. Fig. 10 shows an approximation of the function 1/(1 1 x) on the interval [1/2, 2). This function could approximated directlyÐ that is, without trimmingÐonly for values less than 1 (Fig. 7). Moreover, from Theorem 4 the trimming algorithm applies as well for a [ (2R, 0]. Thus, the rendering domain can also be extended to negative values (Fig. 11). More generally, if we consider F as a function of a complex variable, the trimming algorithm extends the domain when the trimming value t a is further than t 0 from the poles. Since a must be in the convergence
4.3.1. Analytic extension Let us ®rst consider the well-known case of an analytic function de®ned by a Taylor series, and then generalize, in the next section, to the Poisson setting. Let F x P k k P F k 0xk =k! for uxu , R: The function Fa x F ax =k! has a radius of convergence of at least R 2 uau, where a [ (2R, R). But Fa may have a radius of convergence greater than R 2 uau. For example, for F x 1= 1 1 x; R 1, but F1=2 x 1= x 1 1=2 1 1 has radius of convergence 3/2 even though r 2 1/2 1/2 [11]. 4.3.2. Application to the Poisson setting For a function F, having a convergent Taylor series is equivalent to having a convergent Poisson series. Analytic extension applies as well in the Poisson representa-
Fig. 11. Approximation of the function F x 1= 1 2 x over the interval [22, 1] 6(solid8 line) by the control polygons of order 5 of the functions 24 2 2 F22 ; F22 ; F22 (dashed lines). The domain has been extended using two successive rounds of trimming: one at 23/4, one at 22.
822
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
Fig. 12. Trimming successively at a1, a2, a3 and a4, should, in theory, enable us to render the function for values greater than the pole. Unfortunately, in practice, this algorithm is numerically unstable.
domain, several rounds of trimming may be necessary, as in the example in Fig. 11. Now, let us consider the case where t a is closer than t 0 to a pole: ² Suppose a0 is a pole of F in R1. Can we approximate F on an interval [a, b] where a . a0? Since any analytic function on R can be extended on C, the pole can be avoided by going into the complex plane, as shown in Fig. 12. Unfortunately we found that, in practice, these successive trimming steps often do not lead to the expected result due to error accumulation. Therefore, we did not manage to approximate successfully on domains beyond a singularity. ² The closest poles to 0 are in C\R1 {x 1 iyuy ± 0 or x , 0}: By successive trimming, we can ®nd a sequence of Poisson series de®ning the function on [0, b] for any b , 1. Using left subdivision, we can then render the function F on [0, b] (Fig. 13). Fig. 14 illustrates this procedure for the function log(2 2 2x 1 x 2), which has poles at 1 1 i and 1 2 i. Approximations to pthis function are computed, ®rst on the interval 0; 2 where the Poisson series
Fig. 13. By trimming successively at a1, a2, ¼, the function can be approximated by our rendering algorithm on [0, b] for any b , 1. This successive trimming can be done as long as there is no pole at any positive real value.
converges, and then, using successive trimming at 1 p p and 2, for values greater than 2.
4.4. Intersection of analytic functions The recursive nature of the BeÂzier subdivision algorithm for polynomials leads to an ef®cient, divide and conquer, intersection algorithm for polynomial curves [9]. From our piecewise linear uniform approximation scheme for Poisson functions, we can generate a similar, divide and conquer, intersection algorithm for analytic curves. To determine at a precision e the intersection points of two analytic functions F and G de®ned on the intervals [a, b] and [c, d], the algorithm proceeds as follows: Find control polygons CF and CG providing an e -approximation of F on [a, b] and of G on [c, d]. intersect_control_polygon(CF, CG) where the function intersect_control_polygon(C1, C2) does: if convex_hull(C1) > convex_hull(C2) 0¤ return. else if C1 is a line segment
Fig. 14. The left ®gure represents the function log(2 2 2x 1 x 2) over the interval [0, 3] (solid line), with psuccessive piecewise linear approximation of orders 1, 2, 4 and 6 (dashed lines). Note that the approximations converge to the curve only on the interval 0; 2: The right ®gure shows p p approximations to the same curve (solid line) of orders 3 and 6 (dashed lines) for values greater than 2, generated by trimming successively at 1 and 2.
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
823
Fig. 15. On the left, intersection between the Archimedean spiral (tcos t, tsin t) and a degree 3 BeÂzier curve on the interval [p, 3p]. On the right, the control polygons successively discarded when seeking the left intersection point are highlighted.
if C2 is a line segment return C1 > C2 else return intersect_control_polygon(C1, ®rst_ half(C2)) < intersect_control_polygon(C1, second_half(C2)) else intersect_control_polygon( ®rst_half(C1), C2) < intersect_control_polygon(second_half(C1), C2) Note that this algorithm can be improved when computing the intersection of an analytic function and a polynomial by using for the polynomial curve the classical BeÂzier approach, subdividing only when we need to split the curve into two parts (Fig. 15). The right side of Fig. 15 shows the discarded control polygons while searching for the left-most intersection in the left part of the ®gure. Fig. 16 illustrates the use of the intersection algorithm to ®nd the roots of a transcendental equation on a
bounded interval. By using the trimming algorithm with n 2 6 at p/2 and then the subdivision approximation algorithm of order 6, the root is found with three significant digits of accuracy. The drawback of this algorithm compared with Newton's method, is that the approximation order has to be decided in advance, whereas the result of Newton's iterative method can be re®ned. Nevertheless, the intersection algorithm based on subdivision has two advantages over Newton's method: ®rst, we choose the interval where we seek the root. The same example, treated by Newton's method with an initial estimate of the root at 0, does not lead to the same result, but rather to a root greater than ®ve. Moreover, a piecewise linear e -approximation of a Poisson function computed once can be used several times to solve different equations involving the same function.
Acknowledgements We wish to thank Mike Wolf for his help with the proof presented in the Appendix. This work has been partially supported by the NSF grant CCR-9971004. Appendix A Here, we shall prove the absolute and uniform convergence in x of the series: X ak k F x k! k$0 Fig. 16. Intersection between the function cos(x) and the rational function x2 2 x 2 8=3 2 1 x on the interval [p/2, p]Ðthat is, the solution of the equation cos x x2 2 x 2 8=3 2 1 x on this interval.
28
on any closed subinterval of 0; R 2 uau, where R is the radius of convergence of F. To establish this result, we need to consider F not as a real analytic function, but as a complex holomorphic function. By hypothesis, the
824
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
magnitude of the pole closest to 0 is R. Thus, for any x in 0; R 2 uau, the series (28) is absolutely convergent. Let us take an arbitrary r in 0; R 2 uau so that r 1 uau , R. Let a r 1 uau 1 R=2 2 r uau 1 R 2 r 1 uau=2; then a . uau, but r 1 a , R. We are going to establish the uniform convergence of the series (28) for any x of magnitude less than or equal to r , that is, for x in D(0, r ), the closed disk of center 0 and radius r . First, note that on the closed disk D(0, r 1 a ) the function admits a maximum m. That is, for any t [ D(0, r 1 a ), uF tu # m. Moreover, from Cauchy's integral formula for the derivatives [13], we have k! Z uF ju uF k xu # dj ; 2p C x;a uj 2 xu k11
References
This inequality establishes the absolute and uniform convergence of the series, relative to x, on any closed disk D(0, r ), where r is chosen arbitrarily in 0; R 2 uau: In particular, this inequality implies the absolute and uniform convergence of the series for any real x in [0, r ].A
[1] Barbour AD, Holst L, Janson S. Poisson approximation. Oxford Science Publications, 1992. [2] Curry HB, Schoenberg IJ. On PoÂlya frequency functions IV: the fundamental spline functions and their limits. J d'Analyse Math 1966;VXII:71±107. [3] Davis PJ. Interpolation and approximation. 2nd ed. New York: Dover Publications Inc, 1975. [4] Farin G. Curves and surfaces for CAGD, a practical guide. 3rd ed. Academic Press, 1992. [5] Goldman R. PoÂlya's urn model and computer aided geometric design. SIAM J Alg Disc Meth 1985;6(1):1±28. [6] Goldman R. Identities for the univariate and bivariate Bernstein basis functions. In: Paeth AW, editor. Graphics gems V, Academic Press Inc, 1995. p. 149±62. [7] Goldman R, Morin G. Poisson approximation. Proceedings of Geometric Modeling and Processing 2000;2000:141±9. [8] Hermann T. On the Szasz±Mirakian operator. Acta Math Acad Sci Hung 1978;32:163±73. [9] Hoschek J, Lasser D. Fundamentals of computer aided geometric design. Wellesley, MA: A.K. Peters, 1993. [10] Lane JM, Riesenfeld RF. A theoretical development for the computer generation and display of piecewise polynomial surfaces. IEEE Trans on Pattern Anal and Mach Intell 1980;PAMI-2:35±46. [11] Lelong-Ferrand J, Arnaudies JM. Cours de matheÂmatiques, tome 2, analyse. Paris: Dunod, 1977. [12] Morin G, Goldman R. A subdivision scheme for Poisson curves and surfaces. Computer Aided Geometric Design 2000;17(9):813±33. [13] Sarason D. Notes on complex function theory. Department of Mathematics, University of California, Berkeley, 1994. [14] Stone MH. A generalized Weierstrauss approximation theorem. In: Buck RC, editor. MAA studies on mathematics, vol. 1, studies in modern analysis, Prentice Hall, 1962. p. 30±87.
GeÂraldine Morin is a PhD student in computer science at Rice University, Houston. She received an engineer degree from the ENSIMAG (Ecole Nationale Superieure d'Informatique et de Mathematiques Appliquees de Grenoble) and a Master Degree of Applied Mathematics from Joseph Fourier University both in 1994, in France. After being a mathematics teacher in high school, she joined Rice in 1997. Her ®elds of interest are geometric modeling, computer graphics and approximation.
Ron Goldman is currently a professor of computer science at Rice University. He received his BS in Mathematics from MIT in 1968 and his MA and PhD in Mathematics from Johns Hopkins University in 1973. He spent 10 years in industry solving problems in computer aided design and geometric modeling: as an R&D Mathematician at Manufacturing Data Systems Incorporated, as a Senior Design Engineer at Ford Motor Company, and as a Principal Controller at Control Data Corporation. His current research interests lie in the mathematical representation, manipulation, and analysis of shape using computers. He is particularly interested in algorithms for polynomial and piecewise polynomial curves and surfaces and their generalizations. His work includes research in computer aided design, solid modeling, computer graphics, and splines.
where C(x, a ) denotes the circle of center x and radius a . Since for any x of magnitude less than or equal to r :C x; a , D 0; r 1 a it follows that uF ju # m: Morek11 over, on C(x, a ): ux 2 ju ak11 : Thus, for any x in D(0, r ): X k ak m 1 X a k : F x # 2p a k$0 a k! k$0
www.elsevier.com/locate/cad
Trimming analytic functions using right sided Poisson subdivision GeÂraldine Morin, Ron Goldman* Computer Science Department, Rice University, 6100 Main Street, Houston, TX 77005-1892, USA Accepted 15 January 2001
Abstract Just as the Bernstein basis and BeÂzier points give us tools to manipulate polynomials, the Poisson basis and control points provide a new framework for investigating analytic functions. Using the Poisson representation, two different subdivision algorithms for analytic functions can be derived as extensions of the BeÂzier polynomial subdivision algorithm. The ®rst procedure generates an approximation algorithm, the second a trimming algorithm. Three applications of these two subdivision algorithms are developed here: an evaluation algorithm, a method for extending the convergence domain of the Poisson representation of an analytic function, and an intersection algorithm for analytic curves. q 2001 Elsevier Science Ltd. All rights reserved. Keywords: Analytic function; Poisson distribution; Subdivision; Approximation; Trimming; Intersection
1. Introduction The blending functions for BeÂzier curves and surfaces are provided by the binomial distribution. Although the Poisson distribution has long been known in probability theory to be a limiting case of the binomial distribution [1], in geometric modeling the idea of considering limiting cases of the polynomial setting only appeared later on [2,5,12]. In approximation theory, however, the Poisson distribution has already been used in the Szasz±Mirakian operator [8]. Analytic functions provide a richer collection of modeling tools than polynomials. In Computer Aided Geometric Design, the Bernstein basis is acknowledged to be more suitable for representing polynomials than the monomial basis. Analogously, we shall ®nd it more natural to use a Poisson series to describe analytic functions, rather than the usual Taylor expansion. The generalization from polynomials in Bernstein±BeÂzier form to analytic functions in the Poisson representation is quite natural and extensive [7], but convergence issues speci®c to the Poisson representation need to be addressed. Moreover, to provide useful algorithms, we need to consider ®nite procedures. Based on the Poisson representation of an analytic function, we shall present two subdivision algorithms for analytic curves. These algorithms are limiting cases of the BeÂzier subdivision algorithm for polynomial curves [4]. The Poisson representation and these associated subdivision algorithms provide new tools for representing and * Corresponding author. E-mail address: [email protected] (R. Goldman).
manipulating analytic functions intuitively and ef®ciently. The ®rst part of this paper establishes the theoretical foundations and addresses convergence issues. The second part concentrates on several applications of these tools: representing an analytic function both inside and outside the convergence domain of its Taylor series, evaluating an analytic function or its derivative at a point, and, ®nally, intersecting analytic curves. The Poisson representation has already been exploited in [12], where the Poisson representation for analytic functions is introduced along with a subdivision algorithm that generates a rendering procedure for analytic functions inside the domain of convergence of their Taylor series. Here, we shall begin by recalling this work, although our de®nition of a Poisson function differs slightly. We shall then present a new subdivision procedure, together with corresponding analysis algorithms. In the last part of this paper, we shall develop some applications of these new procedures.
2. The Poisson representation 2.1. The Poisson basis The Poisson basis is an in®nite collection of functions that we shall denote by bk tk$0 . If Bnk t denotes the k-th Bernstein polynomial of degree n, then: bk t limn!1 Bnk
0010-4485/01/$ - see front matter q 2001 Elsevier Science Ltd. All rights reserved. PII: S 0010-448 5(01)00099-9
t ; n
k $ 0:
1
814
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
Moreover, the convergence is uniform on any bounded interval of [0, 1) (for a proof see the appendix of [12]). The Poisson basis functions can also be expressed explicitly: for any non negative integer k; bk t e2t tk =k!: These basis functions are positive and form a partition of unity on [0, 1), just as the Bernstein basis functions do on the interval [0, 1]. 2.2. Poisson curves We call F a Poisson function or a parametric Poisson curve on the interval [0, R), if for any x0 [ [0, R) the function F is de®ned by a Poisson series in a neighborhood U of P x0: that is, for any t [ U, F t k$0 Pk bk t 2 x0 converges. 1 The vector of coef®cients p (P0, P1, ¼) is the vector of Poisson control points of F at x0.
3. Poisson subdivision 3.1. BeÂzier subdivision Poisson subdivision is an extension of the well-known de Casteljau subdivision algorithm for BeÂzier curves. In this section, we shall recall some classical results for polynomials, in order to be able to refer to them later on. A BeÂzier curve is a polynomial P parameterized over an arbitrary interval [a, b]. If n is the degree of the curve, the control points (P0, P1, ¼, Pn) of the curve are the coef®cients of P in the basis t2a Bnk ; b 2 a 0#k#n that is: P t
2.3. Representation of analytic functions in the Poisson basis A Poisson curve F is de®ned by its Poisson control points at 0 over the interval [0, R) if and only if the Taylor series of F P at 0 converges on this same interval. Indeed, if F k$0 Pk bk on the interval [0, R), then F is P the product of the function t 7 ! e2t and the Taylor series k$0 Pk tk =k! et F t which converges, by hypothesis, over the P interval [0, R). Conversely, it is easy to see that if F t Pk$0 Ak tk =k!; t k then there P exist (P0, P1, ¼) such that e F t Pk t =k!; so F t k$0 Pk bk t: Moreover, there are simple and ef®cient algorithms to convert between the Taylor coef®cients Ak k$0 and the Poisson coef®cients P Pk k$0 [7,12]. The Taylor series F t k$0 Ak tk =k!; converges uniformly to F on P any subinterval [0, b] of [0, R). Similarly, the Poisson series k$0 Pk bk t converges uniformly to F on any interval [0, b] , [0, R). However, if R 1 and limt!1F(t) 0, then, the Poisson series converges uniformly to F on [0, 1) [14]. We shall illustrate the difference of convergence between the Taylor and the Poisson series in the application section. Any function de®ned by a Taylor series over the interval [0, R) is C 1[0, R) and its derivatives are also de®ned by a Taylor series on [0, R) generated by differentiating the original series term by term. Therefore, a Poisson function, de®ned by a Poisson series on [0, R), is also C 1[0, R) and its derivatives are also Poisson functions de®ned by Poisson series generated by differentiating the original series term by term. Moreover, since b 0k t bk21 t 2 bk t; it is easy to show that the vector of Poisson control points of the m-th derivative, F (m), of a Poisson function F with vector of control point p is the m-th discrete difference of p. Thus, if D denotes the difference matrix, i.e. Di,i 21, Di,i11 1, otherwise Di, j 0, then D mp are the control points of F (m). 1 The authors want to thank Jean Gallier for suggesting this de®nition of a Poisson curve.
n X k0
Pk Bnk
t2a : b2a
The polygonal curve de®ned by the control points is called the control polygon of the BeÂzier curve, and is independent of the interval of parameterization [a , b ]. Therefore, a common choice is to take [a , b ] [0, 1]. Here, however, we shall consider the parameterization over the interval [0, n] for two reasons: ®rst, the Poisson basis is a limiting case of the Bernstein basis of degree n over the interval [0, n]; second, we are going to extend the subdivision algorithm in two distinct ways by considering two different approaches to subdivision, but, over the interval [0, 1], these two approaches merge into one. Geometrically, a step of subdivision corresponds to splitting a BeÂzier curve into the union of two BeÂzier curves, as illustrated in Fig. 1; that is, from the control polygon of P, one round of subdivision generates two new control polygons, Qk 0#k#n and Rk 0#k#n ; corresponding to these two new BeÂzier curves. Iterating the subdivision process produces a sequence of control polygons that converges to the corresponding BeÂzier curve [10]. In the polynomial setting, the de Casteljau algorithm is both an evaluation algorithm and a subdivision algorithm (Fig. 2). This dynamic programming procedure offers a nice way to visualize, and an ef®cient way to compute, subdivision. In the following discussion, we shall consider a BeÂzier curve P of degree n parameterized over the interval [0, n]. The subdivision parameter can be de®ned in two different ways: either we can choose r in (0, 1) and subdivide at nr, in
Fig. 1. A BeÂzier curve of degree 3, with its initial control points (on the left) and its control points after one round of subdivision (on the right). The ®rst four of the new control points are generated by the left subdivision algorithm, the last four by the right subdivision algorithm.
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
815
braic identities for the change of basis are: X k a t a n2j t n a 1 12 Bk2j Bnj ; 0 # k # n: Bk n n n n n j0 4 From these equalities, we obtain an expression for the new control points: nX 2k a Rk a Bn2k 0 # k # n: 5 P ; j n k1j j0 Fig. 2. The de Casteljau algorithm for BeÂzier curves. A point emerging from two arrows is an af®ne combination of the points at the base of these arrows, e.g. Q1 (1 2 .)P0 1 .P1.
which case the coef®cients along the arrows in the de Casteljau algorithm would be 1 2 r and r; or we can choose an absolute parameter a, in which case the coef®cients in the de Casteljau algorithm would be 1 2 a/n and a/n. We shall extend both of these approaches to subdivision to the Poisson setting: the ®rst approach generates an algorithm for approximating a Poisson curve; the second approach yields both an evaluation algorithm and a trimming algorithm. In the next paragraph, we present algebraic identities characterizing both of these aspects of subdivision in the polynomial setting. In the ®rst approach, subdivision occurs at the parameter rn where 0 , r , 1. Left-sided subdivision is equivalent to a scaling of the parameter interval: from the control points (Pk) of the curve over the interval [0, n], we compute the control points (Qk) of the restriction of the curve to the interval [0, rn]. The following algebraic identities, from [6], express the Bernstein basis functions over the interval [0, rn] in terms of the original Bernstein basis functionsÐthat is, the basis functions de®ned over the interval [0, n]: Bnk
rt n
n X jk
Bjk rBnj
t ; n
In the next two sections we shall present generalizations of these two subdivision algorithms to the Poisson setting, and in the remainder of the paper we shall develop several applications of these algorithms for analytic functions. 3.2. Approximation: left-sided Poisson subdivision In this section, we brie¯y recall some of the results presented in [12], which we shall use later on to justify certain convergence properties. In the following discussion, F will denote a Poisson curve de®ned by the control points (P0, P1, ¼) (at 0) over the interval [0, R), that is: X Pk bk t; 0 # t , R: F t k$0
The subdivision scheme that we shall describe shortly, induces an approximation algorithm for any such function F. This scheme is an extension of the corresponding polynomial subdivision algorithmÐthat is, the left side of the BeÂzier subdivision algorithm over the interval [0, n] at the parameter rn, where 0 , r , 1. Eqs. (2) and (3) de®ning left subdivision in the polynomial case generalize to the Poisson setting: X j bk rt Bk rbj t; k $ 0; 6 j$0
and 0 # k # n:
2
F rt
X k$0
From these equalities, follow expressions for the points Qk as functions of the points Pk : Qk r
k X j0
Bkj rPj ;
0 # k # n:
3
In the same way, we can derive algebraic identities modeling right-sided subdivision [6]. In this second approach, the subdivision parameter is a number a in [0, n]. One step of this right-sided subdivision is equivalent to computing the control points of the restriction of the curve to be interval [a, n] from the original control points over the interval [0, n]. The corresponding alge-
Qk rbk t with Qk r
k X j0
Pj Bkj r;
k $ 0: 7
Although recursive subdivision converges only for r in (0, 1), equalities (6) and (7) hold for arbitrary real values r as do Eqs. (2) and (3) in the polynomial setting. Moreover, the new control points (Qk(r))k$0 can be generated by applying an extension of the polynomial de Casteljau algorithm to the original control points (Pk)k$0 (Fig. 3). The following theorem, proved in [12], establishes the uniform convergence, as r goes to zero, of the control polygon Q r Q0 r; Q1 r; ¼ to the analytic curve F. Theorem 1. Let L(F, r) denote the piecewise linear function de®ned by the control polygon Q(r), generated after one
816
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
approximation to the Poisson curve de®ned by F over the interval [a, R). These approximations can be re®ned, as desired, using the approximation algorithm presented in Section 3.2.
Fig. 3. The left subdivision algorithm for Poisson curves.
round of subdivision at r, so that L(F, r)(jr) Qj(r), for any integer j, and L(F, r) is linear on [jr, (j 1 1)r], for all j in N. If the sequence of subdivision parameters (ri) converges to zero, then the sequence of piecewise linear functions L(F, ri) converges uniformly to F over any closed subinterval of [0, R). Thus, left-sided subdivision generates a sequence of control polygons converging to a Poisson curve. This convergence provides a rendering algorithm for Poisson curves on any interval [0, b] , [0, R). 3.3. Trimming: right-sided Poisson subdivision 3.3.1. Extension of the polynomial procedure In this section, we shall extend right-sided polynomial subdivision to analytic functions. Equalities (4) and (5) from the polynomial case, generalize readily to the Poisson setting. Indeed, it follows easily from the binomial theorem that: bk a 1 t
k X j0
bk2j tbj a;
so, for a in (2R, R): F a 1 t
X
k$0
Rk abk t where Rk a
8
X j$0
Pk1j bj a;
3.3.2. Approximation of right-sided subdivision Eq. (9) de®nes the control points (Rk(a))k$0 of the Poisson curve F restricted to the interval [a, R), but no algorithm to ef®ciently compute these control points is provided. In fact, it is not possible to compute these points by extending the de Casteljau algorithm (Fig. 2) directly, as in the case of left-sided subdivision. Since in the Poisson setting we start with an in®nite number of control points, the points corresponding to the points on the right lateral edge of the de Casteljau diagram in the polynomial case recede to in®nity. The in®nite sum de®ning these points in Eq. (9) expresses this same phenomenon. To obtain ®nite sums, for every positive integer k we shall de®ne a computable sequence Rnk an$0 so that the Poisson curves Fan de®ned by the control points Rnk ak$0 converge uniformly to the Poisson curve Fa, that is, to the trimmed curve with control points Rk ak$0 : Since the scaled Bernstein basis converges to the Poisson basis, we de®ne: n X a Rnk a Pk1j Bnj ; k $ 0; n $ 0: 10 n j0 Once again, the de Casteljau algorithm allows us to visualize the construction of the points Rnk ak$0 and de®nes a dynamic programming algorithm to generate these points from the initial control points (Pk)k$0 (Fig. 4). We shall use the following lemma to prove that Fan is well de®ned and that the sequences Rnk an$0 and Fan n$0 converge to Rk a and Fa , respectively. P Lemma 2. Let F t Pk bk t be a Poisson curve over the interval P [0, R). Then for any positive integer k, the series i$0 Pi1k bi t de®nes a Poisson curve over
k $ 0:
9 These equations provide a tool for trimming Poisson curves, since the points (Rk(a))k$0 are the Poisson control points of the curve parameterized by the analytic function Fa F(a 1 .). We shall discuss, in detail, in the applications section, the domain where the Poisson series at 0 of Fa convergesÐthat is, where we can approximate Fa by our left-sided subdivision algorithm. In this section, though, we shall consider only the interval 2R 1 uau; R 2 uau; where the function Fa is de®ned by the control points Rk ak$0 : if t [ 2R 1 uau; R 2 uau; then a 1 t [ (2R, R) and, therefore, Eq. (9) converges. Geometrically, the points (Rk(a))k$0 de®ne a control polygonÐthat is, a piecewise linear
Fig. 4. A de Casteljau-like algorithm to compute the points Rnk a:
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
P
the interval [0, R), and X i$0
Pi1k bi
k X j0
Cjk F j ;
11
where Cjk denotes the binomial coef®cient k!= k 2 j!j!: Proof. This result is P proved by induction on k. First, for k 0: by hypothesis i$0 Pi bi F, which is a Poisson function on the interval [0, R). For k 1, observe that: b 0k bk21 2 bk : Using this identity, we can easily deduce that: X Pi11 bi F 1 F 0 : 12 i$0
From Section 1, both F and F 0 are Poisson functions over the P interval [0, R), and, by addition, so is the function i$0 Pi11 bi : Eq. (12) yields, by induction: X i$0
Pi1k bi
k X j0
Cjk F j :
13
Thus, by induction, for any positive k, the function P i$0 Pi11 bi is a Poisson function on the interval [0, R).A Now consider the points (Rkn(a)) and (Rk(a)) all de®ned for a ®xed parameter value a in [0, R). Theorem 3.
limn!1 Rnk a
Rk a
Proof. Note that, if we call and R0(a) the ®rst points of F, then Rnk an$0 and Rk a are the ®rst points of the function with Poisson coef®cients Pm m$k . From Lemma 2, this function is a Poisson function over the interval [0, R). Thus, it is suf®cient to prove that limn!1 Rn0 a R0 a, for a generic Poisson function, that is, for the function F. For a in [0, R), it follows from Theorem 1 with rn a=n that: a lim L F; a F a: n!1 n P But, by construction L F; a=n aP ni0 Pi Bni a=n Rn0 a and by de®nition F a i$0 Pi bi a R0 a. So limn!1 Rn0 a R0 a.A Fan ;
Before proving the convergence of the sequence we P need to verify that the series Rnk abk t converges, that is, the function Fan t is well de®ned: n X n XX a Rk abk t Pk1j Bnj b t n k k$0 k$0 j0
j0
Bnj
X a P b t n k$0 j1k k
But, from Lemma 2, P Pj1k bk t converges on [0, R). Therefore, the in®nite series Rnk abk t Fan t is well de®ned. The following theorem establishes the convergence of the functions Fan ; de®ned by the sequences Rnk ak$0 ; to the function Fa , de®ned by the sequence Rk ak$0 ; where a [ (2R, R). Note that this proof does not use Theorem 3 in which a had to be positive. Theorem 4. As n goes to in®nity, the sequence Fan converges uniformly to Fa on any closed subinterval of 0; R 2 uau. Proof. Let (P0, P1, ¼) be the vector of control points for the Poisson function F. We apply on this initial set of control points one round of the subdivision algorithm at the parameter a/n as illustrated in Fig. 4. Denote by Qn0 a; Qn1 a; ¼ the sequence of points generated by left subdivision: that is, the points that appear on the left lateral edge of the algorithm. In addition to F and Fan ; we now introduce the Poisson functions Qn and Qna de®ned by the control points Qn0 a; Qn1 a; ¼ and Qnn a; Qnn11 a; ¼: From Section 3.2, on left-sided subdivision, we have: a n t ; 14 Q t F n and Qna t
Fan
a t : n
15
Moreover, from the proof of Lemma 2:
(R0n(a))n$0
n X
817
Qna t
n X k0
Ckm Qn k t:
From Eqs. (14) and (16), we get: n X n! ak k a t : F Qna t k k! n k0 n 2 k!n
16
17
Thus, by Eqs. (15) and (17) and using the change of variable x (a/n)t: Fan x
n X k0
n! ak k F x: n 2 k!nk k!
18
But by Taylor's theorem: Fa x F a 1 x
X ak k F x: k! k$0
19
From Eqs. (18) and (19) it appears already that Fan x converges pointwise to Fa(x), but to establish uniform convergence we must show that k n X n! k a F F a 1 x 2 k! n 2 k!nk k0
goes to zero as n goes to in®nity for every x in
818
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
0; r , 0; R 2 uau. For any integer n0, we have:
This last inequality implies:
n0 n X X n! ak ak k k F x x # F a 1 x 2 F F a 1 x 2 k k! k! k0 n 2 k!n k0
n k X 0 n! a k e F x # ; 12 k k! 3 n 2 k!n k0
n k X 0 n! a k F 1 12 x n 2 k!nk k! k0 X k n n! a k : 1 F x kn 11 n 2 k!nk k!
27
for all x [ [0, r ] and n $ n4.Using Eqs. (21) and (24) with m n0 and Eq. (27) for any n . n4, we can bound each term of the left-hand-side of Eq. (20) by e=3: Thus, Fan converges uniformly to Fa on [0, r ].
0
P
20
Let us ®x e . 0. The series F a 1 x k$0 ak =k!F k x converges uniformly relative to x on any closed subinterval [0, r ] of 0; R 2 uau (see the proof in the Appendix). That is, there exists n1 so that, for all m $ n1 and x [ [0, r ]: m X ak k e F x , : 21 F a 1 x 2 k! 3 k0 and, from the absolute and uniform convergence of this same sum, there exists n2 so that for all n $ m $ n2: n X ak k e uF xu , : k! 3 km
22
Moreover, for any positive integers n and k: 0,
n! , 1; n 2 k!nk
23
so from Eqs. (22) and (23), for all m $ n2: n n X X n! ak k ak k e # F uF xu , : x k k! k! 3 km 1 1 n 2 k!n km 1 1
The following corollary extends Theorem 4 to the derivatives of Fa: Corollary 5. For any integer k $ 0, the functions Fan k converge uniformly to Fa k on [0, r ], where r is an arbitrary number in 0; R 2 uau: Proof. This corollary is a consequence of a general result on convergence for analytic functions, that is, if we have a uniform approximation of an analytic function by analytic functions, then the k-th derivative of the approximants also provides a uniform approximation of the k-th derivative of the original function. See Davis [3], p. 112.A In this section, we constructed the Poisson control points of a trimmed Poisson curve, but, since these points were not computable, they were expressed as the limit of order n approximating control points. We then proved that the Poisson functions de®ned by these approximating families of points converge uniformly to the trimmed functions. This convergence property is true for any function de®ned by its Poisson expansion at zero, as well as for all its derivatives.
24 We know (for a detailed proof see the appendix of [12]): n! # k k 2 1 : 1 2 25 k 2n n 2 k!n P Let n0 max(n1, n2). Since ak =k!F k x converges absolutely and uniformly on [0, r ], there exists a constant M . 0 so that, for all x [ [0, r ]: n0 k X a k 26 F x # M: k! k0 Choose n4 so that n4 $ M
n0 n0 2 1 3 : 2 e
From Eqs. (25) and (26), for any x [ [0, r ] and n $ n4: k n0 a X n! e j F max j1 2 x # : j j 0#j#n0 k! 3 n 2 j!n k0
4. Applications Here, we ®rst compare the approximation of an analytic function F by its Taylor polynomial and its Poisson partial sum with the same number of terms. In the remainder of this section we shall use the two subdivision algorithms introduced in Section 3 to develop practical analysis algorithms for Poisson curves and surfaces. First, we present an evaluation procedure for Poisson curves via an approximation algorithm that comes about as a direct consequence of right-sided subdivision. Then, we shall propose a rendering method for analytic curves, on different intervals. In the case of analytic functions with poles, we shall consider more closely than in Section 3 the domain of convergence of a trimmed Poisson function; the domain where a Poisson curve can be approximated using a subdivision algorithm can actually be extended thanks to the trimming algorithm. Finally, we shall describe an intersection algorithm for analytic functions, similar to the intersection algorithm induced by subdivision for polynomial curves.
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
819
Fig. 5. The left ®gures show the functions e2t sint at the top and sin t at the bottom (solid) together with their Taylor (long dashes) and Poisson (short dashes) approximations of order 10. The right ®gures show the difference between the function and the corresponding Taylor (long dashes) and Poisson (short dashes) approximations.
4.1. Comparison of Taylor and Poisson approximations of analytic functions The convergence of a Poisson series to a function F such that limt!1 F t 0 is uniform on [0, 1) (Section 2.3), whereas the convergence of the Taylor series is uniform only on ®nite intervals. As an illustration, Fig. 5 (top) shows the function e 2t sin t and its approximations by Taylor and Poisson sums of order 10Ðthat is, the ®rst 11 terms. The Taylor approximation is a polynomial and diverges rapidly (in this example around 2.5). Moreover, the fact that the Taylor approximation diverges and the Poisson approximation is bounded also gives an advantage to the Poisson series for approximating bounded functions. The bottom part of Fig. 5 illustrates the approximation of the function sin t by the ®rst 11 terms of the Taylor and Poisson series.
and R0n(a) can be computed by applying the de Casteljau algorithm at a/n (Fig. 4). The derivatives of F can also be approximated using the de Casteljau algorithm. From Section 2.3, the coef®cients of the m-th derivative of F can be derived from the coef®cients of F by taking discrete differences (Fig. 6). Then, using the same procedure for F m as for F, we can ®nd a sequence of values converging to F m a for any a in 0; R: 4.2.2. Approximation on a domain [0, b), where b . 0 The approximation algorithm for a Poisson curve on an interval [0, b], where 0 , b , R, is a direct application of
4.2. Approximation of a Poisson function inside the convergence domain of its Poisson series 4.2.1. Evaluation of a Poisson curve and its derivatives at a given parameter value As in the polynomial setting, the de Casteljau algorithm provides an evaluation procedure for analytic curves. At a given parameter a in the domain of convergence, F a R0 a; but, from Theorem 3, we have: lim
n!1
Rn0 a
R0 a;
Fig. 6. From this dynamic programming algorithm, the Poisson control points P m i of the derivative of order m of an analytic function are computed from the Poisson control points Pi of the function.
820
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
Fig. 7. The function F x 1= 1 1 x (solid) and successive approximations using the left-sided Poisson subdivision algorithm of order 1 (clearly piecewise linear curve), 4 (close to the original curve) and 8 (overlapping with the original curve for values less than 1.5 and then diverging).
the left subdivision algorithm and has already been presented and detailed in [12] where, in particular, a justi®cation that only a ®nite number of control points is necessary to compute the approximation is presented. For all the applications given below the subdivision parameter is r 1/2, and the order of the subdivision is the number of iterations of subdivision at r 1/2. But, one step of subdivision at 1/2 i is equivalent to i steps of subdivision at 1/2 [12] and, if the intermediate results are not needed, is more ef®cient. Therefore, we actually compute i rounds of subdivision at r 1/2 by applying one round of subdivision at 1/2 i. Fig. 7 represents approximations of the function 1/(1 1 x) by several control polygons generated from subdivision. This function has a pole at 21; the convergence result for the approximation algorithm holds only in the domain of convergence of the Taylor seriesÐthat is, in this example, for values less than 1. We shall see in Section 4.2 how to approximate this function for values larger than 1. The approximation algorithm extends to surfaces. Fig. 8 illustrates the use of the subdivision algorithm to generate approximations of the sphere. Note that this algorithm is the only known stationary subdivision scheme providing an exact representation of the sphereÐthat is, such that the sphere is the limit of the subdivision process. To compute the initial control
points (Pjk) of the sphere cosucosv; cosusinv; sinu; u [ [0, 2p], v [ [0, 2p], we ®rst generate the Poisson control points ck and (sk) of the trigonometric functions cos t and sin t from their Taylor coef®cients by the algorithm proposed in [12]. Then, the left-sided subdivision algorithm is applied to Pjk j;k$0 cj ck ; cj sk ; sj j;k$0 in both parameter directions. That is, ®rst, one iteration of subdivision is applied to the points Pjk j$0 for each k, to generate the new control points Q~ jk j$0 . Then, subdividing the family of control points Q~ jk k$0 for each j generates the new control points Qij : As with the curves, the subdivision parameter is 1/2 and the order of subdivision is the number of iterations in each direction. The approximation algorithm presented in this paragraph generates arbitrarily good piecewise linear approximations of a Poisson curve on an interval originating at 0 or of a tensor product Poisson surface on a domain originating at (0, 0). In the next section, we shall explore how the trimming algorithm allows us to relax this domain constraint. 4.2.3. Approximation on a domain not containing zero Using both the approximation and trimming algorithms presented in Sections 3.2 and 3.3, we can construct a more general approximation procedure; that is, we can approximate a Poisson function F on a; b; an arbitrary closed subinterval of 2R; R where R is the radius of convergence at zero. From the control points, Pk ; of F we can ®nd a control polygon Rnk a that de®nes a function Fan approximating the trimmed function Fa : We can then apply the approximation algorithm induced by left subdivision to the control points Rnk ak$0 : From Theorem 4, the convergence of the function Fan is uniform on the interval [a, b], and by Theorem 1 the control polygons generated by the leftsided subdivision algorithm also converge uniformly to the function Fan on the interval [a, b]. Therefore, combining these two algorithms, we can construct a sequence of control polygons converging uniformly to the function F on the interval [a, b]. Fig. 9 illustrates this algorithm. This approximation depends on the two subdivision algorithms: for the trimming algorithm, the relevant parameter is the degree4 n, chosen here to be powers of two; in our example F12 is the trimmed function generated for n 2 4 16. For the approximation algorithm, the order,
Fig. 8. A Poisson subdivision scheme converging to the sphere (cosucosv; cosusinv; sinu), u [ [0, 2p], v [ [0, 2p]: control polyhedra of the sphere of orders 1, 2, 3 and 6.
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
Fig. 9. Approximation of the function F(x) cos x (solid line) over the interval [1,4 4] by 8 the control polygons of order 5 of the functions 1 2 F12 ; F12 ; F12 and F12 (dashed lines).
as in the previous section, is the number of iterations at 1/2. In Fig. 9 the order of the approximation algorithm is kept constant in order to observe the convergence of the trimming algorithm. The algorithm presented above is general enough to approximate entire functionsÐthat is, analytic functions having no polesÐon any interval [a, b]. Nevertheless, for an analytic function having a ®nite domain of convergence (2R, R) around zero, we are limited to intervals [a, b] , (2 R, R). In the next section, we shall see how to approximate outside this domain of convergence. In particular, this new method will solve the divergence problem illustrated in the example in Fig. 7. 4.3. Approximation of a Poisson function outside the initial convergence domain Here, we address the approximation of analytic functions with poles. So far, we are able to approximate functions with poles on any closed subinterval of their domain of convergence around zero. Below, trimming is used to extend the domain where a Poisson function is de®ned by a Poisson series. The rendering process induced by left subdivision then applies on this extended domain.
821
Fig. 10. Approximation of the function F x 1= 1 1 x over the interval [1/2, 2) 2(solid4 line) by the control polygons of order 5 of the functions 20 2 2 F1=2 ; F1=2 ; F1=2 (dashed lines). Note that the singularity appearing in Fig. 7 for x . 1 has been eliminated.
tion; that is, the radius of convergence Ra of the Poisson series de®ning Fa may be greater than R 2 uau. Moreover, a Poisson series is associated with a rendering algorithm on the interval [0, Ra). Thus, the function F can be rendered on [0, R) using the Poisson control points of F at 0 and on [a, a 1 Ra) using the trimming and rendering algorithm presented in Section 4.2.3. Fig. 10 shows an approximation of the function 1/(1 1 x) on the interval [1/2, 2). This function could approximated directlyÐ that is, without trimmingÐonly for values less than 1 (Fig. 7). Moreover, from Theorem 4 the trimming algorithm applies as well for a [ (2R, 0]. Thus, the rendering domain can also be extended to negative values (Fig. 11). More generally, if we consider F as a function of a complex variable, the trimming algorithm extends the domain when the trimming value t a is further than t 0 from the poles. Since a must be in the convergence
4.3.1. Analytic extension Let us ®rst consider the well-known case of an analytic function de®ned by a Taylor series, and then generalize, in the next section, to the Poisson setting. Let F x P k k P F k 0xk =k! for uxu , R: The function Fa x F ax =k! has a radius of convergence of at least R 2 uau, where a [ (2R, R). But Fa may have a radius of convergence greater than R 2 uau. For example, for F x 1= 1 1 x; R 1, but F1=2 x 1= x 1 1=2 1 1 has radius of convergence 3/2 even though r 2 1/2 1/2 [11]. 4.3.2. Application to the Poisson setting For a function F, having a convergent Taylor series is equivalent to having a convergent Poisson series. Analytic extension applies as well in the Poisson representa-
Fig. 11. Approximation of the function F x 1= 1 2 x over the interval [22, 1] 6(solid8 line) by the control polygons of order 5 of the functions 24 2 2 F22 ; F22 ; F22 (dashed lines). The domain has been extended using two successive rounds of trimming: one at 23/4, one at 22.
822
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
Fig. 12. Trimming successively at a1, a2, a3 and a4, should, in theory, enable us to render the function for values greater than the pole. Unfortunately, in practice, this algorithm is numerically unstable.
domain, several rounds of trimming may be necessary, as in the example in Fig. 11. Now, let us consider the case where t a is closer than t 0 to a pole: ² Suppose a0 is a pole of F in R1. Can we approximate F on an interval [a, b] where a . a0? Since any analytic function on R can be extended on C, the pole can be avoided by going into the complex plane, as shown in Fig. 12. Unfortunately we found that, in practice, these successive trimming steps often do not lead to the expected result due to error accumulation. Therefore, we did not manage to approximate successfully on domains beyond a singularity. ² The closest poles to 0 are in C\R1 {x 1 iyuy ± 0 or x , 0}: By successive trimming, we can ®nd a sequence of Poisson series de®ning the function on [0, b] for any b , 1. Using left subdivision, we can then render the function F on [0, b] (Fig. 13). Fig. 14 illustrates this procedure for the function log(2 2 2x 1 x 2), which has poles at 1 1 i and 1 2 i. Approximations to pthis function are computed, ®rst on the interval 0; 2 where the Poisson series
Fig. 13. By trimming successively at a1, a2, ¼, the function can be approximated by our rendering algorithm on [0, b] for any b , 1. This successive trimming can be done as long as there is no pole at any positive real value.
converges, and then, using successive trimming at 1 p p and 2, for values greater than 2.
4.4. Intersection of analytic functions The recursive nature of the BeÂzier subdivision algorithm for polynomials leads to an ef®cient, divide and conquer, intersection algorithm for polynomial curves [9]. From our piecewise linear uniform approximation scheme for Poisson functions, we can generate a similar, divide and conquer, intersection algorithm for analytic curves. To determine at a precision e the intersection points of two analytic functions F and G de®ned on the intervals [a, b] and [c, d], the algorithm proceeds as follows: Find control polygons CF and CG providing an e -approximation of F on [a, b] and of G on [c, d]. intersect_control_polygon(CF, CG) where the function intersect_control_polygon(C1, C2) does: if convex_hull(C1) > convex_hull(C2) 0¤ return. else if C1 is a line segment
Fig. 14. The left ®gure represents the function log(2 2 2x 1 x 2) over the interval [0, 3] (solid line), with psuccessive piecewise linear approximation of orders 1, 2, 4 and 6 (dashed lines). Note that the approximations converge to the curve only on the interval 0; 2: The right ®gure shows p p approximations to the same curve (solid line) of orders 3 and 6 (dashed lines) for values greater than 2, generated by trimming successively at 1 and 2.
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
823
Fig. 15. On the left, intersection between the Archimedean spiral (tcos t, tsin t) and a degree 3 BeÂzier curve on the interval [p, 3p]. On the right, the control polygons successively discarded when seeking the left intersection point are highlighted.
if C2 is a line segment return C1 > C2 else return intersect_control_polygon(C1, ®rst_ half(C2)) < intersect_control_polygon(C1, second_half(C2)) else intersect_control_polygon( ®rst_half(C1), C2) < intersect_control_polygon(second_half(C1), C2) Note that this algorithm can be improved when computing the intersection of an analytic function and a polynomial by using for the polynomial curve the classical BeÂzier approach, subdividing only when we need to split the curve into two parts (Fig. 15). The right side of Fig. 15 shows the discarded control polygons while searching for the left-most intersection in the left part of the ®gure. Fig. 16 illustrates the use of the intersection algorithm to ®nd the roots of a transcendental equation on a
bounded interval. By using the trimming algorithm with n 2 6 at p/2 and then the subdivision approximation algorithm of order 6, the root is found with three significant digits of accuracy. The drawback of this algorithm compared with Newton's method, is that the approximation order has to be decided in advance, whereas the result of Newton's iterative method can be re®ned. Nevertheless, the intersection algorithm based on subdivision has two advantages over Newton's method: ®rst, we choose the interval where we seek the root. The same example, treated by Newton's method with an initial estimate of the root at 0, does not lead to the same result, but rather to a root greater than ®ve. Moreover, a piecewise linear e -approximation of a Poisson function computed once can be used several times to solve different equations involving the same function.
Acknowledgements We wish to thank Mike Wolf for his help with the proof presented in the Appendix. This work has been partially supported by the NSF grant CCR-9971004. Appendix A Here, we shall prove the absolute and uniform convergence in x of the series: X ak k F x k! k$0 Fig. 16. Intersection between the function cos(x) and the rational function x2 2 x 2 8=3 2 1 x on the interval [p/2, p]Ðthat is, the solution of the equation cos x x2 2 x 2 8=3 2 1 x on this interval.
28
on any closed subinterval of 0; R 2 uau, where R is the radius of convergence of F. To establish this result, we need to consider F not as a real analytic function, but as a complex holomorphic function. By hypothesis, the
824
G. Morin, R. Goldman / Computer-Aided Design 33 (2001) 813±824
magnitude of the pole closest to 0 is R. Thus, for any x in 0; R 2 uau, the series (28) is absolutely convergent. Let us take an arbitrary r in 0; R 2 uau so that r 1 uau , R. Let a r 1 uau 1 R=2 2 r uau 1 R 2 r 1 uau=2; then a . uau, but r 1 a , R. We are going to establish the uniform convergence of the series (28) for any x of magnitude less than or equal to r , that is, for x in D(0, r ), the closed disk of center 0 and radius r . First, note that on the closed disk D(0, r 1 a ) the function admits a maximum m. That is, for any t [ D(0, r 1 a ), uF tu # m. Moreover, from Cauchy's integral formula for the derivatives [13], we have k! Z uF ju uF k xu # dj ; 2p C x;a uj 2 xu k11
References
This inequality establishes the absolute and uniform convergence of the series, relative to x, on any closed disk D(0, r ), where r is chosen arbitrarily in 0; R 2 uau: In particular, this inequality implies the absolute and uniform convergence of the series for any real x in [0, r ].A
[1] Barbour AD, Holst L, Janson S. Poisson approximation. Oxford Science Publications, 1992. [2] Curry HB, Schoenberg IJ. On PoÂlya frequency functions IV: the fundamental spline functions and their limits. J d'Analyse Math 1966;VXII:71±107. [3] Davis PJ. Interpolation and approximation. 2nd ed. New York: Dover Publications Inc, 1975. [4] Farin G. Curves and surfaces for CAGD, a practical guide. 3rd ed. Academic Press, 1992. [5] Goldman R. PoÂlya's urn model and computer aided geometric design. SIAM J Alg Disc Meth 1985;6(1):1±28. [6] Goldman R. Identities for the univariate and bivariate Bernstein basis functions. In: Paeth AW, editor. Graphics gems V, Academic Press Inc, 1995. p. 149±62. [7] Goldman R, Morin G. Poisson approximation. Proceedings of Geometric Modeling and Processing 2000;2000:141±9. [8] Hermann T. On the Szasz±Mirakian operator. Acta Math Acad Sci Hung 1978;32:163±73. [9] Hoschek J, Lasser D. Fundamentals of computer aided geometric design. Wellesley, MA: A.K. Peters, 1993. [10] Lane JM, Riesenfeld RF. A theoretical development for the computer generation and display of piecewise polynomial surfaces. IEEE Trans on Pattern Anal and Mach Intell 1980;PAMI-2:35±46. [11] Lelong-Ferrand J, Arnaudies JM. Cours de matheÂmatiques, tome 2, analyse. Paris: Dunod, 1977. [12] Morin G, Goldman R. A subdivision scheme for Poisson curves and surfaces. Computer Aided Geometric Design 2000;17(9):813±33. [13] Sarason D. Notes on complex function theory. Department of Mathematics, University of California, Berkeley, 1994. [14] Stone MH. A generalized Weierstrauss approximation theorem. In: Buck RC, editor. MAA studies on mathematics, vol. 1, studies in modern analysis, Prentice Hall, 1962. p. 30±87.
GeÂraldine Morin is a PhD student in computer science at Rice University, Houston. She received an engineer degree from the ENSIMAG (Ecole Nationale Superieure d'Informatique et de Mathematiques Appliquees de Grenoble) and a Master Degree of Applied Mathematics from Joseph Fourier University both in 1994, in France. After being a mathematics teacher in high school, she joined Rice in 1997. Her ®elds of interest are geometric modeling, computer graphics and approximation.
Ron Goldman is currently a professor of computer science at Rice University. He received his BS in Mathematics from MIT in 1968 and his MA and PhD in Mathematics from Johns Hopkins University in 1973. He spent 10 years in industry solving problems in computer aided design and geometric modeling: as an R&D Mathematician at Manufacturing Data Systems Incorporated, as a Senior Design Engineer at Ford Motor Company, and as a Principal Controller at Control Data Corporation. His current research interests lie in the mathematical representation, manipulation, and analysis of shape using computers. He is particularly interested in algorithms for polynomial and piecewise polynomial curves and surfaces and their generalizations. His work includes research in computer aided design, solid modeling, computer graphics, and splines.
where C(x, a ) denotes the circle of center x and radius a . Since for any x of magnitude less than or equal to r :C x; a , D 0; r 1 a it follows that uF ju # m: Morek11 over, on C(x, a ): ux 2 ju ak11 : Thus, for any x in D(0, r ): X k ak m 1 X a k : F x # 2p a k$0 a k! k$0
