ll Milli T115 11111111le "\l‘ficifl 31293 0109 7"“- “\w-I-t-I— -_.._ , gun _. r- c \ 1 J .7 c. .. ...- E. ,i' '.' o . 3‘ o .. ,- fl 4 ..~. I '- . .. .5L’...’».-.3r.r - 9;..3 biu-‘i h. 0 0.1 Unub‘ Wflwr “’ ‘vK-ln" This is to certify that the dissertation entitled THE DES IGN 0F IATION AL 3- SP1. IN E ALGOR ITBIS INTERACTIVE COLORFggADING or SURFACES presented by Mark Norman Pickelnnnn has been accepted towards fulfillment of the requirements for Doctor of Philoggflu degreeinw Engineering MAS 04 8 Major professor Date 4/M/KS Eric D. Goodman MS U i: an Affirmative Action/Equal Opportunity Institution 0- 12771 MSU RETURNING MATERIALS: Place in book drop to remove this checkout from LIBRARIES . 532:. your record. FINES W11] be charged if book is returned after the date stamped below. ‘WJ‘J‘Z’V‘ I 111 _‘ ,. - t? 2.1 it. 3 tr. via"...- ‘- ' l." ”My ,7 I." 7 ’ v‘Wtc’vTSix 2.) v‘ A ,§’ i v A .3 'I It a", THE DESIGN OF RATIONAL B-SPLINE ALGORITHMS FOR INTERACTIVE COLOR SHADING OF SURFACES By Hark Norman Pickelmann A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mechanical Engineering 1985 THE DESIGN OF RATIONAL H-SPLINE ALGORITHIS FOR INTERACTIVE COLOR SHADING OF Submitted by: lark N. Pickelmanm Ph.D. Candidate Approved by: Erik D BoodCan Proteeeor Diaeertation Adviaor n R. Lloyd hairperaon Department of Ieehanical Engineering ii SURFACES 4—13-34: Date Mat/ta Dete .é/M Date ABSTRACT THE DESIGN OF RATIONAL H-SPLINE ALGORITHIS FOR INTERACTIVE COLOR SHADING 0F SURFACES By Hark Norman Pickelmann This dissertation presents several algorithms developed for use with CAD/CAI systems. The new algorithms allow for more efficient evaluations of the entire range of rational R-spline curves and surfaces. A class of rational B-spline called Enhanced Uniform is defined. The algorithms developed include an algorithm which is used es s preprocessor to transform the definition of a nonuniform rational B-spline surface to an equivalent surface based on the Enhanced Uniform Rational B-spline. The second algorithm deve10ped is used for evaluation of enhanced uniform B-splines and their derivatives. Numerical considerations in the implementation of these algorithms are discussed. Examples of the use of these algorithms to generate color images in a surface assessment program are included. To my loving wife Kristi. For my sons Daniel and Kevin and any others who come along. iii AKNO'LEDGEHENTS The author would like to thank my committee of Dr. Erik D. Goodman. Dr. James E. Bernard. Dr. Ronald C. Rosenberg. and Dr. Jake M. Plotkin for all their help getting me through. A special thanks goes to Dr. James E. Bernard. for his help. friendship. and advice. Although he left during the course of this work he is responsible for starting me in the Ph.D. pragram and he was there to get me through qualifers when I needed him the most. The author would like to thank the General Dynamics Corporation for their funding of the COLORSCOPE project. I would like to thank my family which includes mom and dad. my brother. my mother and father in law. and my brothers and sisters in law for all their encouragement during the course of my studies. Finally. I would like to thank my friends at MSU. for the classes we took together and for hanging out on friday afternoons. but most of all I would just like to thank them for being my friends. iv TABLE OF CONTENTS LIST OF FIGURES ......COOOOOOCOO0.00.0000... LIST OF TABLES ............OOOOOOOOOOOCOOOOO CHAPTER I --- INTRODUCTION 1.0 INTRODUCTION ..................... 1.1 ERROR CHECKING ................... 1.2 SURFACE NODELS ................... 1.2.1 The Rational B-spline Surface 1.3 OVERVIEW OF THE DISSERTATION ..... CHAPTER II --- H-SPLINE CURVES AND SURFACES 2.03‘SPL1NEREVIE' eeeeeeeeeeeeeeeeeeeeee 2.1 B-SPLINE SPACE CURVE ............. 2.1.1 Basis Functions .............. 2.1.2 Rational B-spline Space Curves 2.2 B-SPLINE SURFACE ................. 2.2.1 Rational B-spline Surfaces ... 2.3 [NOT VECTOR RESTRICTIONS ......... 2.3.1 Uniform [not Vectors ......... 2.3.2 Enhanced Uniform Knot Vectors ix xi 19 22 24 24 25 26 CHAPTER III --- CURHS ALGORITHM 3.0 ENHANCED UNIFORI RATIONAIJ B’SPLINB eeeeeeeeeeee 3.1 ADVANTAGES OF UNIFORN KNOT VECTOR ............. 3.1.1 Knot Vector Translations .................. 3.1.2 Enhanced Uniform Knot Vector Compression .. 3.1.3 Knot Vector Compression Algorithm ......... 3.2 THE EVALUATION ALGORITHN . 3.3 LINITATIONS OF ENHANCED UNIFORN KNOT VECTOR ... CHAPTER IV --- CURRS ALGORITHN ANALYSIS 4.0 ALGORITRN COIPARISON ..... 4.1 COX-DEBOOR ALGORITHN ..... 4.2 BIG-0H ANALYSIS eeeeeeeeee 4.2.1 Operation Analysis ... CHAPTER V --- CONVERSION ALGORITHN 5.0 NURHS-TO-EURBS CONVERSION 5.1 INITIAL EXPLORATION OF THE CONVERSION PROHLEN . 5.1.1 Knot Vector Rescaling 5.1.1.1 An Example Of Approximate Knot Vector 5.2 EXACT NURHS-TO-EURBS CONVERSION ............... 5.2.1 Bxp‘ndin‘ Th. lnot vectot ......OOOCOOOOOOC 5.2.2 Conversion By Knot Vector Expansion ....... 3.2.3 Storage Considerations vi 28 29 31 33 33 39 50 51 52 CHAPTER VI--- IMPLEMENTATION CONSIDERATIONS 6.0 NUMERICAL CONSIDERATIONS ...................... 65 6.1 BUILDING OF THE BLENDING FUNCTION MATRICES .... 67 6.2 EVALUATION OF THE COX-DEBOOR AND CURHS ALGORITHMS .................................... 68 6.3 SIMULTANEOUS EQUATION SOLVER .................. 69 6.4 EXAMPLES OF EVALUATION PROBLEMS ............... 70 6.4.1 Evaluation of Normal Components ........... 70 6.4.2 Normalization Of The Normal ............... 73 6.4.3 F‘c. To p.00 OO.......OCOOOOOOCOOOOOOOOO... 75 CHAPTER VII --- EXAMPLES AND CONCLUSIONS 7.0 BXA‘PLES AND CONCLUSIONS 0 O O O O O O O O O O O O O O O O O O O O O 18 1.1 Ex.'pl°s 0.000.000...00.000.000.00.00......0... 89 7.2 concln‘ions ..........OIOCOOOOOOOCOO00.0.0.0... 84 APPENDIX A --- CALCULATION EXAMPLE A. RATIONAL B-SPLINE CALCULATION EXAMPLE .......... 87 APPENDIX B --- MATRIX FORMULATION B. FORMULATION OF BLENDING FUNCTION MATRIX ........ 97 APPENDIX C --- KNOT VECTOR TRANSLATION C. EFFECTS OF KNOT VECTOR TRANSLATION ON BLENDING 104 vii APPENDIX D --- [NOT VECTOR COMPRESSION ALGORITHM D.O THE COMPRESSION ALGORITHM .................... D.1 Example of Knot Vector oompreaaion ........... APPENDIX E --- NURES TO EURBS CONVERSION E. EXAMPLE OF NURHS TO EURBS CONVERSION .......... APPENDIX F --- CONVERSION ALGORITHMS F. CODE FOR THE NURBS TO EURBS CONVERSION ........ LIST OF REFERENCES viii 110 112 115 121 135 Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure LIST OF FIGURES Blending Functions of Order Up to Four Fourth Order Blending Functions Non-Zero Blending Functions for i=6 .. Fourth Order Blending Functions. Knot SOVOII EQUII to ant Sileaeeeeeeeeeaee Fourth Order Blending Functions. Knot Sir. Seven. and Eight Equal Fourth Order Blending Functions ...... Homogeneous Coordinates vs Parametric V‘ti‘ble ..........OOOCOOIOOOOOOOO.... Three-space Plot of Homogeneous Coordinates Plot of X-Y Curve Uniform vs Nonuniform Curves ......... Uniform Fit Using End Slopes Uniform Fit Using First Subsegment ... Uniform Fit Using Least Squsres ...... Uniform Fit By [not Vector Expansion . Screen Space Coordinate System Single Precision Evaluation Double Precision Evaluation .......... Unnormalized 2 Component Calculation . Normalized 2 Component Calculation ... ix 10 12 13 14 15 16 2O 21 22 55 56 57 58 6O 67 73 74 75 Figure 7.1.1 : Shaded Cube. Sphere. Cone. and Cylinder ............................. 80 Figure 7.1.2 : Shaded Torus ......................... 80 Figure 7.1.3 : Shaded Automobile Hood ............... 81 Figure 7.1.4 : Curvature Shading of Automobile Hood . 81 Figure 7.1.5 : Intermediate Stamping of a Styled 'heel ................................ 82 Figure 7.1.6 : Shaded Bumper ........................ 82 Figure 7.1.7 : Shaded Turbine Blade ................. 83 Figure 7.1.8 : Absolute Maximum Curvature Shading of a Turbine Blade ...................... 83 Figure 8.1 : Tree Structure of Non-Zero Blending functions ........................... 98 Figure 3.2 : Tree Structure of Non-Zero Blending functions 0............OOOOOOOOOOOOOO99 Table Table Table Table Table Table Table Table Table Table Table 4.2.1 4.2.2 4.2.3 LIST OF TABLES Analysis of the Operations for the Cox-deBoor Algorithm ................. : Analysis of the Operations for the CURBS Algorithm ...................... : Operations Analysis for Three Coordinates .......................... : Example of Calculated Values of the 2 Component of the Normal .............. Original Control Points ................ Homogeneous Control Points ............. Original Control Points ................ Control Points for Figure 5.1.2 ........ Control Points for Figure 5.1.3 ........ Control Points for Figure 5.1.4 ........ Control Points for Figure 5.1.5 ........ xi 46 46 47 76 116 117 117 118 119 CHAPTER I INTRODUCTION 1.0 INTRODUCTION The goal of this dissertation is to present several algorithms deve10ped for use with CAD/CAM systems. The new algorithms allow for more efficient evaluation of the entire range of rational B-spline curves and surfaces. These algorithms allow a considerable increase in the speed of response of interactive CAD/CAM programs which utilize them. The use of mathematical models to represent sculptured surfaces has become not only a design tool but also a manufacturing tool. Before a model can be used for production or analysis it must be checked for errors. Common errors which occur in these types of models are misplaced corner points. missing patches. slepe discontinuities between patches. and gaps between patches. 1.1 ERROR CHECKING Checking for errors in models can be a very costly process. Two methods in common use are two-dimensional line drawings and proofing runs of trial parts on a numerically controlled machine. In the case of two-dimensional line drawings small errors such as lepe discontinuities are very difficult to detect due to the complexity of the drawings. To check the validity of tool paths for numerical control machining. parts are often milled out of a substitute material to check the model. This is costly because it ties up a mill and an operator. and it increases material costs. Another method for checking the model is to use the computer to generate an accurate shaded image of the model. This method is becoming popular because it saves much of the expense of checking the model [1]. Currently many color shading packages are available for a wide range of surface representations [2]. However. most of these packages sacrifice some of the surface information to produce smoothly shaded pictures. This is done by tiling the surface and utilizing some type of filtering to hide the effects of tiling. These altered surfaces do not represent the true surface definitions and hide many of the errors. Researchers at the A. H. Case Center for Computer-Aided Design have develOped a surface assessment package called MSU COLORSCOPE to produce accurate shaded images. The package allows shading based on various light models as well as surface curvature pr0perties. The MSU COLORSCOPE approach is to calculate and shade the surface at the pixel level using specialized scan line techniques [3] [4] [5]. This creates a very accurate picture of the surface. These pictures can then be used for rapid detection of errors or flaws in the model. 1.2 SURFACE MODELS Rapid growth in the computer-aided design / computer-aided manufacturing (CAD/CAM) field has resulted in the lack of a standard method for defining a three-dimensional surface. In recent years the most popular way has been the bi-cubic patch. But sophisticated applications demand more flexibility than the bi-cubic patch can offer. In many situations the rational B-spline surface is being used [6]. The rational B-spline surface allows a great deal of flexibility and the capability of representing many popular types of surface descriptions in an exactly equivalent form. lhile rational B-spline surfaces cannot precisely represent all mathematical surfaces. they are a very powerful form and are an emerging g; fggtg standard in today's CAD/CAM systems. MSU Colorscope is intended to be a general-purpose surface assessment package. However. prior to this work. it used bi-cubic patch definitions based on the Coons [7] blending functions. This has severely limited applications of the package in today's CAD/CAM community. This dissertation develops both a method for the efficient incorporation of the rational B-spline into a general purpose surface assessment package. and the mathematical forms and algorithms to allow conversion among a variety of surface types. 1.2.1 The Rational B-spline Surface B-splines have only recently been widely used in engineering applications. Perhaps this has been due to the complexity of the basis functions and the computing burden associated with calculating them each time a surface point is to be evaluated. It is well known that the time to calculate the basis functions can be reduced if certain restrictions are placed on the class of B-spline [8]. Is shall show that it is possible to achieve similar time reductions without any constraints placed on the class of B-splines used. An important component of this work is to shade rational B-spline surfaces using scan line methods in an interactive mode. The definition of interactive varies from user to user and task to task; however. for the present purpose. it is taken to mean that one should speed up the shading process as much as possible without sacrificing any of the surface information. The techniques presented here will take advantage of a modern thirty-two-bit computing processor using virtual memory. 1.3 OVERVIEW OF THE DISSERTATION Chapter Two presents a review of the rational B-spline surface. Chapter Three defines a class of B-spline called gghggggd 33112;! and develops an algorithm. called the Converted to Uniform Rational B-Splines (CURBS) algorithm. for calculating with enhanced uniform rational B-splines. A comparison of the commonly used Cox-deBoor algorithm to the CURBS algorithm using the enhanced uniform B-spline is made in Chapter Four. Chapter Five deve10ps an algorithm to transform any rational B-spline into an exactly equivalent enhanced uniform B-spline. Discussed in Chapter Six are some of the numerical problems encountered in the implementation of the B-spline algorithms. Conclusions are drawn and examples given in Chapter Seven. CHAPTER II B-SPLINE CURVES AND SURFACES 2.0 B-SPLINE REVIEW This chapter presents a review of the rational B-spline surface patch. The review starts with a the B-spline space curve and extends to the rational B-spline space curve. Next the B-spline surface patch is defined and extended to the rational B-spline surface patch. Finally. various classifications of B-splines are discussed. 2.1 B-SPLINE SPACE CURVE The B-spline approximation was first presented by I. J. Schoenberg [9] in 1946 for statistical data smoothing. More recent work has been done by deBoor [10]. Cox [11]. Riesenfeld [12]. and Versprille [8]. Cox and deBoor each developed a numerically stable algorithm to evaluate the B-spline points and derivatives. Riesenfeld used the Cox-deBoor algorithm to apply the B-spline to geometric problems in computer-aided design. Versprille used the rational B-spline, but restricted the knot vector so as to produce only a small subset of the B-spline family. 2.1.1 Basis Functions Let P be the Cartesian position along a curve as a function of the parametric variable t. A curve generated using the B-spline basis is given by P(t) - 3 C?1 - Ni.k(t) 2.1.1 i-O where CP are the n+1 control points k is the order of the blending functions n is the number of control points minus one N are the n+1 blending functions The order of the B-spline curve is the degree plus one. If the number of control points exceeds the order of the curve the B-spline will have more than one segment. The segments which make up a B-spline will be called ggbgggggntg. The above blending functions are defined by the recursion formulas [13] "1.1(t) 8 1 if x1 1 t ( xi+1 2.1.2 8 0 otherwise Ni.k(t) 8 (t’xi) Ni,k'1(t) I (11+k-1 ’ 31) 2 1 3 * “1+k't’ N1+1.k-1(t’ ’ (‘i+k ’ *i+1) where xi are the entries of the knot vector 10 The knot vector specifies how the parametric value is distributed along the B-spline. Subsegments start and end when the parametric value crosses the knot values. For example if the knot vector is taken to be X1-[ 8 9 9 9 9 ] i 0 0 0 5 6 7 1 2 3 8 9 10 11 12 13 14 15 0 1 2 3 0 4 5 6 The values of the blending functions for it6 with k=1. k-2. k-3. and k=4 are shown in Figure 2.1.1 as a function of the parametric variable t. a I e n A (a, ? N6.l " T 9 F 3 f N52 i W I O n 1/////~\\\\\_ N v f 5.3 a T | /\ U e a N s I l l l l l l l L *1 6’4 Poromctric Variable Figure 2.1.1 : Blending Functions of Order Up to Four From Figure 2.1.1 we can see that the the blending function is always zero for a value of t from zero to three. As k increases. the blending functions are non-zero for a larger and larger value of t. But the maximum value of the blending functions is decreasing as k is increased. In Figure 2.1.2 the values for all of the blending functions for order k=4 have been plotted. The value plotted for each Ni.4 would be used with with the corresponding CPi to determine values along the B-spline curve. At any value of t there are at most four non-zero blending functions. As t crosses a knot. the blending function Nj'g goes to zero and the blending function Nj+4.4 becomes non-zero. It is at this point that one subsegment ends and another starts. ZZ g—o g—. 45 Scale 1&4 34 E4 Z4 @4 34 M4 14 34 L4 m4 GOG-O< 30-w4-ODC'TI c03-O.30—W ZIZZZZIZZZZIZZZIZZIZZZIZZ Parametric Variable Figure 2.1.2 : Fourth Order Blending Functions Figure 2.1.3 shows the non-zero blending functions for i-6 for orders k=1. k-2. k-3. and k=4 as a function of the parametric variable t. Here t varies from three to four. For order k=4. we see that N3.4 is non-zero at :23 and goes to zero It t=4. while N5'4 starts at zero and ends up non-zero e 13 “DC-O< ao-W-rODCTI $03—°O.30—m Parametric Variable Figure 2.1.3 : Non-Zero Blending Functions for i-6 The effect on the blending functions of repeating interior knots is shown in Figures 2.1.4 and 2.1.5. 14 Scale aac—o< :o—nnscm on-an—m l IZIZZZZIZLZIZZZZZZZZZZZ ? a. b- a 4- ! l l l l L l 3 4 S 6 7 Parametric Variable E’F g—n n1 (D Figure 2.1.4 : Fourth Order Blending Functions Knot Seven Equal to Knot Six The knot vector used for Figure 2.1.4 is 12 - l o o o s 6 7 s s s s 1 1 = o 1 2 9 1o 11 12 13 14 15 (DO 1 2 3 3 4 4 5 6 7 8 The knot vector used for Figure 2.1.5 is :3 - [ o o o 5 6 7 7 7 7 l 1 = o 1 2 1o 11 12 13 14 15 03° 1 2 3 3 3 4 4 5 6 7 8 9 Scale 7A- 2v4 L4 ; l P l l l l J 004' 0 t 2 3 4 5 5 7 Parametric Variable woc—o< DO-nOJCTI oa-aaa—m ZZIZUZEIZIZIZZZZIZZZZZZ ?' 4. Figure 2.1.5 : Fourth Order Blending Functions Knots Six. Seven. and Eight Equal The non-zero blending functions of Figures 2.1.2. 2.1.4. and 2.1.5 for t between two and four have been plotted in Figure 2.1.6. 16 B 1 I T e 3 i .................... , ---------------------------------- m ...",l ......... ”‘\ ‘ ' r".’. ........... .‘.“\ O n K .;.;."""T‘ .‘.> g a “"--o--- .-.."Oe...-...-- ___‘____.-——-” 1 I 11 F u ........ g ........ n .......... ’.’.’.’ \ )II” ‘‘‘‘‘‘ ‘ I’”. c ’J’ ...... , ‘\\ ‘,2~.~ t .’. ....... ”’a ‘. // ‘\“. b m ‘m. ‘ / o d ‘ .....4‘; I, _-_E#. n [m I” ‘\\ v ’I’ \\\ O .......... .- ------- \. I \ /’—~\\ 2”, | .,.,.< ..... ”3g. /\(\ ’/,~<-\\ u P" ...... 3.”’ \‘ \\ o/ \ C e a '-" ........ \‘V/ 4’" “~~~ i‘i.‘ “I. - """" J’ __ Q a L L l l I 1L8 2J5 2L8 :15 440 Parametric Variable Figure 2.1.6 : Fourth Order Blending Functions a) Knot Vector X1 b) Knot Vector X2 c) Knot Vector X3 Iith no repeated knots there are five non-zero blending functions in this range. Iith one repeated interior knot there are six non-zero blending functions. At t-3 two of the blending functions go to zero as Opposed to one with no repeated interior knots. Iith the interior knot repeated two times. three blending functions go to zero at t-3. 17 Each control point is associated with an Ni and each "1 is associated with a knot x1 in the knot vector. However. it is difficult to associate a control point with a knot in the knot vector. As shown later. only a number of the control points equal to the order of the curve affects a given subsegment. But the number of knots which affect the blending functions for that subsegment is greater than the number of control points. A way to think of this is that a given control point affects more than one subsegment. Because of continuity constraints between subsegments. the control point must be able to look ahead and look back to influence the other subsegments. The control point uses the knot vector to do the looking. Repeated knots in the knot vector can then be viewed as walls to limit the number of subsegments the control point can see forward or backward. The only restriction placed on the elements of the knot vector is that xi 2 xi-1. For a given set of control points. an infinite number of curves could be generated. depending on the values assigned to each xi. It is this fact which causes the need to evaluate equations 2.1.2 and 2.1.3 or some version thereof. each time a point is to be found. However. not all of the Ni,k are non-zero for a given subsegment. In fact only a number equal to the order (k) of the curve have non-zero values. These Ni k are the 18 only ones that need be evaluated for that subsegment. The k non-zero Ni.k are the same for all t such that 1j 5 t < xJ+1s points on this part of the curve are said to form the jth subsegment of the curve. It is possible to derive a formula for the Ni,k for a given subsegment without an explicit value of t. provided the knot vector is given. These formulas are valid for any t in the range xj S t ( xj+1. However. there are an infinite number of possible knot vectors. These formulas would have to be derived each time a subsegment was to be evaluated for a new knot vector. It has been shown in [10] and [11] that the use of equations 2.1.2 and 2.1.3 directly can cause numerical instabilities when the entries of the knot vector become irregularly spaced or the order of the curve becomes large. Thus Cox [11] and deBoor [10] independently develOped the same algorithm which is numerically stable. But this algorithm returns values of the Ni.k and/or their derivatives only for particular values of t. and must be recalculated entirely for each new t. No shall return to the stability of equations 2.1.2 and 2.1.3 later in the discussion of the restrictions that can be placed on the knot vector. 19 2.1.2 Rational B-spline Space Curves The extension to rational B-spline curves is handled by assigning a positive non-zero weight to each of the control points. The weights are used to form a homogeneous coordinate system. Each control point CPi(x.y.z.1) is assigned a weight 'i° For calculations with this control point. it 1‘ multiplied by ‘1 to form a new control point CP;(wx,'y,'z,'), The cp' points are used in equation 2.1.1 to calculate a homogeneous point P'(wx.wy.wz.w). The three-space point P(x.y.z,1) is calculated by dividing each coordinate of P. by the calculated value of w. X(t) = IX(t)/l(t) Y(t) ‘ IY(t)/I(t) 2.1.4 Z(t) - urn/7(1) It is hard to think of a curve in four-space. but one way to visualize this is to use a curve in two dimensions. The following three figures show an example of a two-dimensional rational B-spline curve. 2O 30 wa~o=-c1-ooo mcooaoooaoz -10 1 1 1 1 0 2 4 6 8 10 Parametric Variable T Figure 2.1.7 : Homogeneous Coordinates vs Parametric Variable The data for Figures 2.1.7, 2.1.8 and 2.1.9 are given in Appendix A. Figure 2.1.7 shows the values of IX. 'Y and I as a function of the parametric value t. The three ordinates of Figure 2.1.7 are plotted in three-space in Figure 2.1.8. Figure 2.1.9 is a plot of the resulting curve in K-! space; also plotted is the curve which would result if all the weights were equal. 21 WY WX A Figure 2.1.8 : Three-Space Plot of Homogeneous Coordinates A B-spline is then also a rational B-spline with all the veishts set equal to one. A 21512111 11112211 §:anlins will have one or more weights not equal to one or or; a cri. Assigning unequal weights to the control points introduces an additional degree of freedom. This degree of freedom can be used to produce curves or shapes that are 22 not attainable with non-rational curves. circular arcs. such as conics or 3 2 ~ , Y ................ 8-spline C Rational B-spline O O f - <1 1 l : n s 0 : t g e I _ '. - 1 1 - .4 1 1 -t.5 -l.0 -8.5 8.0 8.5 1.8 X Coordinate Figure 2.1.9 Plot of X-Y Curve 2.2 B-SPLINE SURFACE The extension from the space curve to the surface patch is done by adding a second parametric coordinate and associated knot vector. Points on the by surface are given 23 'here Ni.k(t) ‘ lj,h(') 2.2.1 CP are the (n+1)-(m+1) control points k is the order of the curve in the t direction h is the n is the in the t direction m is the in the s direction N are the n+1 blending M are the m+1 blending The control points are new column. The number of rows number of columns. nor does the have to equal the order in function is associated with function is associated with B-spline curve. non-zero. AnalOgous to resulting submatrix of the control point matrix 13121122. t 22122; 2112!- each only a given number of these functions the subsegment of the curve. order of the curve in the s direction number of control points minus one number of control points minus one functions of order k functions of order h associated with a row and does not have to equal the order in the s direction the t direction. A basis row and a different each column. As with the are the defines a The entire B-spline patch will be referred to as 24 2.2.1 Rational B-spline Surfaces More complex surfaces may be constructed if a weight is also assigned to each control point of the B-spline surface. To determine points on the surface. equation 2.2.1 is used in the homogeneous coordinate system as described above. The resulting surface point in three-space is then found by equation 2.1.4. 2.3 KNOT VECTOR RESTRICTIONS The only restriction placed on the knot vector is that each entry must be a real number equal to or greater than the one immediately preceding it. This type of knot vector is referred to as a gggggifggg gag; yggtgg and produces B-splines called Nonuniform Rational B-splines (NURBS). Knot vectors are often gxtggggg by adding knots to the beginning and end equal to the first and last entries of the vector respectively. Repeating the end knots causes the curve to pass through the end control points. In this case the end control points are called a gggggtgig gags. 25 2.3.1 Uniform Knot Vectors Placing a restriction on the knot vector that the entries are only allowed to take on successive integer V91303 (31+1 = x1 + 1). produces a knot vector known as a 23132;; $32; 122325. This knot vector produces what is called a ggifggg bggi; [12]. which leads to Uniform Rational B-splines (URBS). This knot vector can also be extended to form :n 22222222 2222222 2222 222222- The use of uniform or extended uniform knot vectors has some side effects. In particular. the parameterization of the curve is 32; uniform with respect to arc length along the curve. and local refinement is difficult. However. a uniform knot vector makes evaluation of the basis functions easier because equations 2.1.2 and 2.1.3 can be used directly. Thus an equation for each subsegment can be derived without a specific value of the parametric variable. The nonuniform parameterization which often accompanies a uniform knot vector is undesirable in some applications such as numerically controlled machining. in which it is desirable to have a constant step size of the parametric variable produce a constant step size along the curve. However. nonuniform parameterization does not 26 affect the evaluation of surface properties for the purposes of shading. 2.3.2 Enhanced Uniform Knot Vectors To restrict the knot vector to be uniform or extended uniform for the purposes of a CAD/CAM modeler would over-restrict the class of B-splines that could be used. For example. this would not allow repeated interior knots. which can be used to introduce slope discontinuities between the subsegments. Ihile slope discontinuities can instead be introduced by repeating control points in the surface definition. this type of discontinuity has the undesirable preperty of a zero length normal. To allow for s10pe discontinuity and also retain the necessary surface properties. knots must be allowed to be repeated. The effect of repeating knots is to relax the condition on continuity of derivatives at the ends of subsegments -- one less derivative for each time a knot is repeated. However. the surface properties are defined. Repeating knots also makes some subsegments of the curve have zero length. A knot vector with repeated interior knots produces what is called a gpbgplipg bggig [12]. 27 In order to include a more general class of B-spline. the uniform condition ‘i+l - xi + 1, is relaxed to also 4110' 81+1 - 1i- If 3i+1 is equal to xi then the knot is said to be repeated. It can be repeated up to k-l times. where k is the order of the curve. This type of knot vector '111 be called an 22222222 2222222 2222 122222 vhich will produce 22222222 2222222 22222222 222222222 (E0338)- After a review of rational B-spline curves and surfaces. restrictions that can be applied to the knot vector were discussed. The uniform condition of the extended uniform knot vector was relaxed to introduce the enhanced uniform knot vector. The next chapter will discuss the enhanced uniform knot vector and present the CURBS algorithm for calculating with it. CHAPTER III CURBS ALGORITHM 3.0 ENHANCED UNIFORM RATIONAL B-SPLINE The last chapter defined some restrictions on the knot vector. A class of Rational B-spline. the enhanced uniform. was introduced. This chapter describes the advantages of using EURBS and an algorithm for calculating with them. Finally. it presents the limitations of using this class of B-spline. as Opposed to the more general NURBS. 28 29 3.1 ADVANTAGES OF UNIFORM RNOT VECTOR In most engineering applications. surfaces free of ripples are desired. The higher the order of the curve used. the more difficult it is to eliminate ripples. To produce a rippleless surface. the order of the curves used is usually restricted to six or less. Therefore. if a uniform or enhanced uniform knot spacing is used. then equations 2.1.2 and 2.1.3 can be used to calculate the basis functions with acceptable accuracy. This allows a simple equation of degree k-l. for each of the basis functions. to be formulated before a specific surface is to be evaluated. An example of this formulation is given in Appendix B. There are. however. an infinite number of possible knot vectors. resulting in an infinite number of blending functions. Next it will be shown how to translate the parametric value to reduce the number of possible blending functions to a finite tractable set. for any given order k of B-spline. 3.1.1 Knot Vector Translations Making the EURBS restriction to the knot vector allows for the derivation of all possible blending functions of a given order. In order to do this derivation. the 30 parametric variable is restricted to values between zero and one for each subsegment. This is done by translating the knot vector for each non-zero length subsegment so that th. 1th knot is zero and the i+1 knot is one. Repeated knots form zero-length subsegments which are of no interest. Therefore the transformation is possible for any non-zero subsegment. For example. the knot vector ”H II 60 HO NO 950 for the third subsegment. Appendix C shows that this is a linear mapping of the parametric variable and has no effect on the values of the resulting blending functions for any given value of the parametric variable. 31 3.1.2 Enhanced Uniform Knot Vector Compression For a given subsegment of a curve. only a finite portion of the knot vector is actually used in the calculation of the non-zero basis functions for that subsegment. If the recursion of equation 2.1.3 is followed. the portion of the knot vector which affects the blending functions can be determined. If k is the order of the curve. the effective kggt vector for a subsegment reaches k-2 knots back from the starting knot and k-l knots forward of the starting knot. for a total of 2-(k-1) knots. In the above example. if k=4 the effective knot vector would have a length of 2-(4-1) or six. For the first subsegment. it would be X. = [o o o 1 2 2]. For the second subsegment. it would be X0 = [O O 1 2 2 3] which could be shifted to x; 2 {-1 -1 o 1 1 2]. For the third it would be 32 1,-[122344] which could be shifted to x;=[-1oo122]. Using parameter translation and the finite length of the the effective knot vector. the set of possible effective knot vectors has been reduced to a finite set. given a bound on the order of the spline. The number of possible knot vectors in the set is determined by the order k of the curve and is 2(2'k-4). For k-4. there are 2(2‘4-4) or sixteen possible knot vectors for a subsegment. For k=5. there are sixty-four possible knot vectors. Each of the 2(2°k-4) knot vectors yields a unique set of blending functions. In order to determine easily which set of blending functions should be used for a given subsegment. each subsegment is given a label or pointer to the correct set of blending functions. The label is determined by compressing the effective knot vector into a binary-encoded label. 33 3.1.3 Knot Vector Compression.Algorithm The effective knot vector is compressed into a unique binary-encoded label by the gggpggsgigg glggrithg given in Appendix D. The algorithm compares successive entries of the effective knot vector to define a sequence of Zak-4 bits. which are reversed and stored as an integer label. This label will be referred to as a gpbggggggt label. Since a unique set of effective knot vectors for a given order has been identified. a set of blending functions for each unique knot vector can be formulated. as in Appendix Now that all the enhanced uniform blending functions have been identified and labeled we are ready to present the CURBS algorithm for calculating with rational B-splines. 3.2 THE EVALUATION ALGORITHM This section presents an algorithm we call the Converted to Uniform Rational B-Splines (CURBS) algorithm for evaluation of enhanced uniform rational B-splines. It assumes that all the blending functions have been preformulated and stored in matrix fashion according to the 34 labels determined from the effective knot vector. Note that the formulation of the blending functions is done only once and the results are stored for later use. The algorithm as stated below also assumes that the subsegment labels are calculated or read in with the control points that define the surface. The algorithm requires the following information to be given : 1) Control points for the subsegment : CP' vector 2) Subsegment label : LAB 3) Order of the curve : K 4) A strictly rational indicator : RATNOT 5) Parametric value (between 0.0 and 1.0) : T The following informal explanation preceding the algorithm specification may assist the reader in understanding the algorithm. Step 1: Use LAB and K to retrieve the proper set of blending functions : [MATK] Step 2: Construct the t vector : {t vector ] ' [ t"1.t"2.....t.1 ] Step 3: Evaluate blending functions 35 [Bf] ' [t vector] 0 [MATK] Step 4: Calculate homogeneous point up} - {3:} - {cp'} Step 5: Calculate three-dimensional point if strictly rational [P] ' [VP] / I For K-4 the CURBS algorithm can be coded as follows: (the coordinates of the control point CP' for a subsegment are contained in the PK. PT. P2. and PI vectors). T2 - T‘T T3 ' T‘T2 BF(1) ' T3’MAT4(1.1.LAB) + T2‘MAT4(2.1.LAB) + T‘MAT4(3.1.LAB) + MAT4(4.1.LAB) T3‘MAT4(1.2.LAB) + T29MAT4(2.2.LAB) BF(2) + T‘MAT4(3.2.LAB) + MAT4(4.2.LAB) BF(3) = T3‘MAT4(1.3.LAB) + T2'MAT4(2.3.LAB) + T‘MAT4(3.3.LAB) + MAT4(4.3.LAB) BF(4) ' T3‘MAT4(1.4.LAB) IX ' BF(1)‘PX(1) + BF(2)‘PX(2) + BF(3)‘PX(3) + BF(4)‘PX(4) MY = BF(1)‘PY(1) + BF(2)‘PY(2) + BF(3)‘PY(3) + BF(4)‘PY(4) 36 '2 I BF(1)‘PZ(1) + BF(2)‘PZ(2) + BF(3)‘PZ(3) + BF(4)‘PZ(4) IF(RATNOT) THEN I I BF(1)‘P'(1) + BF(2)‘PW(2) + BF(3)‘PU(3) + BF(4)‘PU(4) X I 'X/' Y I III! 2 I 'Z/' X I 'X Y I II Z I '2 END IF If derivatives of the point with respect to the parametric variable are required. the algorithm is modified. The appropriate derivative of the t vector is taken and used in step 3 to evaluate an alternate set of blending functions. This set of blending function derivatives is used in step 4. If the surface is strictly rational. then the chain rule must be employed to replace step 5 to calculate the derivative. For example. to calculate dX/dt it would be necessary to use the following equation: 37 dX/dt - dIK/dt ? s - 7x . dI/dt e ('2) These modifications are shown in the following piece of code. T2 I TIT THREET2 I 3‘T2 T'OT I 2‘T DBF(1) I THREET2‘MAT4(1.1.LAB) + T‘OT‘MAT4(2.1.LAB) + IAT4(3.1.LAB) DBF(2) THREET2‘MAT4(1.2.LAB) + T'OT‘MAT4(2.2.LAB) ... MAT4(3.2.LAB) DBF(3) I THREET2‘MAT4(1.3.LAB) + T'OT‘MAT4(2.3.LAB) + MAT4(3.3.LAB) DBF(4) I THREET2‘MAT4t1.4.LAB) IF(RATNOT) THEN T3 I T‘T2 BF(I) I T3‘MAT4(1.1.LAB) + T2‘MAT4(2.1.LAB) + T‘MAT4(3.1.LAB) + MAT4(4.1.LAB) BF(2) TSIMAT4(1.2.LAB) + T2IMAT4(2.2.LAB) + T'MAT4(3.2.LAB) + MAT4(4.2.LAB) BF(3) T3‘MAT4(1.3.LAB) + T20MAT4(2.3.LAB) + T‘MAT4(3.3.LAB) + MAT4(4.3.LAB) BF(4) I T3IMAT4(1.4.LAB) IX I BF(I)‘PX(1) + BF(2)‘PX(2) + BF(3)‘PX(3) + BF(4)‘PX(4) 38 NY I BF(I)‘PI(1) BF(2)‘PI(2) + BF(3)‘PY(3) BF(4)‘PY(4) 'Z I BF(1)9PZ(1) BF(2)‘PZ(2) + BF(3)‘PZ(3) BF(4)‘PZ(4) I I BF(l)‘P'(1) BF(2)‘PI(2) + BF(3)‘P'(3) BF(4)‘PI(4) END IF D'XDT I DBF(1)‘PX(1) + DBF(2)‘PX(2) + DBF(3)‘PX(3) + DBF(4)‘PX(4) DIYDT I DBF(1)‘PI(1) + DBF(2)‘PY(2) + DBF(3)‘PI(3) + DBF(4)‘PY(4) + D'ZDT I DBF(1)‘PZ(1) DBF(2)‘PZ(2) ... + DBF(3)‘PZ(3) DBF(4)‘PZ(4) IF(RATNOT) THEN D'DT I DBF(1)‘PI(1) + DBF(2)‘PI(2) + DBF(3)‘PI(3) + DBF(4)‘P'(4) '2 I "I DXDT I D'XDT/I - 'X‘D'DT/(IZ) DYDT I D'YDT/U - 'I'DIDT/(VZ) DZDT I D'ZDTI' - 'Z‘D'DT/(U2) ELSE DXDT I DMXDT DYDT I D'YDT DZDT I D'ZDT END IF 39 The above algorithms are given for a curve. To evaluate surface points the CP' vector is replaced by a CP' matrix which is post-multiplied by the second set of blending functions corresponding to the other parametric variable. Unlike the Cox-deBoor algorithm. this algorithm does not rederive the blending functions as they are being evaluated for a surface point and/or derivative. However. the application of algorithm is restricted to rational B-splines defined with enhanced uniform knot vectors. 3.3 LIMITATIONS OF ENHANCED UNIFORM KNOT VECTOR Having the formulas for all of the blending functions predefined in the program appears to be more efficient for evaluating B-splines and derivatives than building and evaluating the blending functions for each point. However. the use of EURBS appears to limit the B-splines which can be described to a subset of the B-splines with a more general nonuniform knot vector. But Chapter Five presents an algorithm to insert knots into a knot vector so that any knot vector can be transformed to the enhanced uniform format. removing the apparent limitation. 40 Ie have discussed the advantages of a uniform knot vector. and reduced the infinite set of enhanced uniform knot vectors to a finite set by knot vector translation. The effective knot vector was defined and then compressed into a subsegment label. An evaluation algorithm was then presented. In the next chapter we shall discuss the performance differences between the above CURBS algorithm and the Cox-deBoor algorithm. CHAPTER IV CURBS ALGORITHM ANALYSIS 4.0 ALGORITHM COMPARISON The last chapter presented the CURBS algorithm for calculating with EURBS. This algorithm relies on the use of preformulated blending function matrices and the ability to determine which matrix should be used for a given subsegment of the curve. In this chapter the CURBS algorithm will be compared to the well-known Cox-deBoor algorithm for calculating with B-splines. The comparison is made by classifying each algorithm according to order of 41 42 complexity and then according to the number of machine operations required. 4.1 COX-DEBOOR ALGORITHM Cox [11] and deBoor [10] separately developed the same numerically stable algorithm for computing the basis functions for B-splines. This algorithm is the standard of today's B-spline algorithms and is in wide use in many CAD/CAM systems. For the purposes of analysis the algorithm as presented by deBoor [10] is reproduced below: Set N(l.1) I l for s I 1 ..... k-l. do: set DP(s) - ti+s-t. Dm(s) - t-ti+1-s. set N(1.s+1) I O: for r I 1 ..... s. do: . set M I N(r.s) I (DP(r)+DM(s+l-r)). . : set N(r.s+1) I N(r.s+1) + DP(r)‘M. .............set N(r+1.s+1) I DM(s+l-r)IM. where k is the order of the blending functions. 43 The columns of N are used to evaluate the B-spline and/or its derivatives as shown below. The B-spline is then evaluated by 2-1 P(t) I 2 cpi - Ni'k(t) 4.1.1 i-o Ihere Ni'k 1. the 2th column of N in the above algorithm To evaluate derivatives. the appropriate column of N would be used. However. N is a triangular matrix. so a set of derivative control points must first be calculated to be used with a specific column of N. The general formula given by deBoor for this evaluation is 44 2-1-1 F(3)(t) a (2-1) ... (k-j) E A15) . Ni.k_j(t) 4.1.2 i=0 Ibere j is the desired derivative 0 S j < k Ni,k-j is the (k-j)th column of N (0) . A1 cpi All) - (Ali-1’ - 11231)) / (ti+k_j - t.) 4.2 BIG-OH ANALYSIS One way to classify an algorithm is the big-oh [14] method. The big-oh analysis determines the order of complexity of an algorithm. This analysis performed on the above Cox-deBoor algorithm results in a classification of 12 or 0122). The same analysis of the canes algorithm given in the previous chapter produces O(k). Big-oh analysis indicates that the computation time using the CURBS algorithm increases linearly with k. It also indicates that computation time using the Cox-deBoor algorithm increases as k2. This savings of an order of magnitude becomes significant when k becomes large. For 45 Our application. the size Of k is usually six or less. so big-Oh analysis is not a very complete measure Of the relative performance Of the algorithms. It does provides an indication of the relative behavior for cases in which k is required to be large (>>6). Ihile the size of k is usually small. the number Of points and/or derivatives to be calculated per picture is usually very large. Therefore the potential for considerable time savings is good if a more detailed operation analysis indicates a significant difference in relative efficiency for the usual range of values for k. 4.2.1 Operation Analysis 'hen the big-Oh analysis is not sufficient to capture the differences between two algorithms. a more detailed analysis is needed to determine performance. A second way to classify algorithms is to determine the number Of machine Operations required. The following tables give the number Of Operations required by each algorithm for kI3 and kI5. where an Operation is defined as an addition. subtraction. multiplication. division. or assignment of value. The tables present the number Of operations required to evaluate a point (nIO). or the 1th derivative. or evaluate the point and n derivatives. 46 Operations for Operations for the jth derivative point and n derivatives l 1-3 I 1-2 I 1-1 I n-o l n-l l n-2 l n-s l l l l l l_ l l _ l l l l l l l l l 2-3 I 1 l 44 l 76 l 76 l 117 l 160 l 161 l l l l l l l l l l l I l l l l l k-s : 146 : 131 l 222 l 212 l 295 l 409 l 529 I Table 4.2.1 : Analysis of the Operations for the Cox-deBoor Algorithm Operations for Operations for the jth derivative point and n derivative S I jI3 I jI2 I jIl I nIO I nI1 I nI2 I nI3 I l _l l l l l l l l l l l l l l l 1-3 l 9 l 9 l 16 l 20 I 36 l 45 l 54 l l l l- _l l 2| l l l l l l l l l | k-s : 30 l 42 : 56 l 54 l 106 i 146 l l Table 4.2.2 : Analysis of the Operations for the CURBS Algorithm This analysis was performed for the evaluation of a single coordinate of a B-spline curve and/or the derivatives of a single coordinate. The implied integer multiplication and addition for matrix subscript indexing are not included. To evaluate other coordinates only 47 requires that the evaluated blending functions be used on the remaining coordinates of the control points. Thus the other coordinates would be quicker to evaluate once the first one is done. Table 4.2.3 contains the Operations required by each algorithm to evaluate three coordinates and three coordinates with derivatives. i=3 l i=5 l l I nIO I nI1 I nI2 I nIO I nIl I nI2 I l l |_ l l l l l l l l l l l c-n l 104 l 227 l 356 l 280 l 529 I 871 l ______ l l_ l l l l =l l l I l l I l cunns : 32 l 60 I 31 l 74 : 146 i 206 : I I I l I I l l l I Table 4.2.3 : Operations Analysis for Three Coordinates This table indicates that the Cox-deBoor algorithm requires more effort to evaluate the additional derivatives. This is because the Cox-deBoor algorithm must calculate the Aij) points for each coordinate for each derivative to be used with the columns of the N matrix. The CURBS algorithm always uses the same set Of control points for all evaluations. 48 For a surface. the algorithms would have to be used a second time on the second parametric variable and the resulting evaluated blending functions used on a matrix of control points. The above tables show that the CURBS algorithm uses fewer machine Operations than Cox-deBoor for B-spline evaluations. Vhile the actual time difference for using each algorithm will vary from processor to processor. the reduction of Operations required does imply an approximate reduction in computation time. The CURBS algorithm does require additional storage for the precalculated matrices. For kI4 the additional storage required is sixteen 4-by-4 matrices. For kI6 the storage required is 256 6-by-6 matrices. The additional storage is a minor factor with the low cost of computer memory and the availability of virtual memory in computers. The CURBS algorithm has been shown to produce a three-to-one reduction in machine Operations over using the Cox-deBoor algorithm. The CURBS algorithm has no numerical stability problems. It does. however. require that the blending function matrices be constructed in a numerically stable manner. The numerical problems of using equations 2.1.2 and 2.1.3 for construction of the the blending 49 functions presented in [10] and [11] are avoided because of the uniform nature of the knot vector being used and because the order of the curves involved is not large. The Cox-deBoor algorithm has been reviewed and a big-Oh analysis performed comparing it to the CURBS algorithm. The big-Oh analysis indicated a major complexity difference between the two algorithms. but was inconclusive as to the importance of the difference in commonly encountered situations. A more detailed Operation analysis was performed which indicated that the CURBS algorithm would give a three-to-one Operation savings over Cox-deBoor. The Only limiting factor to the CURBS algorithm is that the knot vector must be enhanced uniform. which limits its application to a subset Of the NURBS. In the next chapter. we will overcome this limitation by develOping a transformation to transform an arbitrary NURBS surface definition to an equivalent EURBS surface. CHAPTER V CONVERSION ALGORITHM 5.0 NURBS-TO-EURBS CONVERSION In the preceding chapter the CURBS algorithm was shown to have a significant performance advantage over the Cox-deBoor algorithm. However. it was also noted that the surface must be defined with an EURBS knot vector which is a subset of the NURBS knot vector. This appears to be significant in view of the fact that a uniform knot vector does not necessarily produce a curve with uniform arc length parameterization. a desirable attribute for 50 51 applications such as generation of tool paths for numerically controlled machines. In order to produce a curve with uniform arc length parameterization. the knot spacing must reflect the relative arc lengths of the subsegments. which implies a nonuniform knot vector. The above limitations on the use Of the algorithm can be eliminated by converting surface representations based on nonuniform knot vectors into exactly equivalent surfaces based on our enhanced uniform knot vector. This chapter presents the design of an algorithm to do the NURBS-to-EURBS conversion. Discussion of the initial attempts to rescale without expanding the knot vector are followed by the exact reformulation of EURBS based on NURBS. After presentation of the algorithm. the additional storage required to represent the surface in the EURBS form is discussed. 5.1 INITIAL EXPLORATION OF THE CONVERSION PROBLEM Two approaches to converting NURBS to EURBS were tried before a suitable method was derived. These initial approaches were tried to determine whether a conversion was possible which would use the same number Of control points as the original NURBS surface. The next section presents a 52 discussion of one initial approach and includes examples showing the shortcomings. 5.1.1 Knot Vector Rescaling The first approach explored was the use of an enhanced uniform knot vector with the same number Of entries as the nonuniform one to be replaced. This implies that both definitions use the same number of control points (although relocated) in the surface definition. This method results in a rescaling of the parametric variable for each non-zero tnbsos-ont from A I ‘1+1 - xi to A = 1. thus changing to EURBS-based blending functions. Unfortunately. this method is limited in the type Of knot vector to which it can be successfully applied. This is shown in the following proof. Define a nonuniform knot vector x =[x.811eeeelxn] corresponding to the parametric variable t. Let Y be an enhanced uniform knot vector Y = [Y..y1.....yn] corresponding to the parametric variable linear 53 s. The following relationship holds between s and t for the 1th subsegment: 31' (t .. xj)/(xj+1 - x1) ‘I’ (j - I). If continuity of the first derivative is assumed between the 1th and (j+1)th subsegments. then: 'j I (t - xj)/(xj+1 - xi) + (j - 1) 'j+1 ‘ (t ‘ xfill/(814.2 ' 1j+1) + j Equating these relations at the junction Of the two subsegments. t I xj+1 yioldg the relation that must exist for the knots. (3j+1 - xj) a (Ij+2 - 1j+1) Th0 ‘j entries can take on values such thlt (13+1 - xj) I (21+; - xj+1). This means that dsj # dgj+1 gt the junction. which is a contradiction. Thus this method would not be able to yield an equivalent enhanced uniform representation for all nonuniform knot vectors. Another way to view this rescaling is to think of it mapping from subsegments of t-space to a new 54 s-space. Each non-zero distance A between knots in the t-space is mapped to A - l in the s-space. Because the A's in the t-space do not have to be constant. the napping to s-space is. in general nonlinear. Thus to use this form of rescaling. the knots in each effective knot vector must have a constant or zero spacing. which will insure the linear napping. An effective knot vector with a constant or zero knot spacing will be said to htvo the assesses £215; arersrsx- 5.1.1.1 An Exanple 0f Approximate Knot Vector Rescaling The following figures will demonstrate the above rescaling. Figure 5.1.1 contrasts the use of a nonuniform knot vector and a uniform knot vector for the same set of control points. ( Appendix E presents the data for these figures). Figures 5.1.2 and 5.1.3 show the results of solving for new control points to try to fit the NURBS CIIVO. 55 Y C o o r q I n o t e “'2‘ <>- — - - Uniform - Cl Nonuniform H I I I I I I . IUITIVUIUIIIIIIII'UI1UTTITTUI| 0.0 0.5 1.0 1.5 2.0 2.5 3.. X Coordinate Figure 5.1.1 : Uniform vs Nonuniform Curves The control points for Figure 5.1.2 were found by requiring the curve to pass through the ends of each segment and match the slapes on the ends. 56 lo'—— 0 D “08“— Y - C 0.06-- o r - 9 A'o‘—— o i 5 J e 000 “'2‘ O- - - - Uniform ~ D——- Nonuniform 0" VIIIIIIIIITITIIFTTTIIIIFI’UFUI: 0.0 0.5 1.0 1.5 2.0 2.5 3.0 X Coordinate Figure 5.1.2 : Uniform Fit Using End Slepes The control points for Figure 5.1.3 were found by fitting the first segment and then requiring the curve to pass through the end points of the remaining segments. By design there are no discontinuities. but the fitted curves do not represent the original definition. 57 Y C o o r a I n a t e LB. 0- - - - Uniform 0 [3-——-Nonuniiorm 0" TVIT:IIUI:UIIU:Irrrgrll'lliltr{ 0.0 0.5 1.0 1.5 2.0 2.5 3.0 X Coordinate Figure 5.1.3 : Uniform Fit Using First Subsegment In Figure 5.1.4 we have used ordinary least squares [16] to try to fit the nonuniform data to a model using the blending functions based on the enhanced uniform knot vector. The plot is the uniform curve that fits the original nonuniform curve with the least amount of error between input data points and corresponding points on the uniform curve. 58 L0 c1<> n - ‘0 Y C o a r a I n a 1 e D o 0.2- O- — - - Uniform Cl— Nonuniform 0" UIUT I'll [Uifi'll‘ IIIUIIIIUI mm 0J5 1J0 L5 am 2:5 2%! X Coordinate Figure 5.1.4 : Uniform Fit Using Least Squares 5.2 EXACT NURBS-TOIBURBS CONVERSION The following conversion technique was derived after the initial explorations described above showed the necessity for additional knots and control points. 'hile it may produce a surface with more control points than the original surface. it is an exact conversion and the 59 resulting surface can take advantage of all the EURBS properties. 5.2.1 Expanding The Knot Vector The method described in section 5.1.1 would work if each effective knot vector had the constant delta prOperty. In order to obtain the constant delta property. new knots are inserted or repeated to expand the new EURBS knot vector. This implies that additional control points will also be inserted into the surface definition. Only zero-length subsegments are added. so the total number of non-zero subsegments remains the sane. The insertion is done such that the resulting effective knot vector for each subsegment has the constant delta pr0perty. This forces the mapping or rescaling of each subsegment to be linear. The additional enhanced uniform control points give the flexibility needed to match the NURBS curve exactly. The result of this method is shown in Figure 5.2.1. The resulting EURBS curve is identical to the NURBS curve. and of course retains the same level of continuity between subsegments. That is. all continuous derivatives of the 60 nonuniform spline remain continuous in the uniform representation. The algorithm for expanding the knot vector and finding the new set of control points is given below. Y C o o r q I n a 1 e “'2' O- - - - Uniform [3-———-Nonuniform '00 [TlerIlrllllflrrflllrlirllliil 0.0 0.5 1.0 1.5 2.0 2.5 3.0 X Coordinate Figure 5.2.1 : Uniform Fit By [not Vector Expansion 61 5.2.2 Algorithm for Conversion By Knot Vector Expansion ( NURBS-To-EUBBS Algorithm ) To fit a unique curve of a given order k. k pieces of information must be given. The information can take the form of k points along the curve or a point and k-l derivatives or other such combinations. In this algorithm. the point at the end of the subsegment is used along with successive derivatives as necessary. The following outlines the steps involved in the conversion, and the code for this algorithm is given in Appendix F. Step 1 : C0py the NURBS knot vector into the new EURBS knot vector. Step 2 : For each successive effective knot vector in the EURBS knot vector Do: Ihile not constant delta property Do: Duplicate the last knot with constant delta property. End 'hile Label the subsegment per section 3.1.3. End For. 62 Step 3 : For each successive subsegment in the new EURBS knot vector Do: Use label to determine the number of undetermined control points. Evaluate the end point and necessary derivatives of the NURBS curve for this subsegment. Determine the unknown EURBS control points for this subsegment. End For. This algorithm is given for a curve. To use it for a surface definition. the algorithm is first used on each column of the nonuniform control point matrix with the corresponding knot vector. The resulting expanded columns are used to form an intermediate matrix of control points and the algorithm is then used on the rows of this matrix with the other knot vector. The resulting expanded rows contain the uniform control points that define the mathematically equivalent surface. 63 5.2.3 Storage Considerations To use the CURBS algorithm. a submatrix of the control point matrix must be retrieved. The retrieval time will be determined by the size of the submatrix and not the size of the mother matrix. The retrieval of the control points was not included in the operation analysis done because it was assumed to be the same for both algorithms. The additional control points needed for the EURBS definition will have no effect on the performance of the CURBS algorithm. However. there is a possibility of a slight increase in the paging time in a virtual memory environment because of the larger number of control points used is $252. The NURBS-to-EUBBS conversion will represent any surface originally based on a nonuniform knot vector with an equivalent surface based on an enhanced uniform knot vector. The price for doing this is that there will be more control points to store for the enhanced uniform representation. The amount of additional storage depends how many new control points are created. An additional k-l control points will be inserted for each subsegment which does not initially possess the constant delta prOperty. The number of non-zero length subsegments has not increased. and since the knot vector is stored as a set of labels. for non-zero length subsegments in EURBS form. 64 there is no additional storage required for the longer knot vector. The extra virtual memory necessary to store the surface in this form is far outweighed by the time savings resulting from being able to calculate points and derivatives based on prestored blending functions. Notice that the introduction of new zero-length subsegments does not increase the number of calculations required since no evaluations are done in these segments. and the cost of using the new control points is already included in the evaluation of the CURBS algorithm. Thus. the CURBS algorithm offers a fairly direct tradeoff of storage for computational time. The process of transforming a NURBS surface to an equivalent EURBS surface is done only once. after which the EURBS surface may be manipulated quickly and efficiently by the CURBS algorithm. The CURBS and NURBS-to-EURBS algorithms which have been develOped can be used to evaluate any rational B-spline curve or surface. The implementation of these algorithms in a general-purpose surface assessment package requires develOping the routines for the necessary calculations. The next chapter discusses some numerical considerations of the implementation. CHAPTER VI INPLENENTATION CONSIDERATIONS 6.0 NUNERICAL CONSIDERATIONS In the preceding chapters. algorithms have been develOped for calculating with rational B-splines. If these algorithms are to be used on a computer. they must be implemented using floating point arithmetic. Floating point arithmetic with finite word lengths does not always yield accurate or even remotely usable results. 65 66 In this chapter some numerical considerations for implementing these algorithms on a computer will be discussed and some examples given. lhen scan line techniques [3] [4] [5] are used. the patch coordinates are transformed to screen coordinates. Many calculations are done in screen space. which we define with X horizontal. Y vertical. and 2 normal to the plane of the screen. Each pixel is one unit wide in X and one unit high in I. This is illustrated in Figure 6.0.1. The coordinate transformation from world space to screen space can present some problems because of the size of the window assigned to the screen; we will discuss this in terms of the 533; £55525. If the zoom factor is very small. the entire patch may be smaller than one pixel; if the zoom factor is very large. the smallest change that can be made to a parametric variable may represent a difference of more than one pixel. 67 2:1: ers ec : ion Cave x Viewing Screen 2 Figure 6.0.1 : Screen Space Coordinate System 6.1 BUILDING OF THE BLENDING FUNCTION MATRICES The blending functions are used for the evaluation of the surface and its various properties and for the transformation from NURBS to EURBS as well. 'hen building the matrices, care must be taken to avoid introducing numerical errors which change the surface definitions 68 during the transformation. A surface could also be inadvertently altered by using blending functions different from those used by the original designer. The uniformity of the knot vector yields a denominator which is common to each of the blending function matrices of a given order. 'hen the denominator is factored out of the matrix. the entries take an integer values. The matrix can thus be stored as integer values along with the denominator. This will avoids the roundoff or truncation which occurs if the matrix is stored as real numbers. 6.2 EVALUATION OF TEE COX-DEBOOR AND CURBS ALGORITHMS Roundoff error is usually thought of as a problem when the number of Operations becomes large. such as in inverting or reducing a large matrix. However. the examples given later show that roundoff error has a significant effect on the evaluated values of even low-order blending functions. The algorithm for the transformation from NURBS to EURBS requires that the NURBS surface and its derivatives and the blending functions for the EURBS surface be evaluated at specific points. To ensure that the resulting 69 EURBS surface is identical to the NURBS surface within the full precision available. these evaluations should be done as accurately as possible. 6.3 SIIULTANEOUS EQUATION SOLVER The transformation algorithm also requires the solution of a set of linear equations. The number of equations to be solved depends on the number of control points to be found. If only one control point needs to be found for a subsegment. then only one equation is generated. But if multiple control points are to be found. then a stable means must be used to solve the set of linear equations for the control points. Even for the relatively small number of equations to be solved. matrix inversion is not a good means to solve the system of equations because of the roundoff error which occurs and because of the problems which can arise if the matrix is ill-conditioned. There are many good numerical analysis textbooks such as [17] and [18] which deal with the subject in depth. An acceptable procedure is to use row operations on the augmented matrix to reduce the matrix to triangular form so that back substitution can be used to find the control points. 70 If increased precision is desired. an iterative method or iterative improvement can be used. 6.4 EXAIPLES OF EVALUATION PROBLEMS This section presents three examples of problems which occurred when implementing the CURBS algorithm. The central problem in each case was the evaluation of the blending functions. For each example. the order of the surface was 3 or 4 and the evaluation was done on a thirty-two-bit computer. 6.4.1 Evaluation of Normal Components To use scan-line techniques [3] [4] [5]. it is necessary to evaluate the silhouette edges as well as the physical edges of the patch. The silhouette edge of a patch is made up of points on the patch where the 2 component of the normal is zero. The evaluations are done on these edges at points where the edge crosses a scan line. This is called walking the edge. Valking the physical edges from scan line to scan line is relatively easy to do. 71 To walk the silhouette edge. a bi-variate Newton-Raphson method is used to iterate to each scan line that is crossed by the edge. The result of the Newton-Raphson calculation is a search direction and an approximation of the distance to go in that direction. Figure 6.4.1 is a plot of the 2 component of the normal as a function of the distance along the search direction. The basis function evaluation was done in single precision (24-bit mantissa). Depending on the error bounds used for zero. the point could be anywhere along this search direction. Because of the noisy nature of the 2 component. convergence is a problem. Figure 6.4.2 is a plot of the 2 component along the same search direction with the evaluation of the basis function done in double precision. In this case. the method converges in three iterations. ~aaaooaoo N q«O —03IOZ 72 x10'2 0.4 -n.4 [III I I I l ' l ' r 0.0 0.2 0.4 0.6 0.8 Search Direction Figure 6.4.1 : Single Precision Evaluation 73 x10'3 0.966 J l 53 I I l l l "3030'0300 N I 53 i? l l l l l —OBHOZ «.0 Ié 8 L0 4.970 . } . ' 1 ' I T I I I EM IL? 0u4 0J5 0J3 L0 Search Direction Figure 6.4.2 : Double Precision Evaluation 6.4.2 Normalization Of The Normal Calculating the 2 component of the normal alone requires less effort than computing all three components of the normal. However. this can cause problems with defining what is considered to be zero for the given component. In some cases the normal is changing so rapidly that it is 74 difficult to find a small value of the unnormalixed 2 component. This is shown in Figure 6.4.3. “Nggm Ems: Q / l l I J s 2‘: I J l ////J/V//// ///////// 7///////// / /////////// / ”soaovaoo N S 9 S 1 l l —03IOZ -O é t2 1 I -W////////// ‘n-azllt illllillllilll II 10 15 20 Iteration Step Figure 6.4.3 : Unnormalized 2 Component Calculation If all three components of the normal are calculated and the components normalized so the normal has length one. then what is considered to be zero is defined with respect to the value one. In this case. finding the zero point is not a problem, as shown in Figure 6.4.4. In this case the absolute value of the 2 component is between zero and one. 75 Using this method. acceptable limits about zero can be defined. 010'2 0.2 1411 naoacoaon N /////'/ [/// 0A //// / //////////// —03IOZ «O 5'9 l I Iteration Step Figure 6.4.4 : Normalized 2 Component Calculation 6.4.3 Face To Face The 2 component of the normal is also used to decide which face a point is on. If the 2 component is positive. then the point is on the plus face of the surface. If it 76 is negative. then the point is on the minus face. Silhouette edges are then the edges of a face. Even when the 2 component is normalized. if the blending functions are not calculated in double precision, the errors which occur result in misleading results. as shown in the Table 0.3129074EI03 0.7580551EI03 -O.6958499EI05 6.4.1. : S : T I Single Precision : Double Precision I I I I I 0.25 I 0.25 I0.5982137EI04 I 0.241537IEI02 I I_ _ ' I ==__ __ I I I I 0.00 l 0.25 0.00000000+00 l I I I I I I I I — _ I Table 6.4.1 : Example of Calculated Values of the 2 Component of the Normal The errors in the above table may not seem to be significant. and the results are more than adequate for such things as wire frame drawings. outlining patches or subpatches. or finding a pixel. The single precision values. however. are not acceptable for deciding which face a point is on. 77 Doing evaluations in double precision may require more computation time. as will computing all the components of the normal. But the additional effort is more than repaid when iterative methods converge quickly as apposed to taking many iterations or not converging at all. In some modern processors. all instructions are done in double precision and rounded to single precision when the results are stored. In this case the only cost associated with using the double precision is the additional storage for the results in that form. Numerical considerations of the implementation of evaluation algorithms have been discussed. and some examples of problems have been given. Sample output from the surface assesment package. ISU COLORSCOPE. which uses the algorithms developed in this dissertation. is given in the next chapter. In addition to the examples. conclusions are drawn about this work. CHAPTER VII EXAMPLES AND CONCLUSIONS 7.0 EXAMPLES AND CONCLUSIONS The preceding chapters develOped algorithms useful for the interactive shading of surfaces defined by rational B-splines. This chapter presents some examples of surfaces shaded using these algorithms. After the examples. some conclusions about this work are drawn. 78 79 7.1 Examples The following seven figures were generated using the NSU COLORSCOPE surface assessment package. These surfaces represent geometric as well as freeform surfaces. These figures are not intended to demonstrate the primary error detection function of the COLORSCOPE package. but rather that the package can handle the rational B-spline surface. The figures also demonstrate some of the capability of the rational B-spline. Figure 7.1.1 : Shaded Cube. Sphere. Cone. and Cylinder Figure 7.1.2 : Shaded Torus 81 Figure 7.1.3 : Shaded Automobile Hood Figure 7.1.4 : Curvature Shading of Automobile Hood 82 Figure 7.1.5 : Intermediate Stamping of a Styled Wheel Figure 7.1.6 : Shaded Bumper 83 Figure 7.1.7 : Shaded Turbine Blade Figure 7.1.8 : Absolute Maximum Curvature Shading of Turbine Blade 84 7.2 Conclusions This dissertation reviewed the rational B-spline. It defined a class of rational B-splines called Enhanced Uniform. The advantages of using this class of B-spline were discussed. The CURBS algorithm for calculating with Enhanced Uniform Rational B-splines was developed and presented. The performance of this algorithm was compared to the commonly used Cox-deBoor algorithm for B-spline evaluation. This comparison indicated that the CURBS algorithm would typically produce a three-to-one savings over the use of the Cox-deBoor algorithm. even for low-order B-splines. An apparent limitation of the CURBS algorithm. that the B-splines must be of the enhanced uniform rational B-spline type. was overcome by develOping the nonuniform-rational-B-spline-to-enhanced-uniform-rational- B-spline (NURBS-to-EURBS) conversion algorithm to convert any surface based on the NURBS knot vector to a mathematically equivalent surface based on the EURBS knot vector. The converted surface was shown to have a storage penalty but no time penalty for surface evaluation. In fact. the CURBS algorithm can be used on the converted surface so that there is a three-to-one savings in machine 85 operations. The CURBS and NURBS-to-EURBS algorithms combine to offer the tradeoff of lowered computation time for increased storage space. Some numerical consideration of the implementation of the algorithms deve10ped here were discussed. Finally. examples of shaded images produced by these algorithms were given. The algorithms deve10ped in this dissertation are a viable solution to the problem of rapidly rendering accurate shaded images of rational B-spline surfaces. They offer significant advantages over the prevalent algorithms now in use for this purpose. as well as for other tasks requiring calculations with rational B-splines. APPENDIX A CALCULATION EXAMPLE 86 87 A. RATIONAL B-SPLINE CALCULATION EXAMPLE This appendix presents examples of calculating a B-spline and a rational B-spline curve. The results of these calculations were used to generate the figures in Chapter Two of this dissertation. The control points used to define the B-spline curve are given in Table A.l. “011-2020 Each of the CPi is then multiplied by its th set Control point 0 1 2 3 4 5 6 7 8 9 10 11 Table A.l To define the weight is of rational homogeneous X IO.1490382E-Ol IO.1279431E+01 IO.1903875E+01 0.1243925E+01 0.1540538EI01 I0.7496327E+OO 0.6192721E+OO IO.8542724EIO2 IO.4414083E+OO 0.3807828E+OO O.1471788E+OO 0.1153550EIO2 Y 0.2372459E+01 0.2134099E+01 IO.1257925E+01 IO.7112664E+OO O.1074566E+01 IO.4256855E+OO IO.3516818E+OO 0.5969771E+OO I0.2523811E+OO IO.2515893E+OO 0.2455054E+OO 0.1825014E+OO Original Control Points assigned to control BIspline CEIVO. I. positive each control point 0P1. weight points CP;. control points are given in Table A.2. 0 HI‘ pic\oceq "3.4(t) N4’4tt) 102 I (tI1)°(tI1)/(3-1) + (4-tI'0/(4-2) = (t2 — 2-: + 1)/2 = (~:3 + 6-:2 - 12-: + 3)/4 = [(t-O)-(t2-4-t+4)/2]/(2-0) + [IS-t)-(~2-t2+6-t-3)/2]/(3-0) . (7-:3 - 36-:2 + 34-: - 13)/12 = [(t-O)-(-2-t2+6-t-3)/2]/(3-0) + [(4—:)-I:2-2-:+1I/2I/I4-1> = I-3-:3 + 12-:2 - 12-: + 4)/6 = [It-1)-(t2-2-t+1)/2]l(4-1) + (3-:)-0/(3-2) .. (:3 - 3-:2 + 3-: - 1)/6 Equations B.11. 8.12. B.13. and B.14 are then matrix form as: [Ni.4It). N2.4(t). "3'4(t), ”4'4(t)] 3 1/12 - [:3 :2 : 11 I" -3 7 -6 I 13 -36 24 I -36 34 -24 i 24 -13 3 written 05 (:-0)-o/(1-0) + [I2-:)-(:2-4-:+4)/21/(2-0) 3.11 B.13 in APPENDIX C KNOT VECTOR TRANSLATION 103 104 C. EFFECTS OF KNOT VECTOR TRANSLATION ON BLENDING FUNCTIONS Appendix B formulated the blending function matrix for the second subsegment of the fourth order knot vector 4 5 6 6 6 6 ] x = [ 0 0 0 0 l 2 3 i = 0 l 2 3 4 5 6 7 8 9 10 11 12 This appendix formulates the blending function matrix for the second subsegment of the translated knot vector X I [I1 -1 -l -l 0 1 2 3 4 5 5 5 5 ] i = 5 6 7 8 9 10 11 12 The effective knot vector for the second subsegment is 105 The translation has not changed Figure B.2 and the same ten Ni.k must be found. Using equations 2.1.2 and 2.1.3 :1:1 3' N4.1(t') = 1 0.1 N3'2(t') = (:'+1)-0/-o/(3-1) c.6 = I:')2 12 N1.4<:') = (t'+1)-0/(o+1) + [I1-t')-((:')2-2-(:')+1)/2I/(1+1) c.7 106 = (-(:')3 + 3-(:')2 - 3-(:') + 1)/4 N2,4(:'> - [I:'+1)-I(:'I2e2-(:')+1)/2I/(1+1) + [I2-:'I-<-2-(:'>2+2-I:')-1)/2I/(2+1) 0.3 a (7-(t')3 - 15«(t')2 + 3-(t') - 7)/12 N3,4I:'I - II:'+1>-(-2-I:')2+2-(:')-1)/2I/I2+1) + II3-:')-II:')2>/2I/<3-o> - 0.9 - (-3-I:')3 + 3-(:')2 + 3-(:') + l)/6 N4’4(t') - II:'—0)-(I:')2>/2J/I3-o) + (4—:')-0/(4-1) 0.10 .