Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Implementation aspects of Bézier surfaces to model complex geometric forms
(USC Thesis Other)
Implementation aspects of Bézier surfaces to model complex geometric forms
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
IM PLEM ENTATION ASPECTS O F BEZIER SURFACES TO MODEL COM PLEX GEOM ETRIC FORMS by Christian Molt A Thesis Presented to the FACULTY O F TH E GRADUATE SCHOOL UNIVERSITY O F SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirem ents for the Degree M ASTER OF SCIENCE (Applied M athem atics) December 1994 UNIVERSITY O F S O U T H E R N CA LIFO RN IA T H E GRADUATE SC H O O L U N IV ER SITY PA RK LO S A N G ELES, C A L IFO R N IA 0 0 0 0 7 This thesis, written by ...Chris tian__Molt _ ...................... under the direction of A .i? .....T he sis Com m ittee, and approved by all its members, has been p re sented to and accepted by the D ean of T he Graduate School, in partial fulfillm ent of the requirements for the degree of Master of Science (Applied Mathematics) Dean THESIS) COMMITTEE q , a Chairman J&±zh!Lr... Dedication To my parents who made it all possible Acknowledgements I would like to express my sincerest gratitude to my advisor Professor Wlodek Proskurowski. He provided the careful guidance needed throughout the process of writing this thesis. My work benefited significantly from his valuable advice and his encouragement. I would also like to thank Professor Edward K. Blum and Professor Julian A. Domaradzki for their support on my thesis committee. The Rocketdyne Division of Rockwell International Corporation provided the data for the numerical experiments. I want to express my appreciation for their support. A very special thanks to Julia Gromadecki and my parents for their moral sup port. They always kept me up and never let me down although they were far away. I would also like to thank Marc Hombeck and my housemates Stephan Beifuss, Mar tin Eberhardt and Andreas Fugmann. They kept me from working too hard most of the time and helped me to maintain my motivation. Contents D edication ii A cknow ledgem ents iii List O f Tables vi List O f Figures vii 1 Introduction 1 2 B ezier C urves and their P roperties 4 2.1 Example of a Bezier Curve for the Quadratic C a se .................................. 5 2.2 Basic P roperties................................................................................................. 9 2.2.1 Bm lies within the convex hull of its control p o in ts..................... 10 2.2.2 Bm is s y m m e tric .................................................. 10 2.2.3 B m interpolates its first and last control point . . . . . . . . . 11 2.2.4 B m is tangent to its control polygon at the e n d p o in ts............... 12 2.2.5 Bm has quasi local c o n tro l........................................ 14 2.2.6 Bm can be blended by two Bezier curves of degree (m — 1) . . 15 2.2.7 Bm is invariant under affine mappings .................................. 16 2.2.8 Bm has the variation diminuishing p r o p e r t y ............................... 17 2.3 Param etric Representation of Bezier C u r v e s ........................... IS 2.4 M atrix Representation of Bezier C u rv es...................................................... 19 2.5 Bezier S u rfa c e s ................. 21 3 T h e Subdivision A lgorithm and its Im plem entation 24 3.1 The Subdivision T h e o r e m ..................................... ....................................... 25 3.2 The Generation of a Bezier Curve B m ......................................................... 29 3.2.1 The subdivision a lg o rith m ..................................................... 29 3.2.2 Generation of a curve Bm ................................................................ 31 3.2.3 Implementation in C .............................. 35 3.3 The Generation of a Bezier Surface £?m,n .................................................. 36 3.3.1 The subdivision a lg o rith m ................................................................ 36 iv 3.3.2 Generation of a surface Bm > n .................................................. 39 3.3.3 Implementation in C ......................................................................... 43 3.4 The Subdivision Algorithm in Matrix R ep resen tatio n ..............................44 3.4.1 Formulation for Bezier c u r v e s ........................................................ 45 3.4.2 Formulation for Bezier s u r f a c e s ................................ 47 4 Test R esu lts on Im peller B ladings from th e Space Shuttle 52 4.1 General Performance Evaluation of the Subdivision Algorithm . . . . 56 4.2 Testing the Stopping C r i t e r i a ................................................................... 60 5 C onclusions 67 R eference List 71 A p p en d ix A Program Codes of the Implemented A lg o rith m s............................................... 73 A .l Program for the Generation of a Bezier C u r v e ...................................... 73 A.2 Program for the Generation of a Bezier S u r f a c e ................................... 86 A.3 Computation of the Machine P re c is io n ............................ 107 A p pend ix B Test D a t a ............................................................................................. 108 B .l Plots of the Control M e s h e s ..........................................................................108 B.2 Comprehensive Test Results for the Stopping C rite ria ............................. I l l B.3 Samples of Numerical Test Results .............................................................116 v L ist O f T ables 4.1 Performance analysis for mesh 3 .................................................................. 57 4.2 Computational cost of the subdivision algorithm for surfaces . . . . . 59 4.3 Comparison of stopping criteria for mesh 1 and tolerance 0.01 . . . . 61 4.4 Comparison of stopping criteria for mesh 2 and tolerance 0.01 . . . . 61 4.5 Comparison of stopping criteria for mesh 3 and tolerance 0.01 . . . . 62 4.6 Comparison of stopping criteria for mesh 1 and tolerance 0.001 . . . . 64 4.7 Comparison of stopping criteria for mesh 2 and tolerance 0.001 . . . . 65 4.8 Comparison of stopping criteria for mesh 3 and tolerance 0.001 . . . . 65 vi L ist O f F igu res 2.1 Bernstein polynomials of degree fiv e ............................................................ 5 2.2 Quadratic Bezier curve and polar f o r m ...................................................... 6 2.3 Combination of two Bezier curves B4 and B 7 ......................................... 13 2.4 Example of a Bezier surface with control m e s h ........................................ 22 3.1 Subdivision for a quartic Bezier c u rv e .................... 28 3.2 Control mesh of a Bezier surface................................................................... 37 4.1 Blade created by control meshes 1 and 4 .................................................. 53 4.2 Blade created by control meshes 2 and 5 .................................................. 54 4.3 Blade created by control meshes 3 and 6 . . . ................................ . 55 4.4 Comparison of runtime versus output for control mesh 3 58 vii C h a p ter 1 Introduction In the late 1950’s and 1960’s, companies like Boeing in the aircraft industry or Citroen and Renault in the car industry started experimenting with computerized methods in the design process of their products. Those companies had realized that their design process was inflexible, lengthy, costly and error prone. Let us briefly review the design process in the car industry during the 1960’s. Af ter deciding on the main characteristics of a new automobile, the styling departm ent built various models until the shape of the new car was accepted by management. Then, the model was surveyed and translated into drawings. These drawings in cluded all parts of the new car and every piece of equipm ent (like shock absorbers, windscreen, transmission, ...). The latter process was in particular likely to be in accurate. After the generation of the drawings, a m aster of all parts was built, the tools were designed, and the manufacturing process was prepared. Furthermore, 1 most parts had to be fitted manually due to the limited accuracy of the production machines (see Bezier in [3] and [4] for further description). Towards the end of the 1960’ s, computers were sophisticated enough to be used in many areas within companies. Firms like Citroen and Renault started to experiment with computers in the design process. Their goals were to make the design process faster, more flexible, more consistent, but less costly and less error prone. Besides these business management requirements their experiments faced further challenges. First, if computer systems were to be used in the design process, then designers and stylists needed to understand the system easily. These people had no advanced mathematical background, but they had good geometric knowledge about the shapes and forms of surfaces. Second, the data needed to be complete, distortion free and indisputable (see [3]). The latter two requirements were obviously the fundamentals for the success of a computerized design system. In addition, having a numerical representation of the object allowed it various groups to use the data simultaneously. To enable this, the data needed to be exchanged quickly. Hence, they wanted to have a small set of points “controlling” the object instead of the object itself (in fact, 1000 control points can easily represent objects of 10 million data points if a certain accuracy is to be achieved). This speeded up the transaction process and reduced the electronic data exchange. Independent from each other, P. de Casteljau at Citroen and P.E. Bezier at Re nault developed such systems for Computer Aided Design (CAD) and Computer 2 Aided Manufacturing (CAM) in the 1960’s. Yet, most of the credit was given to Bezier due to the fact that de Casteljau was hindered from publishing his results until 1974. Bezier developed a mathematical theory, the so called Bezier curves, to generate curves and surfaces numerically. His theory is based on Bernstein polyno mials (see for example [4]). In the next chapter, we will give an introduction to the theory of Bezier curves and surfaces. We will focus especially on those properties of Bezier curves indicating that this concept is rich and flexible enough to have almost all desired properties for the design process. We will also demonstrate that Bezier curves are highly practicable and applicable in the design process because they have an easy underlying concept. After the discussion of those properties, we will consider an efficient algorithm and its implementation in the C language for rendering curves and surfaces. The basic concept of this algorithm dates back to de Casteljau, and, therefore, it is often referred to as the de Casteljau algorithm. Finally, we will discuss the crucial features of the implemented algorithm. For the testing of the algorithm, data of the bladings for the first stage High Pressure Fuel Turbo Pump in the Space Shuttle Main Engine is used.1 This testing includes the performance of the algorithm in general, but it also concentrates on a specific feature of the algorithm: the stopping criterion for the curve and surface approximation. Hhe data comes courtesy from the Rocketdyne Division of Rockwell International Corporation 3 C h a p te r 2 Bezier Curves and their Properties Let a set of (m -f- 1) points be given, where m > 1. Then the Bezier curve of degree m is defined as m Pk , u € [0, l], 1 (2.1) k=0 where bm^{u) = — u)m~k, k = 0 ,1 ,..., m , is the &-th Bernstein basis polynom ial2 of degree m, and Pf. € 7/2, k = 0,1 , ...,m , are the control points of B m(u). See Figure 2.1 for an example of the Bernstein polynomials of degree five. Before we discuss the properties of Bezier curves and surfaces, a geom etric ap proach to Bezier curves using the concept of polar forms is presented. This is done for the quadratic Bezier curve as described by Seidel in [19]. 1w.l.o.g. assume u G [0,1], otherwise use the substitution u — v € [a> ^ ]> a,b £ Ml, a < b 2(™) denotes the binomial coefficient, i.e. (™) = , 0 < k < m 4 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.2 0.4 0.6 0.8 Figure 2.1: Bernstein polynomials of degree five 2.1 Example of a Bezier Curve for the Quadratic Case Consider a quadratic polynomial curve F : IR -4 JR2. The polar form / : IFF -4 IFF of F is given by 1. the intersection of the tangentline T\ at F(u) with the tangentline T2 at F(v), if 2. the point F(u), if u = v 0.9 0.8 f(u, v) 0.7 T2 0.6 T1 0.5 0.4 f(v, v) 0.3 0.2 f(u, u) 0.4 0.2 0.8 0.6 Figure 2.2: Q uadratic Bezier curve and polar form Formally, f :JR2 IR2, (u,u) — 5 - < 7i H T2 if u 7 ^ v F (u ) if u = v (2.2) Figure 2.1 illustrates the definition. Seidel [19] shows th at a quadratic polynomial F can be represented algebraically by /( u ,u ) — F(u) + li^ LF > (u) using a Taylor expansion for F in u. The polar form / has the following obvious properties: it is sym m etric, i.e. /( u ,u ) = f ( v , u ) 6 • it is multiaffine, i.e. /(w, ett?1-t-(l~a)u2) = ocf(u,vi)+ (l—ct)f(u,v2) , [0, 1] (and analogously for u) • it is diagonal, i.e. F (u ) = /(u , u) One can show that these three properties actually characterize a polar form. The importance of polar forms becomes clear from the “Blossoming Principle”: Let F : IR — ) • IB? be a quadratic polynomial. Then, there exists a unique polar form / : IB? — > ■ IR? which satisfies f(u ,u ) = F(u). Using the polar form, we want to characterize a quadratic Bezier curve B 2(u) on [0,1]. Setting u = u + (l — u)0 and using the fact that / is affine in both coordinates, we get F(u) = f(u ,n ) = (1 - u )2/(0 ,0 ) + 2u(l - « ) /( 0 ,1) + u 2/ ( l , l ) where the control points are given by P0 = / ( l , 1), Pt = / ( 0, 1), P2 = /(0 ,0 ). Thus, we showed how the three control points Pk, k = 0,1, 2, determine the Bezier curve. From the above we can deduce the most im portant properties of F (see Fig ure 2.1): the Bezier curve B2 interpolates the two endpoints P0 and P2', it lies in the convex hull of its control points and-by construction-it is tangent to the control polygon PqP\P2 at its endpoints. That means, if we want to connect two Bezier 7 curves B\ and B% continuously, we need the last control point of B\ to be equal to the first control point of £?|. The polar form / of B 2 is also useful to derive the algorithm for computing Bm(u) for any u € [0,1]. Note that B 2(u) = /(u , u) and u = u + (1 — «)0. Then we get the following algorithm using the symmetry and multiaffinity of /: 1ststep : /(« , u) = u f(u , 1) + (1 - t i) /( 0 , u) 2ndstep : /(«, 1) = u/(l, 1) + (1 - u)/(0,1) (2.3) / ( 0, u) = u /( 0, 1) + (1 — w )/(0, 0) Both steps are based on the same formula. In the first step we combine the two linear functions /(« , 1) and f ( 0,u) to determine /(n ,n ), while in the second step we blend the linear functions by the constant functions / ( 0, 0) ,/ ( 0, 1) and / ( 1, 1). Thus, we can compute the Bezier curve Bm by two Bezier curves Bm-1 , m = 1,2. In addition, we see that B 2(u) can be computed using just the control points. We shall find out that this scheme can be generalized and results in an algorithm for computing an arbitrary point of Bm. Polar forms are a fruitful method for studying curves and surfaces. They were introduced by de Casteljau (see [5] for a comprehensive summary) and they are still heavily used in Computer Aided Geometric Design. Seidel [20] also generalized the concept of polar forms to prove the basic properties for Bezier curves and surfaces of arbitrary degree m. S We will not continue with polar forms, but will work with a more traditional ap proach dating back to Bezier and Forrest (see [4] and [8]). This approach is more suitable for the description of an algorithm for computing points on Bezier curves and surfaces. 2.2 Basic Properties A Bezier curve Bm has the following properties: • Bm lies within the convex hull of its control points • Bm is symmetric • Bm interpolates its first and last control point • Bm is tangent to its control polygon at the endpoints • Bm has quasi local control • Bm can be blended by two Bezier curves of degree (m-1) • Bm is invariant under affine mappings • Bm is variation diminishing Each property by itself does not seem to be very exciting. Putting all those together, however, results in a very powerful tool for CAD applications: Bezier curves have almost all properties that are needed in the design process. 9 Let the Bezier curve Bm be defined as in (2.1) and let PQ Pv,..Pm denote the control polygon of Bm, Now we will prove the properties mentioned above. 2.2.1 Bm lies within the convex hull of its control points Consider YX=o QkPk, where a k = (jf)uk(l — u)m~h. Then o& > 0, V fc = 0,1,..., m , and by the the binomial theorem m m ( rn\ 11 = S 7 , ^ '( i “ “ )m“* = (« + (1 - u))m = 1, A -=0 k=0 \ K/ i.e. B m lies in the convex hull of the control points Po, P i ,..., Pm. The convex hull property is important, because it assures that the shape of the curve is predictable and reasonably nice. In many cases the convex hull property gives a short and efficient test whether two Bezier curves intersect: the designer can be sure that two curves do not intersect if their respective convex hulls do not intersect. 2.2.2 Bm is symmetric We want to show that Bm(u,[P0,P u ...,Pm}) = Bm( l - u , [Pm,P m_ !,.„ ,P 0]) VuG [0,1] (2.4) 10 Proof: B m(u,\pa,P i pm}) = f ; ( ™ ) u i ( i - u r - ‘ p,, = Bm(u, [Pm,Pm_i,...,P0]) Sym m etry guarantees that the control points Po, Pi, Pm result in the same curve as the control points Pm, Pm- i , ..., P0. This is a highly desirable property in Com puter Aided Geometric Design (CAGD), because it supports the intuition of a designer. We will apply this property to save computational work. 2.2.3 Bm interpolates its first and last control point s m(0) = £ ? =0 0*(1 - o)m~kPk = Po and B m{l) = £?=„ l fc (l - l ) m-*P A = Pm, hence, the two endpoints are interpolated. In general, these are the only control points interpolated by the Bezier curve. Due to this property, two Bezier curves B m(u, [Po, Pi, •••, Pm]) and B n(u, [Q0, Q i,..., Qn]) of degree m and n respectively, can be connected continuously by choosing Pm = Qq. Chasen [6] describes a method how one can achieve that a Bezier curve interpolates a specific point on B m that is not an endpoint. 11 2.2.4 Bm is tangent to its control polygon at the endpoints The derivative of Bm is given by ^ B m(u,[P0,P u .., P,n]) — rn(Bm-i(u, [P i,P 2,...,P m]) - Bm-i(u,[P 0, Pu Pm_i])) (2.5) Proof: — B m(u,[P0,P i,...,P m]) = E “ u)m~kPk - g (m “ fc)(m - * ) U‘ (1 - ^ m- k- l Pk d = ‘ m P m_ i(u ,[P i,P 2,...,P m]) - [Po,Pl, -..,Pm-l3) This gives us an algorithm for the computation of the derivative by blending two curves of degree (m -1). Equation (2.5) can be generalized to obtain higher derivatives (see for example [7]). It follows from (2.5) that B,'„(0) = 1 ('" ) P, - m f ™ ) Po = m(P , - Po) , and B L ( i ) = m ( m ) p m - ( 171 ) Pm- i = m ( P m — P m - i) ■ \ m J \ m — l J 12 0.9 V 0.8 0.7 0.6 / / 0.5 / 0.4 0.3 PI 0.2 PO Q2 Q 1 0.2 0.4 0.6 0.8 Figure 2.3: Combination of two Bezier curves and By Interpreting Pk, k = 0,1 ,...,m , as vectors in the plane /R2, it is clear that Bm is tangent to its control polygon at the endpoints. Hence, we are able to obtain Wl con tinuity between two Bezier curves Bm(u,[P0, Pl, P m]) and Bn(u, [Q0, Q j,..., Q,\) if the following holds (see Figure 2.3): 1. Pm = Qo 2. Pm-iP m and Q0Qi are collinear and \Pm- tPm\ = \QoQ\\. 13 If the designer knows these interdependencies, he can connect various Bezier curves smoothly. (Z * 1 continuity fixes two control points of a new Bezier curve. By computing higher derivatives of Bm we see that achieving Wk continuity requires (k+1) control points of of the new Bezier curve to be fixed. Consequently, higher degree continuity results in a lower degree of freedom for choosing control points of a new Bezier curve. 2.2.5 Bm has quasi local control We know from Definition (2.1) that the Bernstein polynomial associated with the control point Pf. is defined over the whole interval [0,1]. Thus, altering one of the control points affects the shape of the whole curve and not only parts of it, i.e. Bezier curves do not enjoy local control. Yet, a change of a control point leads to a modification of the curve largely in the region of that particular control point. We can deduce from (2.5) that — u )m~k~l (k " m u ) (2 -6 ) The Bernstein polynomial bm < k attains its (only) maximum at a = ^ . This means that a move of the control point Pf. has its largest effect on the curve around the param eter value u — Therefore, the designer can guess from past experience how a modification of one of the control points effects the shape of the entire curve. A way to avoid the problem attached with global control is to partition the curve. Then, the change of one part does not affect other parts of the partition. 14 2.2.6 B m can be blended by two Bezier curves of degree (m — 1) Bm can be w ritten as a convex combination of two Bezier curves of degree (m — 1): +uBm- 1(u, [PUP2,,.., Pm] ) , u e [0,1] (2.7) Proof: Bm(u, [Po, Pi,---, Pm]) = E (X)«‘d = E (I)^d - u)"~kP k + (1 - u )”A = E ((m ; ') + (::()) « * (! - «)-‘A + ( 1 - «ra = E ( m f c + E ( " * * 1V + , u - « r - i - 1n + i + ( i - « r f t = E (m * + E (m * 1 )«M - ld - lft+ . = (i - “ ) E ( " 7 1) « f c (i - « ) “ -■-*f t + « E ( m * ^ “ ‘ C 1 - d = ’ (1 - u)Bm-i(u, [P0, P i,.. .,P ro.i] ) + Mflm_i(u, [Pi, P2, Pm]) 15 Equation (2.7) allows us to compute a Bezier curve of degree m by blending two curves of degree (m — 1). We will see that, if we continue doing so, Bm can be written as a sum of its control points. Equation (2.7) is a generalization of (2.3) in Section 2.1. Furthermore, it is the basis for the subdivision algorithm3 described in Section 3.2.1. The above formula is based on degree reduction. Besides degree reduction, de gree raising is sometimes also desirable. In [8], Forrest gives an algorithm for degree raising of a Bezier curve by defining new control points using Fo, Pi, •••, Pm- De gree raising is required if two curves need to be of the same degree for a certain application. 2.2.7 B m is invariant under affine mappings Let < jj : IR , — > 1R be an affine mapping. Then we want to show that [P0,F i,...,P m])) = Bm(u, [^(J P0),^ (F i),...,^ (P m)]) (2.8) Proof: by induction on m (i) Let m=l, then [P0, />,])) = ^ ((o)“ °( 1 “ U)P° + ( l ) “ (1 “ “ )0Pl) = <6((1 - u)P0 + u l\) “= (1 - u)<l,(Po) + u4,(Pt) 3the subdivision algorithm is also referred to as de Casteljau algorithm 16 = B , K [.«/>„), <«/>,)]) (ii) Assume, that the claim is true for m, then («> [P o> Pi, •• • • > Pm+l])) = (Bm(u , [P0, P i, Pn])) + < j > [pu P2, Pm+1])) ^ (1 - u)B m(u , [< /> (P 0), 0 (P ),..., <£(Pm)]) + uB m(u, [^ (P ), (Pm+ 1)]) ^ P m+i(u, [<£(P0), ^ (P ), 0(Pm+l)]) Thus, it is also true for (m+1), q.e.d. If the designer wants to look at the object from another angle or if he wants to change its size and position, he just maps the control points first, and then computes the Bezier curve. This is equivalent to computing the curve first and then mapping it. Yet, the latter m ethod is more costly, because every single point of the curve needs to be mapped, whereas in the first m ethod, only the control points are mapped. Note, that Bezier curves are not invariant under projections (see [7]). 2.2.8 Bm has the variation diminishing property In [4], Bezier defines that “a curve is variation diminishing if any plane intersects the curve in no more points than it intersects its control polygon.” Bezier curves have this property. For proving this we need some results obtained from the derivation 17 of the subdivision algorithm presented in the next chapter. Therefore, the proof is postponed until C hapter 3. T he variation diminishing property guarantees th a t the curve “follows” the path given by the control points, i.e. it makes the shape of the curve predictable (see Figure 2.3 for an example). T hat means, if the control polygon is “increasing”, so is its resulting curve. If the control polygon is convex, so is the Bezier curve associated with those control points. T he properties in 2.2.1 to 2.2.8 show, why Bezier curves are so popular in Com puter Aided Geom etric Design. They enjoy (almost) all desirable properties needed in the design process and they have an underlying concept, which makes their han dling easy. 2.3 Parametric Representation of Bezier Curves So far we could only describe functions from 1R to JR. But for designing it is more im portant to describe curves in the two dimensional plane. Therefore we use the concept of trajectories, where u becomes ju st a param eter th at varies along the length of the curve between P0 and Pm (see [10]). 18 Thus, instead of describing the Bezier curve by (u, £ IR2, u 6 [0,1], we now describe the Bezier curve by (P ^ fn ), B ^ u ) ) £ jK2, where B 'm(u) = 52 h,m{u)P'k , u £ [0,1] , i = 1,2 fc = 0 (2.9) and Pk = (P/1,P2) £ #22. The P *. are now points in the plane. This allows more flexibility and it will serve as the basis for the formulation and use of Bezier curves and surfaces in the following chapters. 2.4 Matrix Representation of Bezier Curves Sometimes the m atrix representation of Bezier curves is preferred over the definition given in (2.1). We need the m atrix representation of Bezier curves for an algorithm described in the next chapter. The following description can be found in [7] and is equivalent to (2.1): / \ Po B m{u) - 526m,*:Pfc— ( 6m0 6m,i ... bm < r k-0 \ and (2.10) P \ / um um- 1 ... u° } A , (2.11) 19 A = (/U ,™ -,)*.. » = U - l T k(m ] ( f \ \ J / \Kj k= 0...in j= 0,n .t m (2 .12) We can verify (2.12) by induction. Thus, we have Bm(u) = uTA p , where u = i um um-l _ 1 P = Po Pi (2.13) (2.14) (2.15) (2.16) For example, for m =3 we get: — I u3 u2 u 1 \ ( \ - 1 3 - 3 1 Po 3 - 6 3 0 Pi -3 3 0 0 Pi o o o \ P 3 ) The m atrix formulation can also be used for a change in representation from Bezier curves to B-splines (see [2]). 20 2.5 Bezier Surfaces Bezier curves in 1R2 can be easily generalized to Bezier surfaces in 1R3 by taking the tensor product of two Bezier curves B m and B n. Then we get 771 Tt B m,n(u, v;P) = J 2 J 2 bmik(u)bnJ(v)Pkij , («, v) £ [0,1] x [0,1] , (2.17) k=0j =0 where Pkj is the control mesh of B m < n and bm > k, bnj are the Bernstein polynomials from (2.1). In addition to Figure 2.4, examples of Bezier surfaces and their underlying control mesh will be given in Chapter 4. Next we want to present the matrix representation of B mtn(u,v] P). We use the same concept and the same notation as in Section 2.4, but we denote A by A m, A m G JRm+l’ m+1, and A n, A n € JRn+1,n+1, respectively. Recall that the control mesh P is a matrix P £ lRm+l'n+l, Then we obtain m m Bm > n{u,v\P ) = Y ^ Y l bm,k{u )K j(v)pk,j k=Qj— 0 m n = '£ .b™Au ) '£ , p^ b”j ( v ) (2-is) fc = 0 j-0 = U T A mP{vT A n f = uTAmPA„v, (u,v) e [0,1] x [0,1] Finally, let us also derive a parameterized version of Bezier surfaces. Instead of taking the point (u,u, Bm i„(u,u)) € IB? on the surface, we consider u and v 21 1v 0.8 v . ; 0 .6 - •• 0.4- ■■■•" 0.2 -1.5 0.4 0.3 -2.5 0.2 0.1 3 0 Figure 2.4: Example of a Bezier surface with control mesh to be param eters. Changing u means moving on the surface in .^-direction (with respect to the Cartesian basis) and, similarly, altering v means moving on the sur face in ^-direction. The control points can be viewed as a control mesh of points Pfij G IR3, k = 0 ,1 ,..., m, j = 0 ,1 ,...,n . Thus, we describe the points on the surface by (B ^ n(u ,u ), B ^ n{u,v), B®t„(w,u)) G IR3, where rn n Blm,n(U^ P ) = ]C ^ Z bm,k(U)bn A V)H ,j , i = 1 ,2 ,3 , (2.19) k~oi= o and Pk,i = (Pkj, P£j, P 3 j) G IR3 are the points of the control mesh of B m,n. 22 Using the param eterized version, we can see that Bezier surfaces inherit (almost) all properties of Bezier curves. We notice two exceptions: the variation diminishing property does not hold (actually the variation diminishing property is not defined in IR3) and the Bezier surface interpolates only the four corner points of the control mesh, but not the whole boundary curve (see Figure 2.4). 23 C hapter 3 The Subdivision Algorithm and its Implementation The last chapter provided all necessary tools to formulate the subdivision algorithm for computing Bezier curves and surfaces. We will now discuss this algorithm and its implementation in the C programming language. In public sources like the subroutine library Netlib one can find various imple mentations of B-spline algorithms. However, algorithms for Bezier curves are not available. Farin includes in [7] an implementation of the Horner scheme to compute a specific point of £fm, similar to polynomial evaluation. This algorithm is not ef ficient for computing an entire curve or surface, because it needs a Horner scheme evaluation for each point of Bm. An efficient algorithm for the computation of Bezier curves and surfaces is based on Property (2.7) in Chapter 2. Equation (2.7) gave a simple formula to write the Bezier curve Bm as the sum of two Bezier curves of degree (m-1). Barry and Goldman proved in [1] that the algorithm based on (2.7) 24 is unique to Bezier curves. This is the reason why all formulations of algorithm s for Bezier curves are based on this property. The uniqueness also explains why many modifications of Bezier curves have hardly any practical im pact (for modifications see, for example, T.V.To and H.N.Phien in [21] or N.M.Ntoko in [16] and [17]): there exists no general efficient algorithm for those m ethods. A successfull generalization of Bezier curves are the rational Bezier curves (see, for example, Piegl in [18]). One can derive an algorithm for rational Bezier curves sim ilar to our subdivsion algo rithm . An efficient formulation of the algorithm for rendering Bezier curves and surfaces was described by Lane and Riesenfeld in [13]. The formulation of the subdivison theorem and the resulting algorithm are based on this paper. Due to the uniqueness of the subdivision algorithm one cannot expect developing an efficient algorithm which is completely new. One can only try to derive a different form ulation of the algorithm. Section 3.4 will describe such a different approach to the subdivison algorithm based on the m atrix representation of Bezier curves. 3.1 The Subdivision Theorem We want to subdivide the Bezier curve £?m based on the (m +1) control points IJi, P2, P m into two curves with control points Q0, Q± ,..., Qm and Rq, R \ , R m. This new set of 2 (m + l) control points shall describe the same curve as Po, P i, ..., Pm- 25 T h e o re m 1 (S u b d iv isio n T h eo rem ) Let B m be defined as in (2.1). Then Bm(u, Pi, P„]) = B m(v, V n s p . i ] B m(v, K , P Z ~ ' P “]) V u€[l,l] (3.1) « e [o, l], where Pi if-) +J3 *"1 Pi k~l, 2, m k=0 P ro o f: we will prove the statement for u € [0, |] by induction on m. The equation for u € [^, 1] follows exactly in the same way. (3.2) Note that u = u £ [0, |], v £ [0,1], (i) Let rn=l: By(u, [P0,Pi]) (1 - u)Po + uPl = (1 - 2u)Pq + u(P0 + Pi) = (1 - «)P„ + Lv(Po + Pi) 3 =2 (1 - v)Pg + vPl = p Q » * ( l - v)'-*P£ = Bi(v,lPg,Pi']) 26 (ii) Assume that the claim is true for m-1, then P m(u, [PQ ,P i,...,P m ]) = (1 - « )P m_i(u,[P 0, P i,..., Pm_i] 4- uBm-i(u, [Pi, P2,...,P m]) ' # (i - \v )B n . ,(v, [p°, p , \ p%:;]) + kvB m _ ,(«, [ / f , p i p ” - 1 ]) = a -v)B,„-,(v,{PS,r;...... p ”.-/]) + (B m(v, [P0 °, P ‘,P ™ .- ,1 ]) + [P°, P2 ‘ .......P ” ->])) 1 (1- »)£„-,(«, [P„°, P,1, P ”.-,']) + [ S + a , a ± a , +g - D i 2 (1 - «)& »_,(«, IPg, P i P ^ , 1 ]) + [P,\ P2 2, P ” ]) is7 £ » ,(» ,P ^ .P ,',...,P S ) This shows that the claim is true for m, q.e.d. If we apply this subdivision iteratively on P0, P i , P m, we shall see that the set of those new control points will converge to the Bezier curve. One can observe that the subdivision property is largely based on (2.7). Figure 3.1 illustrates the subdivision theorem. Note that P m( |) = Pm is computed exactly. Bartels et al. show in [2 ] that this midpoint subdivision can be generalized to compute any point of the curve exactly. Thus, their algorithm is an efficient method if one has to compute only a very limited number of points on a Bezier curve. However, we are interested in an algorithm that computes an approximation of the entire curve associated with 1 B m is linear in P 27 1 0.8 0.6 0.4 0.2 0 the control points P0, P \ , P m. Then, their algorithm is inefficient because it is com putationally too costly. Therefore, we are looking for an efficient algorithm which approximates the Bezier curve up to any desired accuracy. It should be noted, though, th at the com putation of B m(u) using the the definition is unsuitable, because it involves powers and m evaluations of a binomial coefficient. 2note that P J c j in Figure 3.1 corresponds to our notation P£ I I i '----1 I T P_20=P 2 P 21 P 10=P 1 P 32 P_33 P 22 P_43 f t P -3 1 P 44 P_42 P 41 P 40=P 4 P 30=P 3 P_00=P 0 .2 0 0.2 0.4 0.6 0.8 1 Figure 3.1: Subdivision for a quartic Bezier curve2 28 3.2 The Generation of a Bezier Curve B m 3.2.1 The subdivision algorithm Let us form ulate an algorithm ic description of the subdivision theorem as described in [13]): Let P = (Po, P i,...,P m) be the original set of (m +1) control points. Then, we want to com pute the new set of control points Q = (Q0, Q i , Q m) = (P£, Px l, P ™ ) and R = (Po, Ri > Rm) = (Pm i ■ > Pm)- Note, that, due to the subdivision of [0,1] into [0, |] U [ |, 1] and the sym m etry relation (2.4), we can set P£ = Pjf, i.e. R m~k = Qk, V/c = 0, Consequently, we get the following subdivison algo rithm : A lg o rith m 1 (cu rv esp lit(P , Q , R , m )) copy Q = P Pm = Pm for j — 0:m — 1 tm p2=Qj for i — j :m — 1 tm p l= tm p 2 tm p2=|(<3t + Qi+1) Q ,= tm p l en d C?m=tmp2 P m_ i„j= tm p 2 end It follows th a t the algorithm curvesplit needs m — 1 m — 1 771— 1 m Z Z 1 = E ( m “ i) = E i = m (m + l) = 1 ^ 2 + h=Q k=j j= 0 j = l (3.3) 29 multiplications and additions, i.e. the cost of the algorithm is of order 0 (m 2) flops. It is im portant to mention that the algorithm involves only multiplications and additions. Therefore it is numerically stable. The following theorem is the basis of the curve generation algorithm: T heorem 2 Let P- be defined as in (3.2) and let M = m axi-0ii m_ 1|/5 t+1 — P,|, then \ P i - P i i } \ < i M Vi = 0,l,...,m-1 (3.4) Proof: • \P i - e f t ? I = I P i - a = P ± L \ = m + P f » I . * = 1 • I Pi - PUi I = iF;z!-i 5 " - < $ { \ p t l - ^ ; i i + 1 p t \ - p t l i) using this relation on the two terms of the right hand side successively gives I/1 ; - | < §(M + M) = M , where M = m ax{= Q A m_ !|Pi+ 1 ~ />-| • finally, we get |P/ — P-+) | < ~M, q.e.d. In addition this also holds for P™, P™~1, P() because of symmetry. Theorem 2 shows that we divide the maximum distance between any two consecutive control points with each curve splitting by two. After the k-th iteration we have 2k(m + 1) control points and the distance between any two consecutive control points is then at most (^)kM . Consequently, we can conclude that the sequence of control points generated by iterating the subdivision algorithm converges to the Bezier curve. 30 We are now able to prove the variation diminishing property of Bezier curves. Recall th at B m is variation diminishing if any straight line intersects the Bezier curve in no more points than it does intersect its control polygon. Theorem s 1 and 2 proved th at we can approxim ate B m by “refining” the control polygon, i.e. we approxim ate Bm by piecewise linear functions. Therefore, it is clear th at any straight line cannot intersect the Bezier curve in more points than its control polygon. Thus, we showed th at B m is variation diminishing. 3.2.2 Generation of a curve Bm The observations of the previous section yield the following recursive algorithm for approxim ating a Bezier curve: A lg o rith m 2 (c u rv e g e n (P , m , s to p p in g -c rite rio n )) if( stopping-criterion = = yes ) output P else curvesplitfP 1, P 1, m) curvesplit(P2, Q2, R 2, m) curvegen(Q, m, stopping-criterion) curvegen(R, m, stopping.criterion) e n d The algorithm uses the param eterized version of the Bezier curve description as described by (2.9). Remember th at P = (P * ,P 2) € IRm'2. With each step of the algorithm the set of control points is refined. Algorithm 2 splits P into Q and R, and then proceeds to split Q and R in the same way. The refinement of each of the two parts is done until the desired level of accuracy is reached. This means th at the 31 stops in Algorithm 2 correspond to the traversal of a binary tree. In the literature (e.g. in [13] and [2]), we can find various stopping criteria: 1. global error We showed, th at our curve splitting algorithm converges and that |Ti+i — T{\ < (^)hM, Vi = 0,1,..., m — 1, for the control points T) of the fc-th iteration, where M = maxi= Q ti — P* I, |^ + i ~ maximum distance of any two consecutive points in the original set of control points. If we preset the global error on tol, then we can precalculate the num ber of iterations needed to obtain an approximation within this error bound. By solving the equation tol — {~)kM we get stop if k > log(M) — log(tol) log{ 2) (3.5) This criterion requires just one addition per iteration (increase/decrease the counter for depth of iteration). Despite the low cost, this criterion is more of theoretical than of practical importance because it is too inflexible. For more details see Chapter 4. 2. m esh refinem ent Here, we use a preset maximum distance tol between any two consecutive control points in the &-th iteration, i.e. stop if m o ii- 0,1 m-ilTf+i — Til < tol, (3.6) 32 where T{ are the control points in the k-th iteration. This stopping criterion uses up to m additions per iteration. It might be used in programs for plotting Bezier curves, where tol is given by the resolution of the screen or the maximum accuracy of the system. 3. lo cal flatness The subdivision algorithm computes basically a piecewise linear approximation to the Bezier curve. Local flatness means that the control points are within a preset tolerance tol of being linear: slop if distance(Pk, P0Pm) < tol V /c = 1,2, 1 (3-7) where distance(Pk, PaPm) gives the Euclidean distance of P fr o m the straight line PoPm■ Let y = m ix — P j) + Pn, where m = is the slope of FoPm, M r0 and let yt- — — P£) -b P£ be the normal line through P/.. 33 Then we obtain distance(Pk, PoPm) — distance(Pk,yk H y) = y /(.P l ~ “ )2 + ( P I - !>)* where “ = x (Pt + m(mPi + pk ~ p3» b = m(a - P i ) + P i T he im plem entation of distance requires seven m ultiplications and nine addi tions. (see Appendix A .l for the exact im plem entation). Consequently, the local flatness criterion uses up to 7m m ultiplications and 9m additions. In m ost approxim ation problems we are interested in a piecewise linear approx im ation. This criterion gives a linear approxim ation up to a given tolerance. 4. com b in ation o f th e m esh refin em en t and th e local fla tn ess criterion As we will see in C hapter 4, it is most efficient to combine the two latter criteria. In addition we employ a m axim um depth of recursion. If this m axim um depth of recursion is reached, the program stops the curve splitting. This is necessary to avoid an excessive runtim e of the program. 34 Our implementation of the combination criterion uses the same tolerance for mesh refinement and local flatness. The program can be easily adjusted so that two separate tolerances for mesh refinement and local flatness are allowed. In C hapter 4 we will discuss the practical implications of those stopping criteria. 3.2.3 Implementation in C The program for the generation of Bezier curves (and surfaces) is w ritten in the C language and compiled with the Sun C-compiler “cc” (for more details about the C language and the compiler “cc” see [11], [12] and [14]). Double precision is used throughout. The implementation of the subdivision algorithm for generating a Bezier curve can be found in Appendix A.I. The program accepts control points stored in an array of size 2m, where P£ is stored in the first m places and P 2 in the last m places. One can select one of the stopping criteria discussed in 3.2.2 (default is local flatness) and enter the preset tolerance level tol. The approximation algorithm is implemented as described in the previous sections with two m ajor exceptions: the control points are not saved in a m atrix P E lRm’ 2, but in an array P E //?2,n’1 as described above. The curvesplit algorithm needs only two vectors P and Q E IR2m to use less memory. In Algorithm 1 (curvesplit) P is needed only for copying it to Q. Therefore, P is m anipulated directly in the implemented algorithm curvesplit instead of the copied version Q. This means that the program runs with two arrays P and Q, where P is allocated at the program start and Q is allocated dynamically at each step of the recursion. The output is of the same form as the input. However, the function curvegen (Algorithm 2) does also accept any other output function. These are the main implementation characteristics of the program. For more details see Appendix A.I. 3.3 The Generation of a Bezier Surface B m ,n 3.3.1 The subdivision algorithm Now, we have B min(u,v;P) - Y%=aI?j=ot>m,k{u)bn,j(v)Pk,j as in (2.17) with control mesh Pkj (see Figure 3.2). We use the same concept for the approximation of sur faces as we did for curves. Especially, since Bm > n(u,v,P ) = £/T=o ^m,k{u) EjLo ^n,j{v)Bk,j can be written as the product of two “independent” functions, we can look at the coordinates u and v separately, i.e. first, we split u and then we split v. Then, by subdivid- ing it, we split Pj (Pqj, • * ■ , P-m,j) into Qj — (Qo,j, Qi,ji and Rj ~ (/?o,j, Ri,j, ■ ■ ■ ■ > Rm,j), V y = 0,1, Afterwards v is subdivided, and, thus, we split Qk = {Qk,o,Qk,u—,Qk,n) into Sk = (Skto,Sk,i,...,Sk,n) and Tk = (Tkt0, TklU..., 74,n), as well as R k = (Rk,ot Bk,l, B k> n ) into Uk = (t/jt,o, Uk,i , Uk,n) and 14 = (14,o, 1 4 , i , 14,n), V A : = 0,1, ...,m . For all subdivisions we can use exactly the same 36 -3 1 Figure 3.2: Control mesh of a Bezier surface3 scheme as for Bezier curves. As a result we get the following algorithm: A lg o rith m 3 (su rfsp lit(P , S, T , U , V , m , n)) copy Q = P \* split in -u-direction *\ fo r k=0:n Rm,k = Pm,k for j = 0 : m — 1 tmp2=<2ji f c for i = j : m — 1 tm pl= tm p2 tmp2=f(Qi,fc + Qi+i,k) Qijt= tm p l en d Qm,f c =trnp2 3note that P_kj in Figure 3.2 corresponds to our notation P k j of the control mesh 37 R m- i-j,fc=tmp2 end end \* split Q in u-direction * \ copy S = Q for k— 0:m Tk,n ~ Sk,n for j = 0 : n — 1 tmp2=Sk,j for i = j : n — 1 tm pl= tm p2 tmp2=J(5'feI i + Sk,i+i) S k ,i= tm pl end 5, j t,m=tm p2 rjt,n_i_j= tm p2 end end \* split R in u-direction *\ copy U = R for k=0:m Vk,n = Uk,n for j = 0 : n — 1 tm p2=Uk,j for i — j : n — 1 tm p l= tm p 2 tm p2=f(t7fc ii + Uk,i+1) C4,;=tmpl end Uk,m= ^mp2 V fc ,n-i-j= tm p 2 end end In the original formulation by Lane and Riesenfeld in [13] two lines of the algorithm were missing and several indices were incorrect. We corrected these errors in the algorithm described above. 38 The splitting requires the following number of flops: • split in ^-direction: E£=0 E ^ q 1 E S j 1 1 = (n + flops • split in u-direction: 2E£L0 Ej=o 1 = (m + l)n(n + 1) flops Consequently, Algorithm 3 uses X(n + l)m(m + 1) + (m + l)n(n + 1) = (n + l)(m + 1)(?2 -f \m ) « n 2m 4- |n m 2 (3-8) additions and multiplications, i.e. the algorithm is of order 0 (n 2m + m 2n). The subdivision is again numerically stable because it involves only multiplications and additions. The convergence of the algorithm comes from the fact that we use ex actly the same procedure for the splitting in u— and u-direction as we did for the subdivision of Bezier curves. 3.3.2 Generation of a surface B m )n We use the parameterized version (2.19) of J9m,„ to formulate the algorithm for surface generation. It is again a direct generalization of Algorithm 2 for Bezier curves. A lgorith m 4 (surfgen(P, m , n, stopping_criterion)) if( stopping_criterion ==■ yes ) output P else 39 surfsplit( P \ S l , T l,U l , V \ m, n) surfsplit(P2, S 2, T 2, U2, V 2, m, n) surfsplit(P3, S 3, T 3, U3, V 3, m , n) surfgen(S, m, n, stopping.criterion) surfgen(T, m, n, stopping_criterion) surfgen(U, m, n, stopping_criterion) surfgen(V, m, n, stopping_criterion) end The stops in Algorithm 4 correspond to the traversal of a tree with four leaves. W ith each iteration the control mesh is refined and it approximates the actual Bezier surface up to the given tolerance in the stopping criterion. We regard the same stopping criteria as for Bezier curves. The first two criteria global error and mesh refinement remain virtually the same. We just need to replace M = m«.'ct=0,i....m-i|P{+i - Pi\ by M = m a x k = * ,t,...,m -i { | / V h - Pkti I, lA+i,; - j= o r l 1 i.e. we take M to be the maximum distance of any two consecutive control points of the mesh in u and v direction. Therefore we get: 1. global error stop if k > log(M) — log(tol) (3.9) log{2) where M is the maximum distance of any two consecutive control points and tol is the preset global error. 2. m esh refinem ent stop if maxk=o,i m -i {|Tjfcli+1 - Tfcii|, |7&+i,f - Thii\} < to l, (3.10) j=0»l n — 1 40 where T is the computed mesh in the k-th iteration and tol is the preset tolerance for the output mesh. 3. local flatness This criterion is also a straightforward generalization of the local flatness cri terion for Bezier curves. Algorithm 4 computes a piecewise planar approxi mation. Local flatness means that the control mesh is locally planar within a preset tolerance tol. Let PoflPm< a P o ,n denote the plane given by the three points P0,o, Pm,o and o- Also let PkPj denote the straight line through Pk and Pj. Then we get stop i f distance(Pkj, PO, 0,P m,oPo,n) <tol and distance(Pk$, , PafiPm.a) < iol and distance(P0j, , P0,oP0 > n) < tol and ^ Z q’ distance(Pm rj , , Pm,nPm,o) < tol and distance(Pk,n, , Pm,nPo,n) < tol W ith this criterion we test whether all points are within the given tolerance of being planar, and then we check whether the boundary lines of our approxi mation are within the given tolerance of being linear. Easy calculations give the following formulas for the Euclidean distances: distance(P, Q R S ) Remember that P,Q , R, S € Hi3, then denote a = (tii 02 ct3)T = R — Q 41 and b = (bi b2 b3)T = S — Q. The normal vector n of a and b is then given by / \ a2b3 — a3b2 n = \ o ,3bi — a 163 a\b2 ~ The Euclidean distance of P from QRS can be computed using the Hesse- Normal-Form: distance(P, Q R S ) = ( t t — r,n{P — Q) (3.11) distance(P, QR) Denote d = R — Q and e = P — Q, where d, e £ JR3. The distance of a point from a plane in 1R3 can be computed by finding P \ the point of intersection between PQ and the plane through P with normal vector d: distance(P, QR) = distance(P, P') F e dr e, T = e e — cP e dTe F d The implementation of the stopping criterion local flatness requires up to 19777.+ 1677 + 23 multiplications and 20m + 6n + 37 additions (see Appendix A.2 for the actual implementation). 4. com bination of th e local flatness and m esh refinem ent criterion Combining the criteria means stopping if either one of the above criteria re quires stopping the iteration (see 3.2.2). Note that we use again the same tolerance for both criteria. In Chapter 4 we discuss the practical implications of those criteria. 3.3.3 Implementation in C The program code of the surface approximation can be found in Appendix A.2. The implementation is once more a straightforward generalization of the algorithm for generating Bezier curves. The program needs the control mesh to be stored in an array of size 3(m -f l)(?i + 1), where the columns of the control mesh are stored in the array successively for the x, y and z coordinate, i.e. p __ l p i p i p i p i p i p i p i p 2 p 2 p 3 p 3 \ 1 \-* 0,0) 0,1) ■■•) 1 0,m> 1,0) *•') 7 1—1 ,TTi) ■ r n,0i ’**> n,m) 0,0) ■••) 1 n ,mi 1 0,0) •••) 1 n,m)' In addition, the program needs the size of the control mesh, i.e. it needs the param eters m and n. We have the same options for the stopping criteria (default is local flatness) and the preset tolerance tol. The implementation of the surface generation algorithm uses the same storage for the newly generated control mesh as for the original mesh P. The program does not make use of all the variables P, Q, R, S, T, U and V (as described in Algorithm 3), but it uses only two variables for storage and manipulation. Therefore, we do not 43 need to copy any data. Q, R and S are allocated dynamically in one variable for each iteration. The program accepts any output function, where the default output function is the successive printing of the computed planar approximation. This means that the output can be used for the pointwise representation of a surface. The program includes only rudimentary controls. For example, it checks for the control mesh only whether it could read (m + 1) x (n + 1) control points altogether but not whether each it could read (m + 1) points for each column of the control mesh. Further controls can be built in easily. The program is written in modular form so that it is flexible and can easily be adjusted to specific needs. This means that flexibility was more important than efficiency. These are the main implementation characteristics of the program. The perfor mance evaluation is discussed in the following chapter. 3.4 The Subdivision Algorithm in Matrix Representation In Section 2.4 we introduced the matrix representation of Bezier curves. We will now discuss the matrix formulation of the subdivision algorithm for Bezier curves 44 and surfaces. The description is a generalization'1 of an article by V.J. Harrand, J.P. Ziebarth and A. Choudry in [9]. 3.4.1 Formulation for Bezier curves We use the notation of (2.12), i.e. B m(u) = w TAp, but we simplify n(«, [Po, Pi j f ’ m]) to B m(ir, P ). We want to find a new set of control points P ', such that P ) = B m{v\ P') Vu € [0, |] ,u € [0,1]. Note th at the existence of Bm(v, P') is clear from the subdivision theorem. We know th at u = so we get Vu € [0, |] « = ( u™ it™-1 ... u° ) = ( ( |u ) m ( J « r - 1 ... (|u )° ) = S v S = cliag( ... ( i )0 ) and (3.13) v — ( vm v m-l ) We obtain uTAp = v TAp' or v TS A p = vTAp' from the requirement B(m\ P ) = B m(v\ P'). A is invertible, because A = (A *..ifn_j-)fc=o,i m = A J, where A = (Ak,j)k=o.i » > and J = (Ik,n-j) k=o,i m with * j s O , L , . . . , r r t j = 0 , l r n j s s O pl m J 2 = /, where / denotes the identity matrix; A is an upper triangular m atrix with only non-zero elements on the diagonal.5 Hence, AJ is upper triangular and A can be w ritten as the product of two invertible matrices, i.e. A is invertible. Consequently 4 [9] describes the surface generation only for the case m =n=3; our description allows a control mesh of arbitrary size m and n 5remember that A k ,m - j = (— (j.) and (£) = 0 V A * > j 45 we get p' = J lSAp. Let us define A' = A ~lS A (3.14) then we obtain B m(u;P) = Bm(v ,P ,) = v TA 1SA p — vTA'p \fu 6 [0,1], v e [0, 1] (3.15) Thus, we expressed B m{u; P ) , u e M l , in terms of B m(v; P'), v € [0,1]. Now let us consider B m(u ;P ), u € [f,l]. Then we want to find a new set of control points P" such that Bm(u \P ) — B m(v,P ")y u S [5,1], v € [0,1]. For computing P" we apply the symmetry argument (2.4) similar to the subdivision algorithm. Thus, we reverse the order of our control points which means we reverse the order of the columns in A. Let us define p' = A'p (3.16) A" = A 'J, (3.17) then we obtain £?m(u; P) = Bm(v, P") = vT A 'J p = vTA"p, V« € [J, 1], u G [0, 1] (3.18) p" = A”p (3.19) (3.19) 46 where J is defined as above. Since A! and A " do not depend on p, the two matrices can be precalculated, i.e. computed once and for all. The numerical stability of the matrices A! and A " needs to be investigated further, because matrix inversion is involved. If we assume that they can be computed in a stable way, we will obtain a new set of control points pf and p" by two m atrix/vector multiplications. Hence, the algorithm requires6 m +1 2 ^ 2 k = (m + 1 )(m + 2) = m 2 -f 3m + 2 (3.20) k=o multiplications and additions, i.e. the cost of the algorithm is of order 0 (m2). Remember that the subdivision algorithm (Algorithm 1) was also of order 0(?t?2), but it required only \(m 2 + m) additions and multiplications. We conclude that the m atrix formulation needs the precalculation of A ! and A " (including the inversion of a matrix) and then each subdivision requires more flops than Algorithm 1. 3.4.2 Formulation for Bezier surfaces For Bezier curves we used a vector with control points. For Bezier surfaces we have a control mesh, i.e. a matrix P = (Pk,j) fc= o,t m, of control points. The following .7= 0,1 description is based on the matrix description of Bezier surfaces as given in (2.19). Splitting in ^-direction corresponds to splitting each column in the matrix P. Using (3.15) and (3.18) we get: 6note that A J is upper triangular; hence, A is the matrix A turned by 90 degrees 47 • Vu £ [0, i] : B m,n{u,v\P) = Bm > n(s,V] P') = sTA'mP A ^v, s £ [0, 1], V £ [0, 1] where P ‘ = A ^ S A mP = A'mP • V u £ [ | , 1] : B m < n (u ,v;P ) = B m < n(s,v,P ") = sTA , J nP A lv , 5 6 [0,1], v £ [0, 1] where P" = A'mJ P = d" P In this way we split the control mesh P into P' and P". Splitting in u-direction means splitting each row of the matrices P' and P ". Using again (3.15) and (3.18) we get: • Vu € [0,i], v £ [0, J] : Mn,n(«, v; P) = Bm,n(s, i- P W ) = STA'mPA'Jv, 3 £ [0, 1], t £ [0, l] where p(m) = A'mPA'^ • Vu € [0, ±], v £ [£, 1] : B min(u,v; P) = Bmtn(s, t; P (n,)) = sTA'mP A^Tv, s € [0,1], t £ [0,1] where ph") = A'mPA"T • Vu € [J,l], v £ [0,^] : Bm,n(u, v; P) = B m > n(s, t ; P M) = sr A , J nP A , J v , s £ [0,1], t 6 [0,1] where P < * > = A'^PA'? 48 • Vu E [§,1], v € [§,1] : B m,n(ut v]P ) = B m ^ s ^ p W ) = sTA ^ P A ^ v y s E [0,1], t E [0,1] where P^O = A" PA "r Using the above results we derived the following algorithm: A lgorith m 5 (su rfsp lit(P , Q, R , S, T , m , n)) T m p l = A'mP Tm p2 — A '^P Q = T m p l A'T R = T m pl A"r S = Trnp2 A f T = Tm p2 A"7' Algorithm 5 subdivides the control mesh P into the four meshes Q , R, S , and T. Note that we have a m atrix Am E ^ m+1 - rn+1 for the split in w-direction and a m atrix A n E iPn+1,n+1 for the split in u-direction. This means that we have to precalculate the four matrices A ' m, A(^, A^, A", where inverting the m atrices A m and An is included. Then, the new control mesh is given by p (tu\ P^tv\ P^v\ P^vtK This formulation of the subdivision requires7 • for the split in u-direction (computing A ’ mP and A" P) 2(« + 1) E/TJi1 k = (n + 1 )(m + 1 )(m + 2) flops • for the split in ^-direction (computing T m p \A 'n and T m p 2 A ") 2(m 1) X)3 t=i & — (m + 1)(™ + l)(n + 2) flops. 7recall that A'm , A '^ , A ', and A " are upper triangular matrices turned by 90 degrees 49 Altogether, Algorithm 5 uses (n + l)(m -f l)(m + 2) + (m + l)(n + 1)(« + 2) fa n m 2 -f m n 2 (3.21) additions and multiplications. Thus, the algorithm is also of order 0{m.n2 + n m 2) but this means that the m atrix formulation needs more flops than the previous formulation in 3.4.2.8 For the actual implementation of this surface splitting we use the param eterized version of £?m,n, i.e. we use Algorithm 5 for the three coordinates x, y and z of the control mesh independently. In their article (see [9]) Harrand et al. use the m ethod described above for a control mesh of size m = n — 3. They showed th at their m atrix formulation of the subdivision algorithm is significantly faster if the m atrix formulated algorithm is implemented on a supercomputer (like CRAY or IRIS) using a hardware driven m atrix multiplier for matrices of size 4x4. Yet, this implementation loses its advantage if one has a control mesh of arbitrary size in and n. Having a control mesh of arbitrary size means th at the built in m atrix m ultiplier cannot be used directly any more, i.e. one has to go back to a software implemented m atrix m ultiplication.9 8Recall, that the previous formulation needed (n + l)(m + l)(n + ~m) flops 9another possibility would be to partition the problem by “breaking” the control mesh into smaller meshes of maximum size 4 x 4 50 In this section, we showed that a software based m atrix m ultiplication needs more flops than our implementation of the subdivision algorithm. In addition, we have to take into consideration the costs for computing the matrices A 'n, A A ' n and A", where two m atrix inversions are involved. Although the formulation of the algorithm in m atrix form is more elegant than the previous formulation of the subdivision algorithm, the m atrix formulation loses its advantages when arbitrary control meshes are needed. 51 C h ap ter 4 Test Results on Impeller Bladings from the Space Shuttle In theoretical analysis Bezier curves are used more often than Bezier surfaces. This stems from the fact that working with Bezier curves is easier than working with Bezier surfaces and that (almost) all properties of Bezier curves can be generalized to surfaces (see Chapter 2). In (real life) applications surfaces play a dominant role. Therefore, we will solely test our subdivision algorithm using control meshes to generate surfaces. The data used for the tests comes courtesy from the Rocketdyne Division of Rockwell International Cooperation. The data consists of control meshes to produce bladings for the first stage High Pressure Fuel Turbo Pump impeller in the Space Shuttle Main Engine. We have three bladings, each with two control meshes, one for the front and one for the backside of the blade. Figures 4.1 - 4.3 show those bladings. Appendix B .l depicts figures of the respective control meshes. Blade 1 in Figure 4.1 consists of the control 52 -3 1 Figure 4.1: Blade created by control meshes 1 and 4 meshes 1 and 4. Both control meshes have size 50 x 11, i.e. they have 50 x 11 = 550 control points each. Blade 2 in Figure 4.2 is formed by the control meshes 2 and 5, both of size 75 x 11, i.e. each control mesh has 825 control points. Finally, Blade 3 in Figure 4.3 is created by the control meshes 3 and 6. Those two control meshes have size 100 x 11, i.e. they consist of 1,100 control points each. All three figures were generated using the implemented algorithm. We choose a constant recursion of depth three to obtain a uniform distribution of the output points. Figure 4.1 consists of 70,400 output points generated in 62 sec runtime, 53 _4 u.o Figure 4.2: Blade created by control meshes 2 and 5 whereas Figure 4.2 is generated by 105,600 output points that were computed in 98 sec runtime. Figure 4.3 has the largest underlying control mesh. It consists of 140,800 points computed in 132 sec. All figures were plotted in about 15 sec using MATLAB (see [15]). The tests were taken using the program given in Appendix A.2. The program run on a Sun SparcStation 10 workstation with 128MB RAM and machine preci sion 2.22 ■ 10-16 (see Appendix A.3). A sample of the test data that was produced by the algorithm can be found in Appendix B.3. 54 -3 0 Figure 4.3: Blade created by control meshes 3 and 6 First we will discuss the performance of the subdivision algorithm and then we will analyze in detail the results obtained when the four different stopping criteria described in (3.3.2) are employed. Note that the runtime values were obtained using the UNIX tool “prof” . They vary slightly from one program execution to another. 55 4.1 General Performance Evaluation of the Subdivision Algorithm One can observe th at the program needs approximately 0.9 msec per output point independently of the control mesh or the stopping criterion (see Appendix B.2). Consequently we can obtain an a priori estim ate for the runtim e of the program that uses a stopping criterion with a known and fairly constant depth of recursion. Then, we know that there are 4fc (m + l)(n + 1) control points in the k-th iteration. Thus, the program needs approximately 0.Qmn4k msec for a recursion of depth k. We can conclude that the number of output points grows exponentially with a growing depth of recursion and so does the runtim e. This means th at subdivision for surfaces with large control meshes and a high required depth of recursion is com putationally very expensive. To increase efficiency, one should, if possible, subdivide the control mesh into smaller parts or use smaller control meshes in the design process, and then combining the different parts. Next we want to analyze the performance of our algorithm on an example. We choose control mesh 3 and stopping criterion local flatness. The data for the following discussion and results of numerical experiments are combined in Table 4.1. Similar tables for the other control meshes can be found in Appendix B.2. We observe that the runtim e for one call of surfsplit remains constant with changing tolerance. This is obvious because the algorithm surfsplit depends solely on the size of the control 56 control mesh: mesh 3 stopping criterion: local flatness # of control points: 100 x 11 = 1,100 criterion tolerance 0.1 0.05 0.01 0.005 0.001 0.0005 0.0001 # of out put points 30,800 60,500 436,700 941,600 3,782,900 9,171,800 45,844,700 runtime (in sec.) 36 55 396 848 3,374 8,234 40,812 max. depth of recursion 3 3 6 7 8 9 10 min. depth of recursion 2 2 3 4 5 6 7 # of calls surfsplit 27 54 396 855 3,438 8,337 41,676 msec/call surfsplit 78.6 81.0 78.3 78.8 79.3 78.8 78.3 # of calls surfgen 37 73 529 1,141 4,585 11,117 55,569 msec/call surfgen 1.1 1.0 1.4 1.6 1.3 1.5 1.3 Table 4.1: Performance analysis for mesh 3 mesh, but not on the level of accuracy. The runtime of surfsplit, i.e. the CPU time used by the algorithm, accounts for approximately 8% of the total runtime of the program, whereas the algorithm surfgen accounts for less than 1% of the total time. Since surfsplit is the function with the largest proportion of CPU time (see Appendix B.3), it follows that algorithm surfsplit dominates the runtimeof the program. Figure 4.4 comprises a graphical representation of the runtime and the number of output points versus different tolerances as given in Table 4.1. Note that the graph is in a double logarithmic scale. Figure 4.4 shows that the growth of the number of output points and the growth in the runtime is linear with slope approximately 57 # o f output points runtime log(tolerance) Figure 4.4: Comparison of runtime versus output for control mesh 3 — 1 in a double logarithmic scale. Consequently, tolerance and runtime as well as output are inversely proportional. This means, for example, that, whenever the tolerance changes by a factor 10“*, then output and runtime change by a factor of 10fc . Furthermore, since the graphs for runtime and output are parallel, Figure 4.4 shows that the runtime and the output grow at at the same rate, i.e. runtime and output are proportional. We obtain similar results if we plot the number of calls of surfsplit and the number of surfgen versus tolerance. These two factors are also inversely proportional to the tolerance. Each iteration, i.e. each call of surfgen, results in three calls of surfsplit (one call for each coordinate direction) and four additional calls of surfgen with the respective refined meshes. Consequently, the 58 criterion control mesh mesh 1 mesh 2 mesh 3 (m + 1) x (n + 1) 50 X 11 75 x 11 100 x 11 time (in msec) 24.0 47.2 77.7 k flops 16.9 34.8 58.9 time k flops 1.42 1.36 1.32 Table 4.2: Computational cost of the subdivision algorithm for surfaces number of calls of surfgen, surfsplit and the number of output points are directly related. These results show that it is sufficient to concentrate on the change of the number of output points in the discussion of Section 4.2. The runtime, number of calls of surfsplit and the number of calls of surfgen change similarly. Now let us consider the depth of recursion that the algorithm attained for the stopping criteria local flatness. The maximum (minimum) depth of recursion is the maximum (minimum) recursion needed by the algorithm to stop the iteration, i.e. to print a part of the surface. We observe in Table 4.1 that both, the minimum and maximum depth of recursion, grow linearly when the tolerance changes expo nentially from 10-1 to 10- '1 . Hence, we can conclude that the depth of recursion grows logarithmically for the tolerances given in Table 4.1. Comparable results can be obtained using other control meshes or other stopping criteria. Note that one output point consists of three real numbers. If we assume that each real number needs ten bytes for storage, then, for example, we need approximately 13 MB to store 3,782,900 output points for a tolerance of 0.001 in ASCII format (as given in Table 4.1). 59 Now let us verify the computational cost of the subdivision algorithm. For that purpose, we consider the ratio of the average CPU time needed for one call of the surfsplit algorithm and the number of flops needed in one call of surfsplit The number of k flops was computed using Formula (3.8). We can observe that the ratio remains constant within 5% (see Table 4.2). Hence, the runtim e of the surfsplit algorithm depends largely on the computational cost of the algorithm (see Formula (3.8)). Finally we want to mention that the program employed the routine “meshmax” to generate the test data. This routine accounted for approximately 3% of the total runtime of the program (see Appendix B). Hence, the total runtime can be reduced by 3% for the actual use of the algorithm. 4.2 Testing the Stopping Criteria In this section we will compare the efficiency versus the accuracy of the four different stopping criteria described in Section 3.3.2. We use the three control meshes mesh 1 with 550 control points, mesh 2 with 825 control points and mesh 3 with 1,100 control points, i.e. the sizes of the control meshes have the ratio 1.0 : 1.5 : 2.0. First we will discuss the results for a tolerance of 0.01. Tables 4.3 - 4.5 contain the necessary data. The attained error is the maximum distance between any two consecutive points of the approximated surface in the x, y and .s coordinate. In the global error criterion the depth of recursion is determined by the maximum 60 criterion local flatness mesh refinement combination # of out put points 31,900 36,850 26,950 max. depth of recursion 4 4 4 attained error 0.0174 0.0099 0.0174 runtime (in sec) 28 33 25 msec/call stop ping criterion 1.43 4.16 3.27 # of calls stop ping criterion 77 89 65 Table 4.3: Comparison of stopping criteria for mesh 1 and tolerance 0.01 criterion local flatness mesh refinement combination # of out put points 124,575 109,725 84,975 max. depth of recursion 5 4 4 attained error 0.0119 0.0099 0.0119 runtime (in sec) 112 106 79 msec/call stop ping criterion 2.24 6.38 3.68 # of calls stop ping criterion 201 177 137 Table 4.4: Comparison of stopping criteria for mesh 2 and tolerance 0.01 61 criterion local flatness mesh refinement combination # of out put points 436,700 192,500 192,500 max. depth of recursion 6 4 4 attained error 0.00989 0.00993 0.00993 runtim e (in sec) 396 178 177 m sec/call stop ping criterion 2.8 9.0 7.2 # of calls stop ping criterion 529 233 233 Table 4.5: Comparison of stopping criteria for mesh 3 and tolerance 0.01 distance between any two consecutive control points in the original control mesh. If the control points are not equidistant, this criterion is very inefficient, because parts of the control mesh have a smaller error than actually needed. Yet, it is the least expensive criterion, because it involves only two one bit shifts per recursion (increm ent/decrem ent of the counter variable for depth of recursion). Since the global error criterion is so inefficient, we did not include it in Tables 4.3 - 4.5 and we will om it it in our further discussion. Criteria like the mesh refinement criterion are in general more efficient overall despite their higher costs per function call since they result in a nearly uniform error. For a tolerance of 0.01 and a smaller control mesh like mesh 1 we observe that the local flatness criterion achieves an error that is considerably higher than the error attained by the mesh refinement criterion. Yet, the latter criterion yields a slightly slower runtim e and results in more output points. The relationship between 62 these criteria changes if we examine larger control meshes. For both larger meshes, mesh 2 and mesh 3, we obtain better results with the mesh refinement criterion: the mesh refinement criterion results in an up to twice as fast runtime (for mesh 3) compared to the local flatness criterion and it achieves the same error although the mesh refinement criterion needs approximately three times the runtime per function call the local flatness criterion needs. Hence, we can conclude that the cost of the stopping criterion can be neglected compared to the overall cost of computing the approximating surface. The combination of mesh refinement and local flatness achieves the best results for the smaller meshes mesh 1 and mesh 2. This combination criterion results in the same error as the local flatness criterion for both cases but the runtim e is shorter as well as the number of output points is smaller, i.e. the combination criterion achieves the same accuracy as the local flatness criterion but it is more efficient. For the larger control mesh 3 we observe that the combination criterion yields almost the same results as the mesh refinement criterion. It means that, for a tolerance of less than 0.01 between any two consecutive control points, the output mesh is never locally flat within a tolerance of 0.01. This explains why the local flatness criterion results in approximately the same attained error as the mesh refinement criterion but needs twice the output and runtime. Let us consider the stopping criteria for a tolerance of 0.001. Tables 4.6 - 4.8 contain the relevant data for the subsequent discussion. The conclusions of our 63 criterion local flatness mesh refinement combination of out put points 299,200 4,013,350 299,200 max. depth of recursion 6 7 6 attained error 43.5 • 10" 1 i — * © o o 1 43.5 • 10"' runtime (in sec) 262 3817 269 msec/call stop ping criterion 1.43 4.17 1.27 # of calls stop ping criterion 725 9,729 725 Table 4.6: Comparison of stopping criteria for mesh 1 and tolerance 0.001 numerical experiments change if we choose a smaller tolerance of 0.001. Again the global error criterion is the least efficient criterion. However, the local flatness criterion is now advantageous for all three control meshes, with regard to efficiency and compared to the mesh refinement criterion. The attained error of the local flatness criterion for the the control meshes 1 and 2 is approximately five times larger than the error of the mesh refinement criterion. The runtime of the program for the control mesh 1 using the mesh refinement criterion and, hence, the number of output points, is about 15 times larger than the runtime of the program with stopping criterion local flatness. It is still almost four times larger than the runtime of the local flatness criterion if we use control mesh 3. The runtime for one call of the stopping criterion is again approximately three times larger for the mesh refinement criterion than it is for the local flatness criterion. 64 criterion local flatness mesh refinement combination # of out put points 1,243,280 9,274,650 1,233,380 max. depth of recursion 8 7 7 attained error 54.4 • 10"“ 10.0 • 10- “ 54.4 • 10““ runtim e (in sec) 1,108 8,504 1,112 msec/call stop ping criterion 1.99 6.32 2.06 # of calls stop ping criterion 2,009 14,989 1,993 Table 4.7: Comparison of stopping criteria for mesh 2 and tolerance 0.001 criterion local flatness mesh refinement combination # of out put points 3,782,900 14,445,200 3,360,500 max. depth of recursion 8 7 6 attained error 28.0 • 10"“ 7.0 • 10"“ 28.0 • 10" “ runtim e (in sec) 3,374 13,262 3,068 m sec/call stop ping criterion 2.7 8.8 3.9 of calls stop ping criterion 4,585 17,509 4,073 Table 4.8: Comparison of stopping criteria for mesh 3 and tolerance 0.001 65 The combination criterion is again the most efficient of the four criteria but here it is of the same order as the local flatness criterion. This means that, if the ap proximated surface is now planar within a tolerance of 0.001, the maxim um distance between any two consecutive control points is larger than 0.001. Consequently, we get results th at are contrary to the tests with tolerance 0.01. Note that the runtim e of the combination criterion depends on which of the two criteria local flatness or mesh refinement is executed first. Our implementation uses, somewhat arbitrarily, the local flatness criterion first. Finally a word of caution. We compared the stopping criteria with regard to accuracy and efficiency. The mesh refinement and global error criteria are based on a given accuracy that needs to be attained by the approximation. Yet, the local flatness criterion and, similarly, the combination criterion do not take into account the maxim um error. They stop if the local approximation of the surface is fiat. This means that they are influenced by the complexity of the geometry in the approxima tion, whereas, if we use the mesh refinement criterion, it is not. Consequently, we have to be careful with interpreting the results and -besides accuracy and efficiency aspects- we also have to choose the stopping criterion which fits best the the needs of a specific geometry. 66 C h ap ter 5 Conclusions Nowadays CAD/CAM methods are standard in every manufacturing process. Al though Bezier curves are historically the first method in Computer Aided Geometric Design, an algorithm for generating Bezier curves and surfaces is not available in the public domain. Therefore, we decided to implement such an algorithm as a so phisticated exercise with great practical implications to gain experience in modern CAD methods. First we gave an introduction to Bezier curves. We focused on the properties of Bezier curves (and surfaces) which show that Bezier curves are highly practicable and easy understandable. Then, we discussed the subdivision algorithm for generating Bezier curves and surfaces. We saw that the subdivision algorithm is an efficient and numerically stable method for approximating Bezier curves up to any desired accuracy. Thereafter, we generalized the algorithm for generating Bezier curves in a straightforward manner to derive an algorithm for approximating Bezier surfaces. 67 The subdivision algorithm is based on property (2.7), i.e. it is based on the fact that a Bezier curve of degree m can be blended by by two Bezier curves of de gree m — 1. We mentioned that this property is specific to Bezier curves. Therefore, we could not expect to derive a new and efficient algorithm for approxim ating Bezier curves based on property (2.7). However, we introduced an element of novelty in generalizing a different formulation of the subdivision algorithm (see [9]). This new m atrix formulation gave an elegant and concise expression of the subdivision algo rithm . Unfortunately, the new formulation was computationally not advantageous in the current com puter architecture; it should become competitive when the algorithm is implemented on supercomputers with hardware driven m atrix multipliers. Finally, we tested our algorithm on control meshes of the bladings from the first stage High Pressure Turbo Pum p of the Space Shuttle Main Engine. We saw that the num ber of output points for the approximation and the runtim e grew inversely proportional to the required accuracy. Afterwards we concentrated on the efficiency and accuracy of different stopping criteria for the surface approximation. The algo rithm was flexible enough to allow the use of various stopping criteria for the curve or surface approximation. This implied that one could employ a stopping criterion which suited the desired application best. We tested various stopping criteria: we showed that the global error stopping criterion was inefficient for surface approx imation. A more efficient criterion than the global error criterion using the same underlying concept was the mesh refinement criterion. It term inated the recursion if 68 any two consecutive control points are within a given tolerance. We also tested the local flatness criterion which term inates the recursion if the approxim ation is locally flat. The test results showed that comparing the stopping criteria with respect to efficiency depends on the desired accuracy and the size of the control mesh. Yet, the deployment of one of the stopping criteria should largely depend on the definition of tolerance needed for a specific application. To increase the efficiency one can use a combination of various criteria. The use of our implemented algorithm is not necessarily restricted to curve and surface approximation. For example, due to the convex hull property of Bezier curves it can also be used to determine whether two curves or surfaces intersect (see [13]). Furthermore, the subdivision algorithm was first used in the automobile industry where most surfaces could be modeled using rectangular control meshes. To gain flexibility in the objects to be modeled, one could also implement a subdivision algorithm using triangular control meshes. And to broaden the concept even further, we could implement a subdivision algorithm for rational Bezier curves because there is practical evidence that rational Bezier curves are more flexible than traditional polynomial Bezier curves. In addition, rational Bezier curves and surfaces follow their control mesh more closely. Besides Bezier curves, there are other, similar, methods used in Com puter Aided Geometric Design. These methods consist in particular of B-splines and rational 69 splines.1 It should be mentioned that B-splines are the basic method in polynomial approximation and that one can show that Bezier curves are just a special case of B-splines. However, the concept of Bezier curves is easier to understand in practice and flexible enough for most real life applications. Although the Bezier curve method is the oldest in CAGD, Bezier curves are still heavily in use and research in this field is very active. 'for an introduction to splines in Computer Aided Geometric Design see for example [2] 70 Reference List [1] P.J. Barry and R.N. Goldman, “De Casteljau-type Subdivision is Peculiar to Bezier Curves,” in Computer Aided Design, Vol.20, 1988, 114-116 [2] R.H. Bartels, J.C. Beatty, B,A. Barsky, An Introduction to Splines for Use in Computer Graphics and Geometric Modeling, Morgan Kaufmann Publishers Inc., Los Altos, California, 1987 [3] P.E. Bezier , “UNISURF, from Styling to Tool-Shop,” in Computers in Industry, 4, 1983, 115-126 [4] P.E- Bezier , “The First Years of CAD/CAM and the UNISURF CAD System,” in Fundamental Developments of Computer Aided Geometric Modeling, L.Piegl, ed., Academic Press Inc., San Diego, 1993, 13-26 [5] P. de Casteljau, “Polar Forms and Surface Modeling at Citroen,” in Fundamen tal Developments of Computer Aided Geometric Modeling, L.Piegl, ed., Aca demic Press Inc., San Diego, 1993, 1-12 [6] S.V. Chasen, Geometric Principles and Procedures for Computer Graphic Ap plications, Prentice Hall Inc., Englewood Cliffs, N.J., 1978 [7] G. Farin, Curves and Surfaces for Computer Aided Geometric Design, Second Edition, Academic Press Inc., San Diego, 1990 [8] A.R. Forrest, “Interactive Interpolation and Approximation by Bezier Polyno mials,” Computer Journal, 15, 1972, 71-79 [9] V.J. Harrand, J.P. Ziebarth and A. Choudry, “Fast Subdivision of Bezier Sur faces,” in Simulation, May 1990, 241-244 [101 D. Kahaner, C. Moler, S. Nash, Numerical Methods and Software, Prentice Hall Inc., Englewood Cliffs, N.J., 1988 [11] B.W. Kernighan, D.M. Ritchie, The C Programming Language, Second Edition, Prentice Hall Inc., Englewood Cliffs, N.J., 1988 [12] B.W. Kernighan, R. Pike, The UNIX Programming Environment, Prentice Hall Inc., Englewood Cliffs, N.J., 1984 71 [13] J.M . Lane and R.F. Riesenfeld, “A Theoretical Development for the Com puter Generation and Display of Piecewise Polynomial Surfaces,” in IE E E Transac tions on Pattern Analysis and Machine Intelligence, Vol. PAMI-2, N o.l, Jan uary 1980, 35-46 [14] P. van der Linden, Expert C Programming, Prentice Hall Inc., Englewood ClifFs, N.J., 1994 [15] M A TLA B Reference Guide, The M ath Works, Inc., 1992 [16] N.M. Ntoko, “A Formulation for Bezier-Type Quintic Curves,” in Computers in Industry, 12, 1989, 67-72 [17] N.M. Ntoko, “A Formulation for Bezier-Type Curves,” in Computers in Indus try, 15, 1990, 363-368 [18] L. Piegl, “A Geometric Investigation of the Rational Bezier Scheme of Computer- Aided Design,” in Computers in Industry, 7, 1986, 401-410 [19] H.-P. Seidel, “A Geometric Approach to Bezier Curves,” in Computer Graphics in Mathematics, B. Falcidieno, I. Hermann, C. Pienovi, (eds.), Springer Verlag, Berlin, 1992, 31-43 [20] H.-P. Seidel, “An Introduction to Polar Forms,” in IE E E Computer Graphics and Applications, January 1993, 38-46 [21] T.V. To and H.N. Phiem , “Development of Bezier Based Curves,” in Computers in Industry, 20, 1992, 109-115 72 Appendix A Program Codes of the Implemented Algorithms A .l Program for the Generation of a Bezier Curve l * curve.h */ /* Header file with functions for curve generation and stopping criteria */ int n; / * int rec; / * double tol; / * double Min [2] ; /* double Max[2]; / * nate between any two consecutive control points */ variable for maximum distance in the x and y coord: nate between any two consecutive control points */ / * function to generate a Bezier curve as described in Section 3.2; ’P’ contains the control points for the curve generation; ’fp’ contains the file pointer for the output file; ’stoppingscriterion* is self explanatory * / / * USAGE void curvegenCFILE *fp, double *P, int (*func)Q stopping_criterion) */ void curvegenQ ; /* stopping criterion 5local flatness’ as described in Section 3.2; *P’ contains the control contains; function returns 1 if stopping criterion <= tol, 0 otherwise * / /* USAGE int difference.!(double *Px, *Py) * / 73 int difference_l(); /* stopping criterion 'mesh refinement’ as described in Section 3.2; ’P’ contains the control points; function returns 1 if stopping criterion <= tol, 0 otherwise * / / * USAGE int difference_d(double *Px, *Py) */ int difference_d(); /* stopping criterion that stops if depth of recursion is greater than given maximum depth of recursion (needed for the implementation of a fixed depth of recursion or the ’global error’ criterion; function returns 1 if stopping criterion <= tol, 0 otherwise */ / * USAGE int difference_f() */ int difference„f(); / * function that computes the maximum (minimum) distance in the coordinates x and y between any two consecutive control points for the given file with file pointer ’fp’; results are written in variables ’Mind’ and ’Max[] ’ * / /* USAGE void meshmax(FILE *fp) */ void mesh_max(); 74 /* curve.c * f # include 1 1 curve.h" # include "peripher.h" # include <stdlib.h> # define MAX(x,y) ((x) > (y) ? (x) : (y)) # define MIN(x,y) ((y) > 0 ? ((x) > (y) ? (y) : (x)) : (x)) # define ABS(x) C(x) > 0 ? Cx) : -(x)) / * self explanatory * / mesh_max(fp) FILE *fp; • c double xold, xnew, yold, ynew; Max[0]=0.0; Max[l]=0.0; Min[0]=10e8; Min[i]=10e8; rewind(fp); fscanfffp, "7,le '/,le\n", &xold, feyold); while (fscanf(fp, '"/,le 7,le\n", fcxnew, &ynew) > 0) Max[0]=MAX(Max[0], ABS(xnew-xold)); Max[1] =MAX(Max[l], ABS(ynew-yold)); Min[0]=MIN(Min TO], ABS(xnew-xold)); Min[1]=MIN(Min [1], ABS(ynew-yold)); xold=xnew; yold=ynew; > \* function that computes the distance of the point Q=(Qx, Qy) from the straight line given by y=m(x-Px)+Py and returns this value */ double distance(Px, Py, m, Qx, Qy) double Px, Py, m, Qx, Qy; { double tmpl, tmp2; tmp=(Qx+m*(m*Px+Qy-Py))/(m*m+i); Qx-=tmp; Qy-=m*(tmp-Px)+Py; 75 return Qx*Qx+Qy*Qy; > /* stopping criterion ’local flatness as described in Section 3.2 */ int difference_l(Px, Py) double *Px, *Py; ■ C int k=0; double x, y, ra; x=*Px; y=*Py; m=CPy[n-l] -y)/(Px[n-l]-x); while(++k < n-l) if (distance(x, y, m, *++Px, *++Py) >= tol) return 0; return 1; } /* stopping criterion ’mesh refinement’ as described in Section 3.2 */ int difference_d(Px, Py) double *Px, *Py; - C int k; for(k=0; k<n-l; k++) ■ C if ((ABS(*(Px+l) - *Px) >= tol) || (ABS(*(Py+l) - *Py) >= tol)) return 0; Px++; Py++; > return 1; } / * self explanatory */ int difference_f() if (rec < tol) return 0; else 76 > return 1; /* implementation of the ’curvesplit’ algorithm as described in Section 3.2 * / void curvesplitCQ, R) double *Q, *R; int i, k; double tmp1, tmp2; R+=n-l; *R=Q[n-l] ; for(k=l; k<n; k++) ■ c tmp2=*Q++; for (i=-l; i<n-k-l; i++) ■ C tmpl=tmp2; tmp2=Cq[i]+QCi+l] )/2.0; Q[i] =tmpl; y q[i] =tmp2; *— R=tmp2; > > /* ’curvegen’ algorithm as described in Section 3,2 * / void curvegen(fp, P, func) FILE *fp; double *P; int (*func)(); double *R; ++rec; i f CC*func)CP, P+n)) OUTPUTC P) else ■ C R= Cdouble *) callocC2*n, sizeof Cdouble)); curvesplitCP, R); curvesplitCP+n, R+n); curvegenCfp, P, func); curvegen(fp, R, func); free(R); /* peripher.h */ /* Header for input and output functions of control points (implementation dependent) */ / * Input written as sequence of control points in a file is accepted * / / * Output written into a file to meet MATLAB specifications to plot the computed curve */ / * Macro that prints the control points *P* into the file with filepointer ’fp’ */ # define OUTPUT(P) { int k; \ for(k=0; k<n-l; k++) \ fprintfCfp, '"/,le '/,le\nM, \ *(P+k), *(P+n+k)); / * function that computes the number of control points m in file ’infile’; function returns 0 if no error occured, 1 otherwise * / f * USAGE: int countline(char *infile) */ int countline(); /* function that reads the control points P=(Px, Py) from file with filepointer ’fp’; function returns number of control points read * / /* USAGE: int input(FILE *fp, double *Px, *Py) * / int inputO; 79 /* peripher.c *1 # include <stdio.h> # include <string.h> /* implementation as described in peripher.h */ int countline(infile) char *infile; • C FILE *fp; int count; char string[6]; char tmp [L_tmpnam]; / * function utilizes UNIX tool ’wc’ for counting input lines */ tmpnam(tmp); strcpyCstring, "wc -1 "); strcat(string, infile); streat(string, " > "); strcat(string, tmp); system(string); if ((fp=fopen(tmp, "r")) == NULL) return 0; else fscanf(fp, '"/d fecount) ; fclose(fp); remove(tmp); return count; int input(fp, Px, Py) FILE *fp; double *Px, *Py; • C int ct, count=0; while((ct=fscanf(fp, " '/ , l e '/,le\n", Px++, Py++)) > 0) count+=ct; return count; > 80 /* main.c * / /* Program for Bezier curve generation * / / * USAGE: main [ TOL | -d[TOL] | -f[TOL] | -m[TOL]] [infile outfile I in/outfile] */ / * the flags stand for the different stopping criteria: -d ’mesh refinement’ -m ’global error’ -f ’maximum depth of recursion’ default is ’local flatness’ in addition a specification of a tolerance for the stopping criterion is optional as well as of an input and output file * / # include <math.h> # include <stdio.h> # include <stdlib.h> it include <string.h> # defi: # defi: # defi: # defi: OUT_FILE "out.dat" /* INJFILE "in.dat" / * MAX_LE 24 /* T0L1 0.0001 / * TDL2 0.005 / * MAX„REC 5 / * # define MAX(x,y) (Cx) > (y) ? (x) # define ABS(x) (Cx) > 0 ? (x) : criterion ’local flatness’ * / default tolerance for stopping criterion ’mesh refinement’ */ default depth of recursion for stopping criteria ’global error’ and ’maximum depth of recursion’ */ : (y)) -Cx)) main(argc, argv) int argc; char **argv; /* file pointer for input and output file */ /* file name for input file */ /* file name for output file * / /* variable for control points * / int (*func)()=difference_l; / * variable for stopping criterion (default is ’local flatness’) * / f * variable for ’global error’ criterion */ / * flag needed for ’global error’ criterion */ / * variable for last control point (as last and first control point of two consecutive refinements are equal, only the first n-1 points are printed. FILE *fp; char infile[MAX_LE]; char outfile[MAX_LE] double *P; double error; char flag=’x’; double p[2] ; SI Therefore the last control point needs to stored for the final print) */ /* initialization of global variables */ tol=TOLl*TOLi; rec=-l; if(— argc > 0) if(**++argv == ’-’) switch(*++*argv) { case ’d’ : /* stopping criterion fimc=difference_d; ’mesh refinement' */ if ((tol=atof(++*argv++)) <= 0) tol=T0L2; break; case ’f’: / * stopping criterion func=difference_f; ’max. depth of recursion’ */ if ((tol=atof(++*argv++)) <= 0) tol=MAX_REC; break; case ’m’: /* stopping criterion func=difference_f; ’global error’ * / if (Cerror=atof(++*argv++)) <= 0) tol=MAX_REC; else flag=’m’; break; default: perrorC"USAGE: main [ TOL | -d[TOL] I -f[TOL] -m[TDL]] [infile outfile | in/outfile]"); exit(); break; argc— ; > / * read tolerance for default stopping criterion ’local flatness’ */ else if(isdigit(**argv)) • c tol=atof(*argv)*atof(*argv++); argc— ; > / * input and output file handling * / switch (argc) ■ { 82 case 2: strncpyCinfile, *argv++, MAX_LE); strncpyCoutfile, *argv, MAX_LE); break; case 1: strncpyCinfile, *argv, MAX_LE); strncpyCoutfile, *argv, MAX_LE); break; case 0: strncpyCinfile, IN_FILE, MAX_LE); strncpyCoutfile, 0UT_FILE, MAX_LE); break; default: perrorC''USAGE: curve [ TOL | -d[T0L] | -f[MAX_R£C] -mETOL]] Einfile outfile | in/outfile]"); exitC); > /* reading the control points into the variable JP’ * / ifCCfp=fopenCinfile, "r")) == NULL) ■ C perrorC'Inputfile could not be opened11); exitC); > n=countlineCinfile); P= Cdouble *) callocC2*n, sizeof Cdouble)); if CinputCfp, P, P+n) != 2*n) fcloseCfp); perrorC'Data of inputfile could not be read correctly"); exitC); > fcloseCfp); mesh_maxCfp); ifCCfp=fopenCoutfile, "w")) == NULL) perrorC'Outputfile could not be opened"); exitC); > / * computing the maximum depth of recursion for the ’global error’ criterion */ if Cflag -= ’o’) ifCCtol=CceilClogCMAXCMax[0], Max[l])) - logCerror))/logC0.5))) <= 0) ■ c perrorC'Given max. error is greater than max. error of control points"); fcloseCfp); exitO; > /* generation of a Be2ier curve */ p [0] =P [n-1] ; p [1]=P [2*n-l] ; curvegenCfp, P, func); fprintf (fp, " * / t le ’ /,le\n", p[0] , p[l]); / * clean up * / fclose(fp); free(P); /* makefile for main.c */ main: main.o curve.o peripher.o cc -p -dalign -02 -o main main.o curve.o peripher.o -lm peripher.o: peripher.h cc -p -c peripher.c curve.o: curve.h peripher.h cc -p -c curve.c main.o: curve.h peripher.h cc -p -c main.c clean: rm -f main 85 A.2 Program for the Generation of a Bezier Surface /* surf.h */ I * Header file with functions for surface generation and stopping criteria * / int M, ra, n; int rec; double tol; double Max[3] double Min[3] / * variables for mesh size, where M=m*n */ / * counter variable for depth of recursion */ / * variable for tolerance of stopping criterion */ / * variable for maximum distance in the x, y and z co ordinate between any two consecutive control points */ / * variable for minimum distance in the x, y and z co ordinate between any two consecutive control points */ /* function to generate a Bezier surface as described in Section 3.3; ’P’ contains the control mesh for the surface generation; ’fp’ contains the file pointer for the output file; ’stopping_criterion’ is self explanatory * / / * USAGE void surfgenCFILE *fp, double *P, int (*func)() stopping_criterion) * / void surfgen(); /* stopping criterion ’local flatness3 as described in Section 3.3; P=(Px, Py, Pz) contains the control mesh; function returns 1 if stopping criterion <= tol, 0 otherwise */ / * USAGE int difference_l(double *Px, *Py, *Pz) */ int difference_l(); f * stopping criterion ’mesh refinement' as described in Section 3.3; P=(Px, Py, Pz) contains the control mesh; function returns 1 if stopping criterion <= tol, 0 otherwise */ / * USAGE int difference_d(double *Px, *Py, *Pz) * / int difference_d(); / * stopping criterion 'combination criterion’ as described in Section 3.3; P=(Px, Py, Pz) contains the control mesh; 86 function returns 1 if stopping criterion <= tol, 0 otherwise * / /* USAGE int difference_c(double *Px, *Py, *Pz) * / int difference_c(); /* stopping criterion that stops if depth of recursion is greater than given maximum depth of recursion (needed for the implementation of a fixed depth of recursion or the ’global error’ criterion; function returns 1 if stopping criterion <= tol, 0 otherwise */ /* USAGE int difference_f0 * / int difference_f(); /* function that computes the maximum (minimum) distance in the coordinates x, y and z between any two consecutive control points for a given control mesh ’P’; results are written in variables ’Min[]’ and ’MaxG ’ */ / * USAGE void meshmax(double *Px, *Py, *Pz) * / void mesh_max(); 87 /* surf.c */ # include "surf.h" # include "peripher.h" # include <math.h> # include <string.h> /* maximum depth of recursion allowed for the surface approximation in the combination criterion */ # define MAX_RECURSION 10 # define MAXCx.y) (Cx) > (y) ? Cx) : (y)) # define MAX3(x, y, z) (MAX(x,MAX(y,z))) # define MINCx.y) C(y) > 0 ? ((x) > Cy) ? (y) : Cx)) : Cx)) # define MIN3Cx, y, z) CMINCx.MINCy.z))) # define ABSCx) CCx) > 0 ? Cx) : -Cx)) /* macro needed in surfsplit for splitting in v-direction */ #define SPLIT_VCQ, S) forCj=0; j<m; j++) { S+=M-m; *S=Q[M-m]; forCk=l; k<n; k++) ■ C tmp2=*Q; for Ci=0; i<m*Cn-k); i+=m) { tmpl=tmp2; tmp2= CQ Ci]+Q[i+m])/2.0; Q[i]=tmpl; > Q Ci]=tmp2; Q+=m; S-=m; *S=tmp2; 3 * Q-=H-m“l; S++; } • /* macro needed for the local flatness criterion * / # define DISTANCECindex) { tmp[0]=Px[k] - Px[index]; tmp[l]=Py[k] - Py [index] ; tmp[2]=Pz[k] - Pz[index]; t=tmp [0] *D [0] +tmp [1] *D [1] +tmp [2] *D [2] ; if Ctmp [0] *tmp [0] +tmp [1] *tmp [1] + tmp[2]*tmp[2] - t*t/d >= tol) return 0; > /* implementation as described in Section 3.3 */ int difference_lCPx, Py, Pz) double *Px, *Py, *Pz; ■ C double D[3] , tmp [3] ; 88 double d, t; double toll, tol2; int k; tmp[0]=*Px; tmp[l]=*Py; tmp[2]=*Pz; /* check whether surface is locally flat */ /* normal vector */ D [0] = (Py [m-1] -tmp [1] ) * (Pz [M-m] -tmp [2]) - (Pz[m-1]-tmp[2])*(Py[M-m]-tmp[1]); D [1] = (Pz [m-1] -tmp [2] ) * (Px [M-m] -tmp [0] ) - (Px [m-1] -tmp [0] ) * (Pz [M-m] -tmp [2] ); D [2] = (Px [m-1] -tmp [0] ) * (Py [M-m] -tmp [1] ) - (Py[m-1]-tmp[1])*(Px[M-m]-tmp[0]); toll=tol*sqrt(D[0]*D[0]+D[1]*D[1]+D[2]*D[2]); for(k=l; k < M-1; k++) if (ABS(D[0]*(tmp[0] - Px[k]) + D [1] *(tmp[l] - Py[k]) + D[2]*(tmp[2] - Pz[k])) >= toll) return 0; tol2=tol*tol; /* check whether the four boundary lines are locally linear */ D [0] =Px [M-m] -Px [0] ; D[1]=Py[M-m]-Py[0]; D[2]=Pz[M-m]-Pz[0]; d=D [0] *D [0] +D [1] *D [1] +D [2] *D [2] ; for(k=m; k<M-m; k+=m) DISTANCE(0) D [0] =Px [m-1] -Px [M-1] ; D[l]=Py[m-l]-Py[M-l] ; D [2] =Pz [m-1] -Pz [M-1] ; d=D [0] *D [0] +D [1] *D [1] +D [2] *D [2] ; for(k=2*m-l; k<M-l; k+=m) DISTANCE(M-1) D [0] =Px [m-1] -Px [0] ; D [ 1] =Py [m-1] -Py [0] ; D [2] =Pz [m-1] -Pz [0] ; d=D [0] *D [0] +D [1] *D [1] +D [2] *D [2] ; for(k=l; k<m-l; k++) DISTANCE(0) D[0]=Px[M-m]-Px[M-1]; 89 D[l]=Py [M-m]-Py [M-1] ; D[2]=Pz[M-m]-Pz[M“1]; d=D [0] *D [0] +D [1] *D [1] +D [2] *D [2] ; for(k=M-m-l; k<M~l; k++) DISTANCE(M-l) return 1; > /* implementation as described in Section 3.3 */ int difference_d(Px, Py, Pz) double *Px, *Py, *Pz; int i, j; for(i=0; i<n-i; i++) ■ C for(j=0; j<m-i; j++) i if(MAX3(ABS(*Px-*(Px+m)), ABS(*Py-*(Py+m)), ABS(*Pz-*(Pz+m))) >= tol || MAX3(ABS(*Px - *(Px+l)), ABS(*Py - *(Py+l)), ABS(*Pz - *(Pz+i))) >= tol) return 0; Px++; Py++; Pz++; > if(MAX3(ABS(*Px - *(Px+m)), ABS(*Py - *(Py+m)), ABS(*Pz - *(Pz+m))) >= tol) return 0; Px++; Py++; Pz++; > for(j=0; j<m-l; j++) if(MAX3(ABS(*Px - *(Px+l)), ABS(*Py - *(Py+l)), ABS(*Pz - *(Pz+l))) >= tol) return 0; Px++; Py++; Pz++; return 1; > 90 /* self explanatory */ int difference„f() if (rec < tol) return 0; else return 1; > /* implementation as described in Section 3.3 */ int difference_c(Px, Py, Pz) double *Px, *Py, *Pz; ■ C if (rec > MAX_RECURSI0N jI difference_l(Px, Py, Pz) II difference_d(Px, Py, Pz)) return 1; else return 0; > f * ’surfsplit1 algorithm as described in Section 3.3 */ void surfsplit(Q, R, S, T) double *Q, *R, *S, *T; int i, j, k; double tmpl, tmp2; for(j=0; j<n; j++) ■ C R+=m-l; *R=Q[m-l] ; for(k=l; k<m; k++) ■ C tmp2=*q++; for (i=-l; i<m-k-l; i++) ■ C tmpl=tmp2; tmp2=(Q[i]+Q[i+1])/2.0; q[i]=tmpl; > q[i]=tmp2; 91 *— R=tmp2; > q++; R+=m; } q-=M; R-=H; SPLIT_V(q,S) SPLIT_V(R,T) /* self explanatory * / void mesh„max(Px, Py, Pz) double *Px, *Py, *Pz; ■ C int k, j; for(j=0; j<n-1; j++) ■ C for(k=0; k<m-l; k++) Max[0]=MAX3(Max[0], ABS(*Px-+(Px+m)), ABS(*Px-*(Px+l) Hax[l]=MAX3CMax[l], ABS(*Py-*(Py+m)), ABS(*Py-*(Py+l) Max[2]=MAX3(Max[2] , ABS(*Pz-*(Pz+m)), ABS(*Pz-*(Pz+l) Min[0]=MIN3(Min[0], ABS(*Px-*(Px+m)), ABS(*Px-*(Px+i) Min[i]=MIN3(Min[l], ABS(*Py-*(Py+m)), ABS(*Py-*(Py+l) Min[2]=MIN3(Min[2], ABS(*Pz-+(Pz+m)), ABS(*Pz-*(Pz+l) Px++; Py++; Pz++; > Max[0]=MAX(Max[0], ABS C *Px-* (Px+m))); Max[i]=MAX(Max[l] , ABS(*Py-*(Py+m))); Max[2]=MAX(Max[2], ABS(*Pz~*(Pz+m))); Hin[0]=MIN(Min[0], ABS(*Px-*(Px+m))); Min[1]=MIN(Min Ei], ABS(*Py-*(Py+m))); Min[2]=MIN(Min[2], ABS(*Pz-*(Pz+m))); Px++; py++; Pz++; > for(k=0; k<m-l; k++) - C 92 Max[0]: Max[i]: Max[2]: Min[0]= Min[i]: Min[2]: Px++; Py++; Pz++; ’ MAX (Max [0] , : MAX(Max [1] , : MAX (Max [2] , : MIN(Min[0] , : MIN(Min[l] , : MIN(Min [2] , ABSOPx- abs (*Py- ABS(*Pz- ABS(*Px- ABS(*Py- ABSC*Pz- *(Px+i)) *(Py+1)) *(Pz+l)) *(Px+l)) *(Py+l)) *(Pz+l)) /* ’surfgen* algorithm as described in Section 3.3 */ void surfgen(fp, P, func) FILE *fp; double *P; int (*func)(); double *Q; ++rec; if (OfuncXP, P+M, P+2+M)) OUTPUT(P) else ■ c Q= (double *) calloc(9*M, sizeof (double)); surfsplit(P, Q, Q+3*M, Q+6+M) surfsplit(P+M, Q+M, Q+4*M, Q+7+M) surfsplit(P+2*M, Q+2*M, Q+5*M, Q+8+M) surfgen(fp, P, func); surfgen(fp, Q, func); surfgen(fp, Q+3+M, func); surfgen(fp, Q+6*M, func); free(Q); > — rec; 93 /* peripher.h */ /* Header for input and output functions of control meshes (implementation dependent) * / /* Input functions are specifically written to read input from files including the control meshes as provided by Rockwell Int. Corp. */ /* Output written into a file to meet MATLAB specifications to plot the computed surfaces */ /* Macro that prints the control mesh ’P’ into the file with file pointer ’fp’ */ /* function that computes the size (m*n) of the control mesh in file ’infile’; any additional information from file ’infile’ is removed and the sequence of control points is written in the file ’tmpfile’; function returns 0 if no error occured, 1 otherwise * / /* USAGE: int countline(char *infile, char ♦tmpfile) */ int countline(); /* function that reads the control mesh P=(Px, Py, Pz) from file with file pointer ’fp’; function returns number of control points read */ / * USAGE: int input(FILE *fp, double *Px, *Py, *Pz) */ int input(); # define OUTPUT(P) { int k; \ V \ \ for(k=0; k<M; k++) fprintf(fp, "*/,le */,le */,le\n", *(P+k>, *(P+M+k), *(P+2*M+k)); fprintf(fp, "\n"); } 94 /* peripher.c */ # include "peripher.h" # include 1 1 surf.h" # include <stdio.h> # include <string.h> /* implementation of functions as described in "peripher.h" */ int countline(infile, tmpfile) char *infile, *tmpfile; - C FILE *fp; char string[42]; char tmp[L_tmpnam]; tmpnam(tmp); /* write # of columns of control mesh into temporary file * / strcpy(string, "egrep -c ’begin curve ’ "); strcat(string, infile); strcat (string, 1 1 > ") ; strcat(string, tmp); if(system(string) != 0) return 1; / * extract data from inputfile and write it into temporary file for further handling by function ’input’ * / strcpy(string, "egrep -e ’-?[0-9]+\.[0-9]+[eIE-|+[0-9]+]?* ") strcat(string, infile); strcat(string, " > "); strcat(string, tmpfile); if(system(string) != 0) return 1; / * write # of rows of control mesh into temporary file */ strcpy(string, "wc -1 "); strcat(string, tmpfile); strcat(string, " » "); strcat(string, tmp); if(system(string) != 0) return 1; /* read parameters m and n from temporary file * / if ((fp=fopen(tmp, "r")) == NULL) return 1; fscanf(fp, "'/id11, &n) ; fscanf(fp, "’ /,d &M); fclose(fp); remove(trap); if (M'/,n == 0) ■ C m=M/n; return 0; > else return 1; /* read control mesh from temporary file into variable 3PJ * / int input(fp, Px, Py, Pz) FILE *fp; double *Px, *Py, *Pz; • C int ct, count=0; while((ct=fscanf Cfp, 1 1 */»le '/le '/,le\n", Px++, Py++, Pz++))>0) count+=ct; return count; > 96 I * data.h * / / * Header file with functions for numerical experiments */ /* variables for numerical experiments * / int Max.rec; /* maximum depth of recursion for stopping the iteration */ int Min_rec; /* minimum depth of recursion for stopping the iteration */ int Out_count; /* number of calls of output macro * / / * function to generate a Bezier surface as described in Section 3.3; but gathers additional information as needed for the numerical experiments described in Chapter 4 */ / * USAGE void surfgan(FILE *fp, double *P, int (*int)() stopping_criterion) * / void test_surfgen(); /* function that prints the data obtained by the numerical experiments (for a sample of the output see Appendix B.3); 'flag* is the flag for the given stopping criterion * / / * USAGE void output.comp(char *infile, flag) */ void output_comp(); 97 /* data.c */ # include "data.h" # include "surf.h" # include "peripher.h" it include <math.h> # include <string.h> # define MAX(x,y) CCx) > (y) ? Cx) : (y)) # define MAX3(x, y , z) CMAX(x,MAX(y,z))) # define MIN(x,y) CCy) >0 ? CCx) > (y) ? (y) : (x)) : (x)) # define MIN3(x, y, z) (MINCx,HIN(y,z))> # define ABS(x) CCx) > 0 ? Cx) : -Cx)) void output_corapCinfile, flag) char *infile, flag; ■ C */,s\n", infile); printf C'DATA ANALYSIS\nM); printf C"--------- *---\n") ; printfC"control mesh: switchCflag) { case ’1J: printfC"\nstopping criterion: printfC"given tolerance: break; case *dJ: printfC"\nstopping criterion: mesh refinement\n"); printfC"given tolerance: break; case ’cJ: printfC"\nstopping criterion: combination criterion\n"); printf Cgiven tolerance: break; case ’m’: printfC"\nstopping criterion: printf Cgiven absolute error: '/tg\n" , tol) ; break; case ’f' : printfC\nstopping criterion: depth of recursion\n"); printf Cgiven depth of recursion: '/.iXn", Cint) tol); break; > printf C\nnumber of control points: */,5d **/,4d = '/d\n", m, n, M) ; printf C"number of output points: '/,5i *‘ /,4d = '/,g\n" , Qut_count, M, Cfloat) 0ut_count*M); local flatness\n"); */.g\n", tol); */,g\n", tol); combination < ’ /,g\n", tol); absolute error\n"); 98 if (flag != ’m’) - C printf ("\nmax. depth of recursion: '/,d\n", Max_rec); printf ("rain, depth of recursion: '/,d\nn, Min_rec); > else printf("\ndepth of recursion: printf("\nmax. difference: printf(" printf(" printf("\nmin. difference: printf(" printf(" printf (1 1 \nELAPSED EUNTIME\n") ; in x: y.g\n", Max[0]) in y: XgVi", Max[l]) in z: '/,g\nM , Max [2] ) in x: */.g\n", Min[0]) in y: '/.g\n", Min[l]) in z: '/.g\n", Min [2]) '/,d\nM, Max_rec) ; printf ("--------------- \n"); \n"); > /* ’surfgen’ algorithm as described in Section 3.3 (with additional information for numerical experiments) */ void test_surfgen(fp, P, func) double *P; int (*func)(); { double *q; ++rec; if ((*func)(P, P+M, P+2+M)) OUTPUT(P) mesh_max(P, P+M, P+2*M); Max_rec=MAX(Max_rec, rec); Min_rec=MIN(Min_rec, rec); Out_count++; > else • c Q= (double *) calloc(9*M, sizeof (double)); surfsplit(P, Q, q+3*M, Q+6+M); surfsplit(P+M, Q+M, q+4*M, Q+7+M); surfsplit(P+2*M, q+2*M, q+5*M, q+8*M); surfgen(fp, P, func); surfgen(fp, Q, func); surfgen(fp, q+3*M, func); surfgen(fp, Q+6*M, func); free(Q); /* main.c */ /* program for Bezier surface generation */ /* USAGE: main [ TDL | -d[TOL] I -c[TOL: | -m[TDL] I -f [TOL]] [infile [outfile]] */ /* the flags stand for the different stopping criteria: -d ’mesh refinement’ -c ’combination criterion’ -m ’global error’ -f ’maximum depth of recursion’ default is ’local flatness’ in addition a specification of a tolerance for the stopping criterion is optional as well as of input file and output file */ # include "surf.h" # include "peripher.h" # include "data.h" # include <stdio.h> # include <stdlib.h> # include <string.h> # include <math.h> # define MAX_LE # define T0L1 # define T0L2 # define MAX_REC # define IN_FILE "/dev/tty" # define OUT_FILE "/dev/null # define MAX(x,y) C(x) > (y) # define MAX3(x,y,z) (MAX(MAX(x 24 / * maximum length for file names * / 0.0001 /* default tolerance for stopping criterion ’local flatness’ */ 0.001 / * default tolerance for ’mesh refinement’ and ’combination criterion’ * / maximum depth of recursion for ’global error’ and ’maximum depth of recursion’ */ /* default inputfile */ " /* default outputfile */ ? Cx) : (y)) ,y),(z))) /* mainCargc, argv) int argc; char **argv; file pointer for input and output file */ file name for input file * / file name for output file */ file name for temporary file * / variable for control mesh * / / * variable for stopping criterion (default is ’local flatness’) * / variable for flag ’stopping criterion’ */ FILE *fp; /* char infile[MAX_LE]; / * char outfile[MAX_LE]; /* char tmpfile[L_tmpnam]; /* double *P; /* int (*func)()=difference. -1; char flag=’l’; /* 101 double error; int i; /* variable for tolerance in ’global error’ criterion */ /* counter variable */ /* initialization */ tol=TOLl; for(i=0; i<3; i++) Min[i]=i.0e+08; Max[i]=0.0; } /* handling of command line arguments */ if(— argc >0) if(**++argv == ’-’) ■ C switch(*++*argv) • { case ’d*: /* stopping criterion flag=’d’; ’mesh refinement’ */ func=difference_d; if C(tol=atof(++*argv++)) <= 0) tol=T0L2; break; case ’c’: /* stopping criterion flag=’c’; 'combination criterion’ */ func=difference_c; if ((tol=atof(++*argv++)) <= 0) tol=T0L2; break; case ’m’: / * stopping criterion flag=’m’; ’global error’ */ func=difference_f; if ((error=atof(++*argv++)) <= 0) tol=MAX_REC; else flag='o’; break; case ’f ’ : /* stopping criterion flag=’f’; 'max. depth of recursion’ * / func=difference_f; if C(tol=atof(++*argv++)) <= 0) tol=MAX_REC; break; default: perrorC’ USAGE: main [ TOL I -d[T0L] I -c[T0L] | -m[TOL] I 102 -f[TOL]] [infile [outfile]]"); exit(); break; > argc— ; > else ifCisdigit(**argv)) /* read tolerance for default stopping criterion ’local flatness’ */ tol=atof(*argv); argv++; argc— ; } f * inputfile and outputfile handling */ switch(argc) - [ case 2: strncpyCinfile, *argv++, MAX_LE); strncpy(outfile, *argv, MAX_LE); break; case i: strncpyCinfile, *argv, HAX_LE); strncpy(outfile, OUT.FILE, HAX_LE); break; case 0: strncpyCinfile, IN_FILE, MAX_LE); strncpy(outfile, QUT„FILE, MAX_LE); break; default: perror("USAGE: main [ TOL I -d[T0L] I -c[T0L] | -m[TDL] | -f[T0L]] [infile [outfile]]"); exitO ; > /* reading the the control mesh out of the input file into the ’P’ * / tmpnam(tmpfile); if (countline(infile, tmpfile)) ■ C remove(tmpfile); perror("Data of inputfile has not correct size"); exit() ; } if((fp=fopen(tmpfile, "r")) == NULL) remove(tmpfile); perror("Inputfile could not be opened"); 103 exitO; > P= (double *) calloc(3*M, sizeof (double)); if (input(fp, P, P+M, P+2*M) != 3*M) £ free(P); fclose(fp); remove(tmpfile); perror("Data of inputfile could not be read correctly'1); exit(); > fclose(fp); /* computing the maximum depth of recursion for the ’global error’ criterion * / if(flag == ’o’) -c mesh_max(P, P+M, P+2+M); if ((tol= (ceil((log(MAX3(Max[0] , Max[l!l, Max[2])) -log(error)))/log(2.0))) <= 0) free(P); fclose(fp) ; remove(tmpfile); perror("Given max. error is greater than max. error of input data"); exitO ; > for(i=0; i<3; i++) - C Min[i]=1.0e+08; Max [i]=0.0; > f l a g ^ m ’ ; > /* open the output file * / if((fp=fopen(outfile, "w")) == NULL) free(P); remove(tmpfile); perror("Outputfile could not be opened"); exitO; > 104 /* generation of the Bezier surface * / surfgen(fp, P, func); /* generation of the Bezier surface for numerical experiments /* test_surfgen(fp, P, func); if (flag == 5m’) tol=error; output_comp(infile, flag); */ I * cleaning up * / free(P); fclose(fp); remove(tmpfile); /* makefile for main.c * / main: main.o surf.o peripher.o data.o cc -p -dalign -02 -o main main.o surf.o peripher.o data.o -lm peripher.o: peripher.h cc -p -c peripher.c surf.o: surf.h peripher.h cc -p -c surf.c data.o: data.h surf.h peripher.h cc -p ~c data.c main.o: surf.h peripher.h data.h cc -p -c main.c clean: rm -f main 106 A.3 Computation of the Machine Precision /* program to compute the machine precision */ main() • c double eps=1.0; while(eps+l.0 F= 1.0) eps/=2.0; printf("\nmachine precision is eps = '/,lg\n\n", eps*2.0); > /* answer */ machine precision is eps = 2.22045e-i6 107 Appendix B Test Data B .l Plots of the Control Meshes control mesh 1 6 5- 4- 6 5- 4 - -2.5 1 1.4 1.2 -3 1 control mesh 4 - 2.5 „ V \ •N> i 1.4 1.2 CO ( O S S ® Jb A -3$ • * « * * * ’ * * , \ , * * » « , * • * * , * • * * , \ control mesh 3 10 5 0 0 1.5 0.5 0 control mesh 6 10 5 0 0 1.5 0.5 0 110 B.2 Comprehensive Test Results for the Stopping Criteria This section contains a summary of the most important numerical results for the three control meshes mesh 1, mesh 2 and mesh 3 as discussed in Section 4.1 control mesh: mesh 1 stopping criterion: local flatness # of control points: 50 X 11 — 550 criterion tolerance 0.1 0.05 0.01 0.005 0.001 0.0005 0.0001 # of out put points 2,200 5,500 31,900 48,400 299,200 724,900 3,208,150 max. depth of recursion 1 2 4 4 6 7 8 min. depth of recursion 1 1 2 3 4 4 6 runtime (in sec.) 3 5 28 44 262 639 2,871 # of calls surfsplit 3 9 57 87 543 1,317 5,832 msec/call surfsplit 26.7 26.7 24.4 25.5 24.0 23.9 23.9 # of calls surfgen 5 13 77 117 725 1,757 7,777 msec/call surfgen 0.0 0.0 0.3 0.5 0.6 0.7 0.7 111 control mesh: mesh 2 stopping criterion: local flatness # of control points: 75 x 11 = 825 criterion tolerance 0.1 0.05 0.01 0.005 0.001 0.0005 0.0001 # of out put points 8,250 23,100 124,575 275,550 1,243,280 2,364,450 11,858,600 max. depth of recursion 2 3 5 6 8 8 9 min. depth of recursion 1 2 3 3 4 5 6 runtime (in sec.) 8 21 112 311 1,108 2,604 10,580 # of calls surfsplit 9 27 150 333 1,506 2,865 14,373 msec/call surfsplit 50.0 50.8 48.7 48.1 47.4 47.6 47.4 # of calls surfgen 13 37 201 445 2,009 3,821 19,165 msec/call surfgen 0.8 0.8 0.7 2.4 1.2 2.2 0.9 control mesh: mesh 3 stopping criterion: local flatness # of control points: 100 X 11 = 1,100 ci'iterion tolerance 0.1 0.05 0.01 0.005 0.001 0.0005 0.0001 # of out put points 30,800 60,500 436,700 941,600 3,782,900 9,171,800 45,844,700 max. depth of recursion 3 3 6 7 8 9 10 min. depth of recursion 2 2 3 4 5 6 7 runtime (in sec.) 36 55 396 848 3,374 8,234 40,812 # of calls surfsplit 27 54 396 855 3,438 8,337 41,676 msec/call surfsplit 78.6 81 78.3 78.8 79.3 78.8 78.3 # of calls surfgen 37 73 529 1,141 4,585 11,117 55,569 msec/call surfgen 1.1 1.0 1.4 1.6 1.3 1.5 1.3 112 control mesh: mesh 1 stopping criterion: mesh refinement # of control points: 50 x 11 = 550 criterion tolerance 0.1 0.05 0.01 0.005 0.001 0.0005 0.0001 # of out put points 550 2,200 36,850 144,100 4,013,350 15,898,300 576,717,000 max. depth of recursion 0 1 4 5 7 8 10 min. depth of recursion 0 1 3 4 6 7 10 runtime (in sec.) 1 3 33 129 3,817 14,014 508,669 # of calls surfsplit 0 3 66 261 7,296 28,905 msec/call surfsplit _ 23.3 24.1 25.0 24.0 24.3 „ # of calls surfgen 1 5 89 349 9,729 38,541 msec/call surfgen 0.0 0.0 0.5 0.7 0.9 0.7 _ control mesh: mesh 2 stopping criterion: mesh refinement # of control points: 75 x 11 = 825 criterion tolerance 0.1 0.05 0.01 0.005 0.001 0.0005 0.0001 # of out put points 825 3,300 109,725 409,200 9,274,650 36,910,500 865,075,000 max. depth recursion 0 1 4 5 7 7 10 min. depth of recursion 0 1 3 4 6 6 10 runtime (in sec.) 1 3 106 373 8,504 33,566 1,055,048 # of calls surfsplit 0 3 132 495 11,241 44,739 msec/call surfsplit . 43.3 47.1 47.9 47.5 47.7 # of calls surfgen 1 5 177 661 14,989 59,653 msec/call surfgen 0.0 0.0 1.2 1.1 1.0 1.3 113 control mesh mesh S stopping criterion: mesh refinement # of control points: 100 x 11 = 1100 criterion tolerance 0.1 0.05 0.01 0.005 0.001 0.0005 0.0001 # of out put points 4,400 11,000 192,500 700,700 14,445,200 57,539,990 1,341- 106 max. depth of recursion 1 2 4 5 7 8 11 min. depth of recursion 1 1 3 4 6 7 10 runtime (in sec.) 6 11 178 641 13,262 53,901 1,213,605 # of calls surfsplit 3 9 174 636 13,131 52,308 msec/call surfsplit 83.3 68.9 77.7 78.6 78.2 66.5 # of calls surfgen 5 13 233 849 17,509 69,745 msec/call surfgen 4.0 1.5 1.0 1.6 1.3 1.7 control mesh mesh 1 stopping criterion: combination criterion # of control points: 50 x 11 = 550 criterion tolerance 0.1 0.05 0.01 0.005 0.001 0.0005 0.0001 # of out put points 550 2,200 26,950 48,400 299,200 724,900 3,208,150 max. depth of recursion 0 1 4 4 6 7 8 min. depth of recursion 0 1 2 3 4 4 5 runtime (in sec.) 1 2 25 43 269 642 2838 # of calls surfsplit 0 3 48 87 543 1,317 5,832 msec/call surfsplit _ 23.3 24:6 24.6 24.9 24.6 24.1 # of calls surfgen 1 5 65 117 725 1,757 7,777 msec/call surfgen 0.0 0.0 0.8 0.6 0.7 0.7 0.6 114 control mesh: mesh 2 stopping criterion: combination criterion # of control points: 75 x 11 = 825 criterion tolerance 0.1 0.05 0.01 0.005 0.001 0.0005 0.0001 # of out put points 825 3,300 84,975 196,350 1,233,380 2,364,450 11,858,600 max. depth of recursion 0 1 4 5 7 8 9 min. depth of recursion 0 1 3 3 4 5 6 runtime (in sec.) 1 4 79 182 1,112 2,110 11,564 # of calls surfsplit 0 3 102 237 1494 2,865 14,373 msec/call surfsplit 53.3 47.7 49.3 47.7 47.2 47.5 # of calls surfgen 1 5 137 317 1,993 3,821 19,165 msec/call surfgen 0.0 0.0 0.7 0.9 0.8 0.8 0.9 control mesh: mesh 3 stopping criterion: combination criterion # of control points: 100 x 11 = 1,100 criterion tolerance 0.1 0.05 0.01 0.005 0.001 0.0005 0.0001 # of out put points 825 3,300 84,975 196,350 1,233,380 2,364,450 11,858,600 max. depth of recursion 1 2 4 5 6 7 10 min. depth of recursion 1 1 3 4 4 5 7 runtime (in sec.) 1 4 79 182 1,112 2,110 11,564 # of calls surfsplit 0 3 102 237 1494 2,865 14,373 msec/call surfsplit 53.3 47.7 49.3 47.7 47.2 47.5 # of calls surfgen 1 5 137 317 1,993 3,821 19,165 msec/call surfgen 2.0 3.1 1.2 1.4 1.8 1.7 1.3 115 B.3 Samples of Numerical Test Results This section contains samples of the test data obtained by the program in Appendix A. DATA ANALYSIS control mesh: stopping criterion: given tolerance: mesh 3 local flatness 0.1 number of control points: number of output points: 100 * 11 = 1100 28 *1100 = 30800 max. depth of recursion: min. depth of recursion: 3 2 max. difference: in x: 0.0224518 in y: 0.0198772 in z: 0.016889 min. difference: in x: 3.8147e-10 in y: 2.54745e-09 in z: 2.84405e-07 ELAPSED RUNTIME real 2:46.4 user 36.3 sys 1.1 RUNTIME ANALYSIS '/.time cumsecs #call ms/cal name 5.8 24.01 27 78.52 _surfsplit 2.5 28.09 29 31.72 _mesh_max 0.4 35.22 37 4.05 _diff erence_l 0.3 35.72 30828 0.00 _fprintf 0.1 36.45 37 1.08 _surfgen 0.0 36.77 1 0.00 _countline 0.0 36.77 18 0.00 _free 0.0 36.77 1 0.00 _input 0.0 36.77 1 0.00 _output„comp 116 DATA ANALYSIS control mesh: stopping criterion: given tolerance: number of control points: number of output points: max. depth of recursion: min. depth of recursion: max. difference: in x: in y: in z: min. difference: in x: in y: in z: ELAPSED RUNTIME real 6:38.8 user 55.0 sys 1.3 RUNTIME ANALYSIS '/.time cumsecs #call ms/cal name 7.9 22.89 54 80.93 _surfsplit 3.2 40.76 56 31.25 _mesh_max 0.3 53.96 73 2.60 _differencs_l 0.3 54.50 60555 0.00 _fprintf 0.1 54.84 73 0.96 _surfgen 0.0 55.28 1 0.00 _countline 0.0 55.28 27 0.00 _free 0.0 55.28 1 0.00 _input 0.0 55.28 24 0.00 _malloc 0.0 55.28 1 0.00 _output.comp mesh 3 local flatness 0.05 100 * 11 = 1100 55 *1100 = 60500 3 2 0.0221959 0.0139245 0.0141005 7.45182e-13 2.54745e-09 2.84405e-07 DATA ANALYSIS control mesh: stopping criterion: given tolerance: mesh 3 local flatness 0 .0 1 number of control points: number of output points: 100 * 11 397 *1100 1100 436700 max. depth of recursion: min. depth of recursion: 6 3 max. difference: m x: in y: in z: 0.00988732 0.00639739 0.00667467 min. difference: m x: in y: in z: 1.33227e-15 9.68495e-10 6.98195e-09 ELAPSED RUNTIME real 43:37.9 user 6:36.1 sys 6.6 RUNTIME ANALYSIS Xtime cumsecs: #call ms/cal name 7.8 161.64 396 78.33 _surfsplit 3.1 291.87 398 31.56 _mesh_max 0.5 385.70 437097 0.00 _fprintf 0.4 394.12 529 2.76 _difference_l 0.2 398.97 529 1.36 _surfgen 0.0 400.19 1 10.00 _input 0.0 400.20 1 0.00 _countline 0.0 400.20 144 0.00 _free 0.0 400.20 138 0.00 _malloc 0.0 400.20 1 0.00 _output_comp DATA ANALYSIS control mesh: stopping criterion: given tolerance: mesh 3 local flatness 0.005 number of control points: number of output points: 100 * 11 856 *1100 1100 941600 max. depth of recursion: min. depth of recursion: 7 4 max. difference: in x: in y: in z: 0.00561025 0.00377534 0.00404464 min. difference: m x: in y: in z: 1.33227e-15 2.86949e-10 6.98195e-09 ELAPSED RUNTIME real 1:35:09.1 user 14:08.3 sys 17.5 RUNTIME ANALYSIS '/.time cumsecs #call ms/cal name 7.9 353.66 855 78.80 _surfsplit 3.3 634.41 857 32.67 _mesh_max 0.4 835.71 942456 0.00 _fprintf 0.4 852.08 1141 2.66 _difference_l 0.2 853.93 1141 1.62 _surfgen 0.0 857.67 291 0.03 „malloc 0.0 857.70 1 0.00 _countline 0.0 857.70 298 0.00 _f ree 0.0 857.70 1 0.00 _input 0.0 857.70 1 0.00 _output_comp DATA ANALYSIS control mesh: stopping criterion: given tolerance: mesh 3 local flatness 0.001 number of control points: number of output points: 100 * 11 = 1100 3439 *1100 = 3.7829e+06 max. depth of recursion: min. depth of recursion: 8 S max. difference: in x: in y: in z: 0.00280208 0.00166624 0.00190732 min. difference: m x: in y: in z: 4.44089s-16 4.21463e-ll 4.51479e-ll ELAPSED RUNTIME real 4:13:59.7 user 56:13 .9 sys 1:06 .9 RUNTIME ANALYSIS */,time cumsecs ttcall ms/cal name 8.0 1418.06 3438 79.27 _surfsplit 3.2 2420.77 3440 32.21 _mesh_max 0.4 3336.063786339 0.00 _fprintf 0.4 3390.40 4585 2.70 _difference_l 0.2 3408.62 4585 1.29 _surfgen 0.0 3421.32 1160 0.02 _fres 0.0 3421.43 1 0.00 „countline 0.0 3421.43 1 0.00 _input 0.0 3421.43 1152 0.00 _malloc 0.0 3421.43 1 ' 0.00 „output_comp DATA ANALYSIS control mesh: stopping criterion: given tolerance: mesh 3 local flatness 0.0005 number of control points: number of output points: 100 * 11 = 1100 8338 *1100 = 9.1718e+06 max. depth of recursion: min. depth of recursion: 9 6 max. difference: in x: in y: in z: 0.0014025 0.00132548 0.00106269 min. difference: m x: in y: in z: 2.22045e-16 4.27036e-12 6.52611e-12 ELAPSED RUNTIME real 15:52:37.8 user 2:17:14.3 sys 2:20.3 RUNTIME ANALYSIS */,time cumsecs #call ms/cal name 7.9 3427.69 8337 . 78.81 „surfsplit 3.2 5878.07 8339 32.19 _mesh_max 0.4 8142.699180138 0.00 _fprintf 0.4 8238.95 11117 2.80 _difference_l 0.2 8285.94 11117 1.52 _surfgen 0.0 8318.57 2785 0.01 _malloc 0.0 8318.70 2794 0.00 _f ree 0.0 8318.71 1 10.00 _input 0.0 8318.76 1 0.00 _countline 0.0 8318.76 1 0.00 _output_comp DATA ANALYSIS control mesh: stopping criterion: given tolerance: number of control points: number of output points: max. depth of recursion: min. depth of recursion: max. difference: in x: in y: in z: mesh 3 local flatness 0.0001 100 * 11 = 1100 41677 *1100 = 4.58447e+07 10 7 0.000701251 0.000504002 0.000520102 min. difference: in x: 2.22045e-16 in y: 1.15019e-13 in z: 8.94113e-12 ELAPSED RUNTIME real 51:45:38.5 user 11:20:12.4 sys 14:45.7 RUNTIME ANALYSIS '/.time cumsecs #call ms/cal name 7.9 17105.19 41676 78.33 _surfsplit 3.2 29269.68 41678 31.96 _mesh_max 0.4 40567.5945886377 0.00 _fprintf 0.4 41047.39 55569 2.80 _difference_l 0.2 41351.48 55569 1.33 _surfgen 0.0 41430.25 13898 0.01 _malloc 0.0 41430.38 13908 0.01 _free 0.0 41430.76 1 0.00 _countline 0.0 41430.76 1 0.00 _input 0.0 41430.76 1 0.00 .output_comp INFORMATION TO USERS This manuscript has been reproduced from the microfilm master. UMI films the text directly from the original or copy submitted. Thus, some thesis and dissertation copies are in typewriter face, while others may be from any type of computer printer. The quality of this reproduction is dependent upon the quality of the copy submitted. Broken or indistinct print, colored or poor quality illustrations and photographs, print bleedthrough, substandard margins, and improper alignment can adversely affect reproduction. In the unlikely event that the author did not send UMI a complete manuscript and there are missing pages, these will be noted. Also, if unauthorized copyright material had to be removed, a note will indicate the deletion. Oversize materials (e.g., maps, drawings, charts) are reproduced by sectioning the original, beginning at the upper left-hand comer and continuing from left to right in equal sections with small overlaps. Each original is also photographed in one exposure and is included in reduced form at the back of the book. Photographs included in the original manuscript have been reproduced xerographically in this copy. Higher quality 6” x 9" black and white photographic prints are available for any photographs or illustrations appearing in this copy for an additional charge. Contact UMI directly to order. A Bell & Howell Information Company 300 North Z eeb Road. Ann Arbor. M l 48106-1346 USA 313/761-4700 800/521-0600 UMI Number: 1376486 UMI Microform 1376486 Copyright 1995, by UMI Company. All rights reserved. This microform edition is protected against unauthorized copying under Title 17, United States Code. UMI 300 North Zeeb Road Ann Arbor, MI 48103
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
A finite element approach on singularly perturbed elliptic equations
PDF
Independent process approximation for the coupon collector's problem
PDF
A Fuzzy Controlled Nonuniform Diffusion Model For Competing Products
PDF
Complementarity problems over matrix cones in systems and control theory
PDF
A model reference adaptive control scheme for linear delay systems
PDF
The analysis of circular data
PDF
Trifluoromethanesulfonates (triflates) for organic syntheses
PDF
Using The Cayley-Hamilton Theorem To Solve Linear Systems Of Difference Equations
PDF
If only I could remember
PDF
Adaptive subdivision algorithm in postprocessing of compressed three-dimensional graphic data
PDF
Altered interaction of human endothelial cells to the glycosylated laminin
PDF
Structural equation modeling in educational psychology
PDF
An analysis of nonresponse in a sample of Americans 70 years of age and older in the longitudinal study on aging 1984-1990
PDF
The study of temporal variation of coda Q⁻¹ and scaling law of seismic spectrum associated with the 1992 Landers Earthquake sequence
PDF
An improved iterative procedure for system identification with associative memory matrices
PDF
Pulse oximetry failure rates
PDF
Structural geology of the Chiwaukum schist, Mount Stuart region, central Cascades, Washington
PDF
A physiologic model of granulopoiesis
PDF
Diminishing worlds?: the impact of HIV/AIDS on the geography of daily life
PDF
The effects of dependence among sites in phylogeny reconstruction
Asset Metadata
Creator
Molt, Christian
(author)
Core Title
Implementation aspects of Bézier surfaces to model complex geometric forms
School
Graduate School
Degree
Master of Science
Degree Program
Applied Mathematics
Degree Conferral Date
1994-12
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
mathematics,OAI-PMH Harvest
Language
English
Contributor
Digitized by ProQuest
(provenance)
Advisor
Proskurowski, Wlodek (
committee chair
), Blum, Edward K. (
committee member
), Domaradzki, Julian A. (
committee member
)
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c18-3357
Unique identifier
UC11356858
Identifier
1376486.pdf (filename),usctheses-c18-3357 (legacy record id)
Legacy Identifier
1376486-0.pdf
Dmrecord
3357
Document Type
Thesis
Rights
Molt, Christian
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the au...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus, Los Angeles, California 90089, USA
Tags
mathematics