INTERACTIVE LEAST SQUARES SURFACE FITTING ANTHONY HARM SAMSOM B.Sc. U n i v e r s i t y of B r i t i s h Columbia, 1978 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF THE FACULTY OF GRADUATE STUDIES (Department of Computer Science) We accept t h i s t h e s i s as conforming to the required standard. THE UNIVERSITY OF BRITISH COLUMBIA May, 1980 (c) Anthony Harm Samsom, 1980 by MASTER OF SCIENCE i n In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make i t freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department nf Co ^ p t>f 6v The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 n . t P W « W 2, /°i*0 D E - 6 B P 75-51 1 E i i ABSTRACT This thesis i s concerned with the design and implementation of a surface f i t t i n g package i n an interactive graphics environment. Surface f i t t i n g techniques are used to generate a smooth looking, easy to evaluate, bivariate function given a set of data points on some domain i n the plane, and are thus useful for a variety of applications. We consider the implementation of a surface f i t t i n g technique using weighted least squares with tensor products of B-splines on regular data grids ( i . e . the position of the data points can be represented as the cross product of two vectors). While somewhat more r e s t r i c t i v e than other surface f i t t i n g methods, t h i s technique, when applicable, i s extremely e f f i c i e n t . Knot placement and weight placement are discussed as methods of adapting the spline surface to rapidly varying regions on the domain. A disadvantage of the o r i g i n a l method used to solve for the coe f f i c i e n t s of the spline surface i s that the domain of the function to be approximated must be rectangular. An algorithm to extend the surface f i t t i n g method to non-rectangular domains, thus removing t h i s r e s t r i c t i o n , i s presented. An interactive surface f i t t i n g package i s provided, which allows a user to f i t a spline surface to a set of data points on a regular g r i d . This provides a powerful tool which may be used to e f f e c t i v e l y modify the spline surface and indicate the accuracy of the approximation. i i i TABLE OF CONTENTS Chapter 1 I n t r o d u c t i o n 1 Chapter 2 Survey o f methods f o r Surface approximation 7 Methods based upon p o i n t s 7 Methods based upon t r i a n g l e s 10 Methods based upon squares 15 Chapter 3 B-Spline Surface Approximation 18 One dimensional s p l i n e s 18 G e n e r a l i z a t i o n to two dimensions 23 Chapter 4 Adapting the s p l i n e s urface to q u i c k l y v a r y i n g data and Non-Rectangular Domains 27 Knot Placement 27 Weight Placement 31 Relaxing Domain R e s t r i c t i o n s 31 Chapter 5 R e s u l t s and E v a l u a t i o n 35 Chapter 6 Conclusions 48 B i b l i o g r a p h y 51 Appendix A A Surface F i t t i n g Package 53 The Command I n t e r p r e t e r 53 The Support Routines 61 Appendix B Generation and E v a l u a t i o n 62 TABLES Table I: E r r o r s on f u n c t i o n f l using evenly spaced knots 39 Table I I : E r r o r s on f u n c t i o n f2 using evenly spaced knots 40 Table I I I : E r r o r s on f u n c t i o n f2 using a good knot placement 41 Table IV: E r r o r s on f u n c t i o n f3 41 Table V: E r r o r s on f u n c t i o n f2 by methods from (Franke 79) 42 Table VI: Timings from generation and e v a l u a t i o n o f s p l i n e 43 Table V I I : Timings on s e l e c t e d methods from (Franke 79) 43 Table V I I I : E r r o r s from f u n c t i o n f4 over an even mesh 45 Table IX: S p l i n e approx. of f u n c t i o n f4 over an optimized mesh 46 Table X: S p l i n e approx. of f u n c t i o n f4 over an L shaped domain 47 FIGURES Figure 1: Max-min t r i a n g l e o p t i m i z a t i o n 13 Fi g u r e 2: Smooth s p l i n e 18 Figure 3: Banded matrix 24 Figure 4: Local t r o u b l e area 27 Fi g u r e 5: G l o b a l t r o u b l e area 28 Figure 6: Non-rectangular domain 32 Figure 7: Function f l (Saddle function) 36 Figure 8: Function f l ro t a t e d 90 degrees 36 Figure 9: Function f2 (Exponential t e s t function) 37 Fi g u r e 10: Function f2 (Rotated 90 degrees) 37 Figure 11: Function f 3 ( C l i f f t e s t function) 38 Figure 12: P l o t of f u n c t i o n f4 44 Figure 13: Function f4 rotated 90 degrees 45 v i ACKNOWLEDGEMENT I would l i k e to s i n c e r e l y thank Dr. U r i Ascher f o r h i s advice and f i n a n c i a l support i n the development and w r i t i n g of t h i s t h e s i s . 1 Chapter 1_ : I n t r o d u c t i o n Many problems a r i s e where the value of a b i v a r i a t e f u n c t i o n i s known o n l y a t some d i s c r e t e p o i n t s on i t s domain. Surface f i t t i n g techniques use t h i s data t o generate an a n a l y t i c f u n c t i o n t h a t approximates the unknown f u n c t i o n . A l t e r n a t i v e l y , i f the a n a l y t i c r e p r e s e n t a t i o n of the o r i g i n a l f u n c t i o n i s known but expensive t o evaluate, s u r f a c e f i t t i n g techniques may be used to generate a f u n c t i o n t h a t a c c u r a t e l y approximates i t but i s much cheaper to evaluate. Surface f i t t i n g techniques are used i n a v a r i e t y of f i e l d s ; f o r i n s t a n c e , i n Graphics these techniques can be used to generate compact, accurate and smooth re p r e s e n t a t i o n s of the surfaces to be d i s p l a y e d . The techniques can be s i m i l a r l y a p p l i e d i n such d i v e r s e areas as Geography, P h y s i c s , Computer Aided Design and i n d u s t r i a l m i l l i n g a p p l i c a t i o n s (Spath 7 4 ) . In t h i s t h e s i s we consider the implementation of one such technique using l e a s t squares f i t s w i t h tensor products of B - s p l i n e s on r e g u l a r data g r i d s . In g e n e r a l , there are three b a s i c b u i l d i n g blocks used f o r generating s u r f a c e s : P o i n t s , T r i a n g l e s and Squares. Methods using p o i n t s to b u i l d s u r faces are u s u a l l y based on some v a r i a n t of Shepard's formula ( B a r n h i l l 7 7 ) . The idea behind these methods i s to generate the value of an a r b i t r a r y p o i n t on the domain as a weighted average of the data p o i n t s around i t . Although Shepard's formula u s u a l l y generates a poor f i t , s e v e r a l v a r i a n t s ( B a r n h i l l 77) which improve on the approximant's c o n t i n u i t y e t c . g e n e r a l l y do q u i t e w e l l . Methods based upon p o i n t s are u s u a l l y simple to code, have minimal preprocessing requirements and need r e l a t i v e l y l i t t l e space t o represent the approximant; however, p o i n t e v a l u a t i o n i s g e n e r a l l y 2 quite expensive. Methods based upon triangles usually can provide good surface approximations. They have r e l a t i v e l y cheap point evaluation cost and reasonable preprocessing time. However they require a considerable amount of space to store the triangles and basis functions and are usually d i f f i c u l t to code. Methods based upon squares are the least f l e x i b l e of the three. They adapt poorly to non-rectangular domains and may also have trouble adapting to i r r e g u l a r l y spaced data. Although lookup may be faster on a rectangular g r i d than on a triangular one, the basis functions over rectangles are more complex than the corresponding ones over triangles so point evaluation may be more expensive. However, such methods do have d i s t i n c t advantages i n more re s t r i c t e d applications. I f the domain i s rectangular i n shape and the data points are given on g r i d intersection points, i t i s possible to b u i l d a tensor product piecewise polynomial surface over the domain. Using B-splines, the spline surface can be constructed e f f i c i e n t l y (de-Boor 79). The B-spline representation of the surface i s compact and the point evaluation cost i s r e l a t i v e l y inexpensive. This method of surface f i t t i n g can be considered as the generalization of one dimensional spline approximation to two dimensions using tensor products. A one dimensional polynomial spline approximation of order k (degree < k) may be represented as a lin e a r combination of l i n e a r l y independent basis functions where each basis function i s a piecewise polynomial of degree < k on the domain of the function. B-splines form such a basis with each basis function having a l o c a l support. I f the number of basis functions i s the same as the number of data points 3 ( i . e . interpolation) a l i n e a r system of equations can be solved to obtain the c o e f f i c i e n t s of the spline approximant (a detailed explanation i s given i n Chapter 3). Using the B-spline basis results in a banded matrix which can be solved quickly. I f there are more data points than basis functions a least squares f i t can be generated by solving the corresponding overdetermined system of equations. I t i s advantageous i n t h i s case to form the normal equations i n spite of conditioning problems. The matrix A has a zero structure which depends upon the placement of of the data r e l a t i v e to the knots specifying the B-splines and the order of the spline. Thus the normal equation matrix A TA has a single band whose width i s dependent only on the order of the spline. One dimensional interpolation and approximation can be extended to two dimensions using tensor products. (1.2) s(x,y) = I g i j Bi (x) Bj (y) where s(xj,yj) = ( ( f ( x i , y j ) , i = 1, .., M ) j = 1 , .. , N) This construction forces the data to l i e on a regular g r i d over a rectangular domain. Although the c o e f f i c i e n t s can be solved d i r e c t l y as i n the one dimensional case, i t i s possible to solve for the < c o e f f i c i e n t s by solving a series of one dimensional problems. This reduces the computational cost considerably (from 0((MN)3) to 0(M 3 + MN2 + NM2 + N 3) i f the matrices are f u l l ) . Both systems are banded i f the B-spline basis i s used but the d i r e c t approach results i n a large matrix with (2(k - 1) + 1) bands while the one dimensional approach produces 2 small matrices with one band each. Therefore, as discussed in Chapter 3, the d i r e c t approach produces a matrix that i s harder to solve e f f i c i e n t l y . 4 Again a least squares f i t can be generated i f there are more data points than basis functions. The co e f f i c i e n t s can be found by solving the system d i r e c t l y or (as i n the interpolation case) solving a series of one dimensional problems. In both cases i t i s advantageous to form the normal equations as A^A has a much nicer zero structure than the matrix A. I f the system i s solved d i r e c t l y the data may be somewhat scattered (there must be at least one data point i n the support of each basis function) and weights may be applied to each data point. However, i f the system i s to be solved as a series of one dimensional problems, the data must be on a regular g r i d and the weights may only be applied to these g r i d l i n e s instead of the individual data points. The one dimensional approach imposes some r e s t r i c t i o n s on the adaptability of the surface to trouble areas (rapidly changing data i n some region) and i t s a p p l i c a b i l t y to non-rectangular domains. Knot placement i s re s t r i c t e d to changing knot l i n e s (because of the way the basis functions are constructed). As a result the adaptability of the spline to trouble spots on the domain by changing the knots i s also l i m i t e d . There are some reasonable knot placement algorithms for one dimensional problems. An algorithm proposed by Dodson (Dodson 72) i s discussed and an extension to the two dimensional problem i s proposed. The method i s applied to two dimensions by having the user i n t e r a c t i v e l y choose a "scan l i n e " on the domain which the one dimensional knot selection technique i s applied to. The way that weights can be applied to the domain implies s i m i l a r problems i n using weights to force the spline surface to adapt to trouble areas. No method i s proposed for the automatic application of weights. Weights are probably most effec t i v e when applied c a r e f u l l y by the user i n an 5 interactive environment where the results can be c a r e f u l l y monitored. In both knot placement and the application of weights the use of a graphics display i n an interactive environment greatly improves the effectiveness of the method. This approach to surface f i t t i n g i s re s t r i c t e d to rectangular domains. Since many problems are posed over non-rectangular domains, i t would be nice to extend the method to include them. The approach suggested i s to f i t a spline surface over a rectangular domain which encloses the domain of the function to be approximated. This leaves holes i n the domain which should be " f i l l e d " with data that has a minimal effect on the f i t of the spline. The method proposed to generate t h i s data consists of two stages: a) F i l l i n g the holes using a Coons' patch defined on a rectangle around each hole i n the domain. b) Generating a (possibly weighted) surface over the domain and using the values of the surface over the holes to update the values generated by the Coons' patch. A surface f i t t i n g package that allows the user to i n t e r a c t i v e l y f i t a polynomial spline surface to a set of regularly spaced data points i s presented. The graphics c a p a b i l i t i e s provide a powerful tool a user can use to adapt the spline f i t and evaluate the results. The package i s divided into two sections: the command interpreter and the support routines. The command interpreter i s written i n Pascal and the support routines are i n Fortran. This guarantees p o r t a b i l i t y r e s t r i c t i o n s . However, as the package i s highly dependent upon UBC graphics routines, i t i s unreasonable to expect p o r t a b i l i t y anyway. The thesis i s divided into f i v e chapters and three appendices, 6 Chapter 1 being the introduction. Chapter 2 contains a survey of some of the popular surface f i t t i n g techniques. A development of the generalization of one dimensional splines to two dimensions i s presented i n Chapter 3. Chapter 4 contains two short sections: Adapting the spline by moving the knots and/or changing the weights and adapting the spline to non-rectangular domains. Chapter 5 contains some results from the surface f i t t i n g package and some limited comparisons with methods presented i n (Franke 79). The conclusions are stated i n Chapter 6. Appendix A contains a descriptions of the surface f i t t i n g package and Appendix B contains the source code for the two routines that generate and evaluate the spline surfaces. 7 Chapter 2 : Survey of methods for surface approximation The problem of surface approximation arises frequently i n practice. Given a set of data points = fCx^y^) , i = 1 .. N i n some domain D, one may wish to construct a smooth-loo king surface s(x,y) which approximates the function f on D. The approximation s(x,y) should be easy to store, evaluate, d i f f e r e n t i a t e , integrate etc. Even i f the function f i s known everywhere i n D, i f f i s complicated one may want to replace i t by the nicer function s. There are three basic building blocks for generating surfaces: POINTS, TRIANGLES and SQUARES. This chapter contains a b r i e f survey of some popular methods of generating surfaces using these building blocks with some discussion of the properties, advantages and disadvantages of each method. a) Methods based upon points. The basic idea behind Shepard's formula and i t s generalizations (Barnhill 77) i s to l e t the value of the approximating function at an ar b i t r a r y point (x,y) be defined by some weighted average of the data points i n D. The data points close to (x,y) influence s(x,y) more than points farther away. Shepard's formula can be written as: N 2 i = l (2.1) SF = N 2 l (x,y) j- ( x i f y i ) (x,y) = (Xi,yi) 8 Where = the distance from the i ' t h data point to (x,y). (2.2) d i = ^ ( ( x - X i ) 2 + ( y - y i ) 2 ) yu = constant > 0 This can be rewritten as (2.3) SF = N 2 Wi F i 1 N I Wi 1 where Wi = f l j=l The slope, smoothness and continuity of the surface depends upon u. I f 0 < u < 1 then there are cusps around the data points. I f u=l there are breaks (jump d i s c o n t i n u i t i e s i n f'(x,y) at ( x i , y i ) ) and for u > 1 there are f l a t spots where (2.4) S_f_ = 6_f_ = 0 d x d y B a r n h i l l i n (Barnhill 77) suggests u=2 as a good compromise (for an a s t h e t i c a l l y pleasing surface). One possible generalization of Shepard's formula i s to force the surface to interpolate the f i r s t p a r t i a l derivatives at each data point as well as the value at the data points: 9 (2.5) S]F • N y i=l L i F N Z 1/di 1 a? (x,y) (X-irYi) \ F i (x,y) = (Xi,yi) This function provides a continuous f i t to the data. Since every data point contributes to the value of s(x,y) the cost of one point evaluation of s(x,y) is quite high relative to other methods for surface approximation (for reasonable N). Franke and L i t t l e (Barnhill 77) have proposed a simple scheme to localize support. The scheme consists of replacing 1/dj in (2.1) or (2.5) by (2.6) (Ri -Ridi where Rj = radius of support around (x^yj) ( Ri - di (Ri - di) > 0 (Ri " <3i)+ = ( 0 otherwise The cost per point evaluation is s t i l l high (for a reasonable f i t ) but i f the number of data points in the domain is large, this considerably reduces the number of operations per evaluation. SF and S^F are rational functions. They do not however f i t polynomial functions very well: since d f = d f = 0 d x d y for SF, i t cannot model a linear function exactly. Using a similar argument i t can be shown that S^F cannot model a quadratic polynomial precisely. Various researchers have modified Shepard's formula to 10 increase function precision relative to polynomials. (Maclean 74) (Barnhill 77). In (Barnhill 77) Barnhill shows how to use the following theorem to increase function precision. The boolean sum PffiQ = P + Q-PQ has at least the interpolatory properties of P and the function precision of Q. Barnhill uses the theorem to construct two hybrid surface interpolants in his paper (Barnhill 77). P is replaced by Shepard's formula and Q is replaced by either a least squares f i t or a 9 degree triangular interpolant. In both cases, the form of Q did not place restrictions on the domain. The cost of increasing the function precision this way is high, however, and the choice of Q may place restrictions on the shape of the domain or on the position of the data points or both. Shepard's formula and most of i t s variants do not have any restrictions on the shape of the domain or placement of the data points and the methods are generally easy to implement. However, approximant evaluation is expensive as is increasing i t s precision (higher derivatives etc.). Information about the shape of the domain and i t s boundaries is not implicit to the method and must be kept separately. b) Methods based upon triangles. Using methods based on triangular elements to build a surface approximation is probably more complex to implement than other methods. Therefore, much more programming effort i s required to build 11 a s y s t e m u s i n g t h i s a p p r o a c h . M e t h o d s b a s e d u p o n t r i a n g l e s a r e u s u a l l y d i v i d e d i n t o t w o s t a g e s : 1) P r e p r o c e s s i n g w h e r e t r i a n g u l a t i o n o f t h e d o m a i n o p t i m i z a t i o n o f t h e t r i a n g l e s a n d g e n e r a t i o n o f t h e f u n c t i o n d e f i n e d o v e r e a c h t r i a n g l e a r e h a n d l e d . 2) F u n c t i o n e v a l u a t i o n w h e r e f o r a n a r b i t r a r y p o i n t ( x , y ) t h e t r i a n g l e w h i c h c o n t a i n s t h e p o i n t i s f o u n d a n d t h e f u n c t i o n d e f i n e d o v e r t h a t t r i a n g l e i s e v a l u a t e d t o o b t a i n s ( x , y ) . 1) P r e p r o c e s s i n g T h e r e a r e s e v e r a l a l g o r i t h m s a v a i l a b l e t o c o n s t r u c t t h e t r i a n g u l a r g r i d ( L a w s o n 7 7 ) ( F o w l e r a n d L i t t l e 79). L a w s o n s u g g e s t s t h e f o l l o w i n g a l g o r i t h m ( f o r D c o n v e x ) . a ) C h o o s e a p o i n t p * o n t h e b o u n d a r y o f D h a v i n g t h e s m a l l e s t X c o o r d i n a t e . b ) S o r t t h e r e m a i n i n g p o i n t s i n D i n i n c r e a s i n g E u c l i d i a n d i s t a n c e f r o m p * . L e t p * = p j a n d d e n o t e t h e r e m a i n i n g p o i n t s a s P2, P3, • • • ' P n • c ) T h e f i r s t e d g e i s d r a w n c o n n e c t i n g p ^ a n d P2. T h e n e x t p o i n t i n t h e s e q u e n c e n o t c o l l i n e a r w i t h p j a n d P2 i s c o n n e c t e d w i t h p ^ a n d P2 t o f o r m t h e f i r s t t r i a n g l e . I f t h e t h i r d p o i n t o f t h i s t r i a n g l e i s n o t P3 t h e n r e l a b e l i t a s p 3 a n d s h i f t a l l o t h e r s t o p ^ u p b y 1. d ) B u i l d a n i n i t i a l b o u n d a r y l i s t c o n s i s t i n g o f p ^ , P2, P3, a n d a s e c o n d o c c u r r e n c e o f p i a l o n g w i t h t h e i r a n g u l a r c o o r d i n a t e s . L e t C b e t h e c e n t r o i d o f t r i a n g l e ( p j _ , P2, p 3 ) ; t h e n t h e a n g u l a r c o o r d i n a t e s o f a p o i n t p ^ i s m e a s u r e d a s t h e a n g l e b e t w e e n t h e h a l f r a y s t a r t i n g a t p j t h r o u g h C a n d t h e h a l f r a y f r o m C t o p ^ . T h e 12 angle i s measured counter clockwise around C and i s always between 0 and 2 f l . The second occurrence of i n the boundary l i s t i s assigned the angular coordinate 2 f l e) For each point p^, k = 4, ..., N i) Determine the angular coordinate of p^. i i ) Find a pair of boundary points whose angular coordinates bracket t h i s angle. i i i ) Attach p^ to these points and place edge opposite pt of new tri a n g l e on a stack. iv) I f the stack i s not empty unstack one edge and apply the l o c a l optimization procedure to i t . I f the edge i s to be swapped, do so and place the two edges opposite p^ on the stack. Iterate u n t i l the stack i s empty. v) t r y to connect pj< to another neighboring boundary point. I f possible goto i i i otherwise k := k+1 and goto i . The processing cost i s estimated by Lawson to be 0(n^/3) + 0(n log n). There are several l o c a l optimization techniques available. The purpose of these i s to prevent any angles i n the triangular grid becoming small. I f some angles are allowed to become small t h i s causes accuracy problems i n the derivatives: See (Strang and Fix) or (Mitchell and Waite). A comparison of several l o c a l optimization methods i s given i n (Lawson 77). The max-min angle c r i t e r i o n i s discussed here. The basic idea i n the max-min angle c r i t e r i o n i s to increase the size of the smallest angle i n the triangular g r i d . This method s p e c i f i c a l l y attempts to maximize t h i s minimum angle l o c a l l y . 13 Consider the following figure: Figure 1 Max-min tr i a n g l e optimization I f the smallest angle i s increased by replacing the edge between and T2, the swap i s made creating two new triangles as shown above. Once the triangular g r i d has been b u i l t on the data points, generation of the functions defined over the tria n g l e can proceed. Since, for interpolation , the basis functions are usually chosen so that they are independent over each triangle ( i . e . the basis i s local) no system of equations need be solved to obtain the surface. As a r e s u l t , i f the functions are simple (linear etc.) the generation of the surface may be postponed u n t i l the value of s(x,y) i s asked for in that p a r t i c u l a r t r i a n g l e . The following i s an example of the construction of a piecewise l i n e a r surface over the triangular g r i d . The surface can be considered to be constructed of N piecewise l i n e a r basis functions N (2.7) s(x,y) = ? F i B-i(xfy) j=l where f 1 i = j B j ( X i , y i ) = . v 0 otherwise Any t r i a n g l e i n the gr i d w i l l have three non-zero basis functions defined on i t . These basis functions are defined on the canonical triangle as follows: 14 - u - v) So the value of s(x,y) on t h i s t r i a n g l e i s given by (2.9) s(x,y) = F]_ (1 - u - v) + F2 u + F3 v ((x,y) i s i n with values (Fl,-,F2,-,F3,-) at the vertices) Since any tr i a n g l e in the domain can be mapped onto the canonical triangle using the following transformation: (2.10) u = (x - X!) (y 3 - yi) - (y - yi) (x 3 - xj) det(J) v = (y - y 1) (xp - X ] ) - (x - xj) (y 2 - yi) det(J) where / ( X 2 " * i ) ( x 3 - xx) \ \(Y2 - yi) (y3 - y i ) / The value of any point i n the domain can be calculated by: a) Finding the tr i a n g l e which contains i t . b) Mapping the point onto the canonical triangle and using (2.3), where F]_, F2, F3 are the values of the function at the vertices of the t r i a n g l e . A number of functions suitable for surface approximation over triangles are given i n (Barnhill 77). To get continuity, f i r s t derivatives must be approximated. Since for higher degree CO and approximation the c o e f f i c i e n t s are expensive to generate, they are usually calculated and stored during the preprocessing stage. The cost per function evaluation depends on two factors: how fast 15 the tr i a n g l e which contains the point to be evaluated can be found and the complexity of the basis functions. A reasonable search time (Lawson 77) i s 0(log2 N ) 2 for an ar b i t r a r y point i n the domain. The piecewise l i n e a r surface involves around 15 multi p l i c a t i o n s per function evaluation plus search time. The cost per function evaluation i s generally considerably better here than most variations of Shepard's formula (for reasonable accuracy). I f not a l l the data points are used to generate the surface, t h i s approach can be made to adapt to l o c a l "trouble spots" in the domain. In (Fowler and L i t t l e 79) t h i s was done by f i r s t generating a triangular g r i d using only the l o c a l maxima and minima on the domain. Then i f the surface over a given tr i a n g l e did not meet i t s accuracy requirements a point was chosen i n t e r i o r to t h i s triangle and the trian g l e was replaced by three new tria n g l e s . This process i s repeated u n t i l each tri a n g l e meets i t s accuracy requirements. When compared to the piecewise l i n e a r f i t generated by using a l l the data on the rectangular g r i d good results were obtained using t h i s method i n terms of data compaction. I f the user has a package available or can afford to build one th i s i s the method of choice when working with large numbers of data points on ar b i t r a r y domains and/or the data points are not regularly spaced. c) Methods based upon squares, (rectangles) Rectangles are the least f l e x i b l e of the three methods. They have domain r e s t r i c t i o n s and are not e a s i l y adaptable to ir r e g u l a r l y spaced data. However i f the domain i s suitable and the 16 data points are regularly spaced (or can be moved to gr i d l i n e s using another method), t h i s method has d e f i n i t e computational advantages. Most interpolation methods defined over rectangles can be obtained from " t r a n s f i n i t e " Coons' patches developed by S. A. Coons i n (Coons 64). The b i l i n e a r l y blended Coons patch defined on the unit square i s given by: (2.11) s(u,v) = [l-u,u] + [F(u,0),F(u fl)] 1-v v F(0,v) F(l,v) - [l-u,u] TF(0,0) F(0,l )Ul-v F(1,0) F ( l , l ) where (0 < u,v < 1) This patch interpolates to the data a l l around the unit square. I f F(0,v) i s replaced by the d i s c r e t i z a t i o n (2.12) F(0,v) « (1-v) F(0,0) + v F(0,1) a b i l i n e a r interpolant results. S i m i l a r l y quadratic, cubic, and other interpolants can be generated using the appropriate d i s c r e t i z a t i o n . The d i s c r e t i z a t i o n of t h i s patch generate only C° interpolants however. Patches which when discretized give cl or better surfaces e x i s t . References may be found i n (Barnhill 77). If the domain of f(x,y) i s rectangular i n shape and the data points e x i s t on gr i d intersection points then i t i s possible to build a tensor product polynomial surface. Although the construction and solution of a li n e a r system of equations i s necessary, the l o c a l support offered by an appropriate basis such as B-splines results i n a system of equations which can be stored 17 and solved e f f i c i e n t l y . Using t h i s approach, obtaining C* continuity or better i s cheap as i s point evaluation of s(x,y) and storage of the surface representation. I t i s also possible to e f f i c i e n t l y generate s(x,y) as a weighted least squares f i t to data which i s usually d i f f i c u l t to accomplish using other approaches to surface f i t t i n g . / 18 Chapter 3: B-Spline surface approximation. This approach to surface approximation can be considered as the generalization of one dimensional spline approximation to two dimensions via tensor products. One dimensional polynomial splines are quite popular as they are both computationally e f f i c i e n t and quite f l e x i b l e for most applications. They are also well understood mathematically. a) One Dimensional splines: Let f l be a p a r t i t i o n a = t 0 < t]_ < ... < t M = b, then a (polynomial) spline function of order k (degree k-1) with N - 2 internal knots i s a function s(x) where s(x) i s a polynomial of degree < k on each i n t e r v a l ( ( t j ^ t i + j ) i=0 , M-l). s(x) i s said to be a smooth spline i f s(x) i s i n C k - 2 [ a , b ] . Let Pj < f ir denote the f i n i t e dimensional l i n e a r space of piecewise polynomials defined on [a,b] having breakpoints f l and having degree < k. Using the d e f i n i t i o n of a smooth splin e , i t i s f a i r l y easy to construct a spline interpolant to a given set of data. Example: Smooth cubic spline interpolation. 4^ ~te t i - - - - - — — "t M Figure 2 Smooth spline s(x) i s a piecewise cubic polynomial which means: (3.1) s(x) = si(x) = a i ( x - t i ) 3 + b i ( x - t i ) 2 + c i ( x - t i ) + d i 19 since s(x) i s i n C 2[a,b]: (3.2) 1) S i ( t i + 1 ) = s i + 1 ( t i + 1 ) = f ( t i + 1 ) 3 i => a-jhi + b^hj + Cjh£ + = d j + j i 2) s ' i ( t i + 1 ) = s ' i + 1 ( t i + 1 ) => 3ajhi + 2bjhj + c j = c j 3) s»i(t i + 1) = s " i + 1 ( t i + 1 ) => 6aihj + 2bj = 2bj where hj = (tj + j - t j ) This implies that we need to solve a system of equations with N+2 unknowns i n N equations to obtain a l l the c o e f f i c i e n t s . In practice, the two additional unknowns are usually provided by either setting s"(a) = s"(b) = 0 (natural spline) or setting s 1(a) and s 1(b) to some approximation of the slope at the end points of the i n t e r v a l . Note that the above equations form a tridiagonal system which i s to be solved for the c o e f f i c i e n t s CJ, i=l,M. This approach, although simple for t h i s example i s not e a s i l y applied to least squares problems or generalized to two dimensions. Al t e r n a t i v e l y the problem can be re-stated as M (3.3) s(Xi) = T a^B^x-i) i = 1,...,M j=l where each i s a piecewise polynomial i n P,g,7rand B^ i s i n C 2[a,b]. Solving the resulting system of equations for aj would generate s ( x ) . Clearly any l i n e a r l y independent polynomial basis i n P4, ir would s u f f i c e , but i f the basis functions are chosen c a r e f u l l y the system could become as computationally a t t r a c t i v e as the previous approach. A poorly chosen basis has, along with the added computation involved in solving a f u l l matrix, conditioning problems as w e l l . One basis 20 w i t h the d e s i r e d l o c a l support i s the B- s p l i n e b a s i s . (3.4) B i ( x ) = [ t i f t i + k ] ( t - x)+' where (t - x ) + = {max (t-x,0) } and [ t i , . . . ,ti+|<] f denotes the k-th d i v i d e d d i f f e r e n c e of f u n c t i o n f defined r e c u r s i v e l y here as (3.5) [ t i , t i + 1 ] f = f ( t i ) - f ( t i + 1 ) t i " t i + 1 [ t i , . . . , t i + ] J f = [ t j , . . . , t i + k - i 1 f - [ t j - n , . . . , t j + k ] f t i ~ ti+k I t i s easy t o show t h a t B i has l o c a l support, i f x > t i + k then (t - x ) + = 0 => B i ( x ) = 0 i f x < t i then (t - x ) ^ - ! i s a polynomial of degree k - 1 and s i n c e [ t i , t i + k ] f = f k ( f ) where 2, i s some p o i n t on the l i n e B i ( x ) = 0. Using L e i b n i z ' formula which s t a t e s t h a t i f f ( x ) =g(x)h(x) f o r a l l x then i+k (3.6) [ t i , . . . , t i + k ] f g = T [ t i , . . . , t r ] f [ t r , . . . , t i + k ] g r=i we can d e f i n e the recurrence r e l a t i o n (3.7) [ t i , . . . , t i + k ] f = x - t j [ t i , — , t i + k ] f t i+k ~ fci + fci+k ~ x C t i , — / t i + k ] g ti+k ~ fci which can be used to e f f i c i e n t l y generate B - s p l i n e s when needed f o r a r b i t r a r y k (de-Boor 79). B- s p l i n e s over an even mesh have some very n i c e p r o p e r t i e s which lead t o a w e l l conditioned system of equations to solv e (de-Boor 79). 21 For (k=4) the corresponding B-spline basis, B j ( x ) , on an even mesh i s given by S((x - Xj)/h) (h i s the distance between mesh points). / 0 (3.8) S(x) = (2 - x ) 3 - (1 - x ) 3 -24 6 (2 - x ) 3 - (1 - x ) 3 -x^ + (1 + x ) 3 4 6 24 (2 - x ) 3 -(2 - x ) 3 (1 - x ) 3 5 ~2T x < -2 -2 < x <-l -1 < x < 0 0 < x < 1 1 < x < 2 x > 2= 2 For an uneven mesh and a given k (3.7) can be expanded to form the basis functions e x p l i c i t l y . However each basis function cannot be generated by scaling a "canonical basis function" as i s done here for an even mesh. For an a r b i t r a r y k, equation (3.7) i s used i m p l i c i t l y to generate the basis functions. There are some problems defining the B-splines at the end points of the i n t e r v a l i n question using the current p a r t i t i o n . For the remainder of the thesis a p a r t i t i o n on [a,b] w i l l be defined as (3.9) ti < ... <tk = a < t k + 1 < ... t M = b< ... < t N + k Note that no more than k knots may coalesce or the system becomes rank de f i c i e n t . I f two knots coalesce at x j the spline approximation loses another degree of smoothness at that knot, i f k knots coalesce at X j there i s a jump discontinuity at that point. The knot sequence above applied to N data points s p e c i f i e s a NxN banded system of equations to be solved. Solving the system 22 (3.10) Bi (xi) B i ( x M ) B M ( x i ) BM (XM) — ' F l = for the c o e f f i c i e n t s (from 3.4) provides the interpolating spline s(x) This approach can be modified quite readily so that s(x) approximates the data instead of interpolating i t . Frequently i n applications there are more data points than B-splines (an overdetermined system of equations). l e t 1^ = number of data points and l e t Ljf, > M then for the least squares approximation, the system (3.11) (A TA)c = A Ty where c = [ai,a2,...,a^] T define the normal equations which can be solved quite e f f i c i e n t l y using a d i r e c t method for solving banded symmetric positive d e f i n i t e systems of equations. Since from (3.6) Bj has support over k intervals: (3.12) BiBj ^ 0 i f f |i - j | <" k the system i s banded (for k = 4 => bandwidth = 5) and diagonally dominant (for an even mesh). Forming the normal equations i n t h i s case, although providing a nice matrix, could cause some conditioning problems for uneven meshes as K(A^A) = K(A)2 (where K(B) denotes the condition number of a matrix B) and K(A) can get quite large i f an abruptly varying mesh i s used. A knot sequence that results i n a high condition number can be e a s i l y specified. I f the knot sequence i s spaced so that some B-splines have many data points i n t h e i r support and some have very few (and these 23 barely within the support of that B-spline), then there w i l l be large and very small diagonal elements i n ATA. However, the zero structure of the matrix A depends on the placement of the data and knots, so solving and storing i t as e f f i c i e n t l y as (ATA) i s d i f f i c u l t . This (for most problems) makes the solution using the normal equations quite a t t r a c t i v e . For B-splines the p r a c t i c a l experience indicates that an uneven mesh w i l l not cause accuracy problems. I t appears that K(A) i s a poor indication of what i s happening for t h i s system, b) Generalization to two dimensions Generalization to two dimensions using tensor products i s the next step i n the process. A two dimensional basis function i s constructed quite e a s i l y from B-splines using tensor products. (3.13) ( V i j ( x , y ) = Bi(x) Bj(y) Splines of d i f f e r e n t orders may be used i n the x and y d i r e c t i o n and we denote the orders of the spline i n the x and y directions as k x and ky respectively. A spline surface approximating the data could be characterized using t h i s basis as: MxN (3.14) s(x,y) = 7 g i j B i M B j t y ) i f j=l Note that when interpolation i s used, t h i s construction implies that the data l i e s on a "wire frame" g r i d over a rectangular domain. Obtaining the c o e f f i c i e n t s for s(x,y) using interpolation i s quite straightforward. There are MN equations i n MN unknowns so a l l that has to be done i s to solve a system of equations MN x MN i n s i z e . Because of the l o c a l support the B-spline basis o f f e r s , the matrix i s again banded: 2(k - 1)-1 bands 24 Figure 3_ Banded matrix k x and ky are assumed to be the same (=k) to simplify the description but because there i s more than one band, (the number of bands depend on k x and ky, and the distance between the bands depends on M and N ) , i t i s harder to store and solve the system e f f i c i e n t l y than i n the one dimensional case. I t i s possible, however, to solve t h i s system of eguations by solving a series of corresponding one dimensional problems i n each di r e c t i o n . Let G M X N be the matrix of c o e f f i c i e n t s (Gjj = g j j ) and l e t L be the matrix containing the values of the function to be interpolated (L^j = f ( x i , y j ) ) . The matrices A and B are the matrices corresponding to the B-splines on each axis ( i . e . A j j = B J ( X J ) and B i j = Bi(yj)) . From (3.14) (3.15) s(xk,yx) M N 2 2 i=i j=i 9 i jB i B j = L k l N M = 2 ( 2 9 i j B i ( x k ) ) Bj(y 2) i=l j=l N = 2 (A<R) . G^ 1)) BiCy-.) j=l row c o l (3.16) s ( x k , y i ) = (A G B T ) k l Solving A G B^ = L for G i s done i n two steps: 25 (3.17) 1) A X(J) = L(J) j = 1,...N => X = G B T =>XT = B GT 2) B G T 0) = X T(J) j = 1,...,M where x ( J ) i s the j ' t h column of X etc. Since A and B are banded (1 band) matrices of size M x M and N x N, the storage and computation requirements are reduced considerably over d i r e c t l y solving the system. De-Boor i n (de-Boor 79) estimates the complexity of the d i r e c t approach as 0((MN) 3) and the separated system as 0(M 3 + MN2 + NM2 + N 3) for a general basis ( f u l l matrices) but the difference (for small k x and ky) i s even larger using a l o c a l basis, as the second system can be solved and stored more e f f i c i e n t l y . We consider again the problem where there are more data points than approximating functions. I f the d i r e c t approach i s used, an MN x MN system of equations again has to be solved to obtain the coe f f i c i e n t s when the normal equations are used. We no longer have the s t r i c t constraints on the placement of the data, only the constraint that there must be a data point i n the support of each basis function (where the function i s non-zero). I t i s again advantagous to form the normal equations, A TA x = A Tb, i n spite of possible conditioning problems as the zero structure for A TA i s dependent only on M, N, and the support of the basis functions which makes the matrix easier to store and manipulate. The matrix i s banded with bands separated by either a distance of M or N depending on the ordering. I f the data points are on a rectangular g r i d , the problem can again be reduced to solving a series of one dimensional problems as i n 26 equation (3.17). (3.18) A?A G B B T = ATL B which can be solved (as f o r i n t e r p o l a t i o n ) i n 2 steps: 1) (A TA) xO) = (ATL)(J) j = l , . . . L n 2) (BT B) G T(3) = (BX T)(J) j =1,...,M which i s eq u i v a l e n t t o (3.16) given the same data. From (3.16) A TA and B TB are banded matrices (1 band each) a property which can a l s o be used to advantage i n c o n s t r u c t i n g A TA e f f i c i e n t l y . Although t h e o r e t i c a l l y there could be problems because o f the l a r g e c o n d i t i o n number t h a t t h i s system of equations could have w i t h an uneven mesh, no d i f f i c u l t i e s have a r i s e n i n p r a c t i c e . Weights can be added to t h i s process as i n the l a r g e r system (3.19) A TW XA G B W y B T = A TW XL B Wy B T where Wx and Wy are diagonal matrices s p e c i f y i n g the "confidence" i n each g r i d l i n e i n x and y. This i s c o n s i d e r a b l y d i f f e r e n t than ATWAX = ATWb i n the f u l l system where each weight a p p l i e s to a s i n g l e data p o i n t . 27 Chapter 4: Adapting the Spline Surface to Quickly Varying Data and Non-Rectangular Domains Least squares spline surfaces generated using the separation method described in Chapter 3 are quite inexpensive to construct and evaluate. However, they have some limitations on their a b i l i t y to adapt to "trouble" areas (where the surface to be approximated changes abruptly) and their a b i l i t y to adapt to non-rectangular domains. To change the shape of the spline surface either weights may be added or the knot spacing may be changed. We discuss the knot placement f i r s t , a) Knot Placement Changing the knot sequence in each direction can in some cases dramatically improve the surface approximation. In other cases, however, since only grid lines can be moved, changing the knot sequence may be ineffective. If we have a surface to be approximated which contains a single local trouble spot, knots can be added until the desired accuracy is reached in that reaion. — ( \ — \ \ 7 — Figure A: Local trouble area Notice that more than just the trouble spot receives added attention. However, i f the trouble spot is not local to one area and not parallel to one of the axes, i t may be impossible to adapt the mesh efficiently. 28 X Figure 5_: Global trouble area Only increasing the number of knots everywhere w i l l improve the accuracy i n t h i s example. This time the whole surface gets a great deal more attention, including areas where i t i s not needed. Despite the l i m i t a t i o n s , the approximation to many surfaces may be improved by a careful placement of knots. Since only knot " l i n e s " can be moved, i t i s possible to use one dimensional knot optimization methods along these l i n e s for certain problems. There are two basic approaches to knot placement i n one dimension: a) Optimal knot placement solved as a nonlinear optimization problem (Jupp 78). This i s expensive and the results are not much better than b) Choosing a good knot placement using an approximation of the error across the i n t e r v a l . One approach to finding a good knot placement i s outlined below, for more d e t a i l s see (Dodson 72). The distance from f(x) to s(x) in the L2 norm i s given by: b (4.1) I If - s| | p = ( J l f - s|P dx)Vp a 29 Dodson presents the following approximation to t h i s distance: b (4.2) | ! f - s l i p < C $ | f k ( t ) f dt a where cv = ^ [k + l ] ^ 1 Equalizing the distance from f(x) to s(x) across each subinterval so that ti+1 b (4.3) $|f k| f f dt = 1 $|f k f dt n - k t i a should provide a good knot placement. Dodson suggests the following as an estimate for f k ( t ) as i t i s unavailable for most applications: a) Generate a spline approximation s(x) to the data with some preliminary knot placement with the desired' number of knots (usually provided by the user). b) Di f f e r e n t i a t e s(x) k-1 times to obtain the piecewise constant function s k _ l ( x ) . c) J o i n the midpoints of each adjacent piecewise constant i n s k~^(x) to form a new piecewise l i n e a r function f k - 1 ( x ) . d) Di f f e r e n t i a t e f k - l ( x ) to obtain the piecewise constant function f k ( x ) . This approximation to f k ( x ) i s used i n (4.3) to generate the desired mesh. If the o r i g i n a l spline f i t s(x) approximates the function f(x) reasonably well then using the new knot placement obtained above w i l l 30 generally improve the f i t of the spline function. However i f s(x) does not f i t the data well to begin with, the new knot sequence may not be an improvement over the old and the error I|f - s l i p may grow. Applying t h i s or any other one dimensional knot selection technique to two dimensional splines i n an automatic manner i s d i f f i c u l t i f not impossible to do e f f e c t i v e l y . I f the problem has a l o c a l trouble spot (Figure 4) then i f the scan l i n e s for the one dimensional optimization are chosen cor r e c t l y , the knots suggested by the optimization technique (along these lines) w i l l improve the overall surface approximation. I f there i s more than one trouble spot i t i s not so obvious how to apply the optimization techniques and i f the s i t u a t i o n resembles Figure 5 then the techniques are useless. I t appears u n l i k e l y that an e f f i c i e n t and r e l i a b l e algorithm can be developed to choose good knot sequences for a two dimensional surface f i t t i n g problem with spline products. I t i s much more r e l i a b l e and e f f i c i e n t to have a user c o n t r o l l i n g the position of the knots and deciding when the f i t i s acceptable (or when i t s time to give up). In the package described i n Appendix A, a user can opt to apply the one dimensional optimization techniques to suggest knot sequences to be used for his problem. A l l that i s intended i s to provide a tool to suggest knot sequences which can be altered i n t e r a c t i v e l y by the user or (more l i k e l y ) completely ignored. The most powerful tool available to generate a good spline f i t i s the a b i l i t y to i n t e r a c t i v e l y a l t e r the spline and to see the r e s u l t s . 31 b) Weight Placement Weights can also be applied to the surface to a l t e r the shape of the spline approximation. They can be used i n two ways, either to force a closer f i t i n one area by increasing the weight i n that area or to ignore the effect of some "impossible" area ( s a c r i f i c i n g accuracy there to increase accuracy on the rest of the surface). Weights must be applied conservatively, however. Improperly applied they can cause severe conditioning problems i n the system which, unlike the large condition number caused by an irregular mesh, may have a severe effect on the accuracy of the computed s(x,y). As stated i n Chapter 3, weights cannot be applied to individual data points i f equation (3.19) i s used to obtain the c o e f f i c i e n t s . The weights can only be applied to g r i d l i n e s . This causes problems si m i l a r to those encountered with the knot placement. I f the problem resembles Figure 4 i t i s not possible to apply weight to the trouble spot without applying i t along the g r i d l i n e s outside the trouble area (a nonlocal weighting r e s u l t s ) . I f the problem i s not l o c a l as i n Figure 5, then the addition of weights w i l l not help the f i t . Providing a system that automatically generates the weights would not be very e f f i c i e n t . Allowing a user to i n t e r a c t i v e l y add, delete and change weights, however, i s a very e f f e c t i v e approach to the problem. c) Relaxing Domain Restrictions Because of the nature of our approximating function s(x,y), i t has a very severe r e s t r i c t i o n on the shape of the domain. Since the method has d i s t i n c t computational and/or space advantages over more 32 general surface f i t t i n g techniques, i t would be nice to be able to effici e n t l y adapt this method to non-rectangular domains. Unfortunately, for a general domain D, no such technique exists but i t is s t i l l possible to extend this method to nearly rectangular domains. For many problems, i f a region does hot conform to the shape of the basis functions, the general technique is to apply some transformation T to the the domain to map i t onto a more suitable domain. The method has a major drawback for our application however, as i t also moves the data points off the grid lines. The data points must therefore be resampled using a more general surface f i t t i n g method before our approach may be taken. The transformation, i f complex, could also add considerable point evaluation and preprocessing cost. If instead, a rectangle D' is superimposed over D, i t i s possible to generate a surface over D' providing that data is given in D' D and the data l i e s along mesh lines. The cost per point evaluation is increased however as checking whether or not a point (x,y) is in D must be done for each point evaluation. D' = D U R Figure (5: Non-rectangular domain The problem now is how to generate data points in R which w i l l not adversely affect the spline f i t on D. If i t were possible to weigh each data point individually, then R could be ignored but because of the way s(x,y) is generated weights can only be applied to 33 mesh l i n e s . Good data must be generated to f i l l R i f a reasonable f i t i s to be expected. To f a c i l i t a t e a good f i t the boundary between D and R should be smooth (possibly higher than C°). While there may be many di f f e r e n t methods of attempting t h i s goal we have chosen the following a) For each region R in D find the smallest rectangle which encloses R. b) Use a variant of a b i l i n e a r Coons' patch over t h i s rectangle to generate data inside R. c) Generate a (possibly weighted) surface s(x,y) over D. Use the values from s(x,y) to replace the data values i n R i f the o r i g i n a l f i t i s not s a t i s f a c t o r y (otherwise q u i t ) . d) Repeat (c) (possibly reducing the weights) u n t i l the i t e r a t i o n does not improve the f i t or the f i t i s satisfactory. If the region R i s not rectangular, a weighted surface approximation to reduce the effect of generated data on D may be undesirable since i t w i l l also affect the f i t along the boundary between D and R inside D. The rectangles enclosing R and the discretized Coons' patch are obtained i n the following manner. i) F i r s t every point on D' i s v i s i t e d , i f (xj,yj) i s not i n D then i t s position i s recorded i n two weight vectors at WXj and WYj. i i ) For each region the following i s done: a) I f one side of R l i e s on the boundary of the rectangle, data on that side i s given by the straight l i n e interpolation between the two end points on that side i n D just outside R. b) I f two adjacent sides are on the boundary of the rectangle, 34 the v e r t e x (corner of D') value i s f i r s t computed as the average of the two nearest corners of the r e c t a n g l e around R and (a) i s a p p l i e d t o both edges. c) The Coons' patch i s then d i s c r e t i z e d by l e t t i n g the f u n c t i o n values on the edge of the u n i t square be the data values on the r e c t a n g l e e n c l o s i n g R. The d i s c r e t i z e d Coons' patch i s then evaluated a t each p o i n t i n s i d e R. Generating data t h i s way has v a r i e d r e s u l t s . Since f o r a non-rectangular R there are r e a l f u n c t i o n values i n s i d e the r e c t a n g l e e n c l o s i n g R there may be a jump d i s c o n t i n u i t y at the boundary between R and D i f the r e a l data i n the area does not agree w i t h the Coons' patch approximation. I f the surface around R i s smooth the s p l i n e f i t may not be a d v e r s l y a f f e c t e d s i n c e the d i s c o n t i n u i t y should be s m a l l . However i f the surface moves q u i c k l y around R then the d i s c o n t i n u i t y w i l l probably be l a r g e and the o s c i l l a t i o n s caused by i t w i l l be hard to damp out. 35 Chapter 5_ : Results and Evaluation To i l l u s t r a t e the effectiveness of the interactive spline approximation, four test functions were chosen. The f i r s t three were taken from Franke who used these functions to test a large number of methods for interpolating scattered data. Some accuracy and timing results from his paper have been included. However, since the methods discussed i n (Franke 79) are for scattered data, d i f f e r e n t data points had to be used for the generation of our spline surface and consequently, comparisons should be made cautiously as the methods have very d i f f e r e n t areas of application. The fourth test function i s included to demonstrate the knot selection algorithm and the data insertion technique. The domain of the three test functions from (Franke 79) i s the unit square. To improve the scaling of the graphics, the sampled data was translated and scaled along the x and y axes to the square [-1,1] X [-1,1]. The three test functions from Franke are: Saddle function: fl(X,Y) = (1.25 + cos(5.4Y))/(6(1 + (3x-l) 2) Exponential function: f2(X,Y) = .75 EXP ( -(9X-2) 2 + (9y-2) 2)/4 ) + .75 EXP ( -(9X+l) 2/49 - (9Y+1)2/10 ) + .5 EXP ( -((9X+7) 2 + (9Y-3) 2)/4 ) - .2 EXP ( -(9X-4) 2 - (9Y-7) 2 ) C l i f f function: f3(X,Y) = 1/9 (TANH(9Y - 9X) +1) The test functions are plotted below. In Figures 7, 9 and 11 the X axis i s horizontal, p a r a l l e l to the page and the y axis i s pointed away from the reader. In Figures 8 and 10 the picture has been 36 Figure 7: Function f l (Saddle function) .37U88lJ8yUF.T0O .83327S297E-02 Figure 8: Function fl rotated 90 degrees The second test function contains two hill s and a depression: SURFACE FITTING . lOUt97S93Ef01 -.18500471 IE + 00 .939S97893E+ D ; ( 1 1 1 ! i i i 1 1 i Figure 9 Function f2 (Exponential test function) SURFACE FITTING . 10U197598E + 01 Figure 10 Function f2 (Rotated 90 degrees) 38 The t h i r d t e s t f u n c t i o n c o n t a i n s a c l i f f running d i a g o n a l l y across the domain: Figure 11 Function f3 ( C l i f f t e s t function) 39 The accuracy of the spline f i t for these examples i s measured using 625 points on the domain of the function. The f i v e worst errors (|f(xi,yi) - s ) ( x i , y j ) | ) and their positions are given as well as an approximation to the least squares error. The spline surface for each example i s generated using 2500 data points evenly spaced domain (unless otherwise specified). Table I : Errors on function f l using evenly spaced knots KNOTS ABS ERROR L.S. ERROR 10 X 10 1.08E-3 at (-.265,-1 ) 2.60E-4 1.06E-3 (-.265,-.837) 1.05E-3 (-.265,-.918) .935E-3 (-.265,-.755) •852E-3 (-.265,-.673) 30 X 30 1.43E-5 (-.265,-1 ) 3.46E-7 1.19E-5 (-.265,-.918) 1.19E-5 (-.265,-.837) .894E-5 (-.347,-1 ) .835E-5 (-.265,.959 ) 40 The second test function i s more d i f f i c u l t to f i t accurately and was used by Franke as his primary test function. This function i s examined more extensively here as seen i n the following results: The number of knots, the knot spacing and the number of data points are a l l varied. Table I I : Errors on function f2 using evenly spaced knots KNOTS ABS ERROR L.S. ERROR 10 X 10 (even) 2.85E-2 2.59E-2 1.89E-2 1.87E-2 1.85E-2 (-.102,.551 ) (-.204,.551 ) (-.102,.469 ) (-.184,.306 ) (-.265,.714 ) 4.26E-3 (225 data points) 10 X 10 (even) 2.87E-2 2.58E-2 1.94E-2 1.89E-2 1.84E-2 (-.102,.551 ) (-.204,.551 ) (-.102,.469 ) (-.265,.714 ) (-.184,.306 ) 4.43E-3 (100 data points = 0> interpolation) 10 X 10 (even) 8.76E-2 8.17E-2 7.30E-2 7.27E-2 6.82E-2 (-.102,.878 ) (-.183,.878 ) (-.204,.878 ) (-.102,.959 ) (-.183,.959 ) 1.11E-2 (2500 data points) 30 X 30 (even) 8.39E-5 6.70E-5 6.37E-5 6.16E-5 6.08E-5 (-.102,.551 ) (-.204,.551 ) (-.102,.469 ) (-.143,.551 ) (-.102,.633 ) 9.79E-6 41 With a set of good knots, which were chosen using the package, the improved result given i n Table I I I are obtained. Table I I I : Errors on function f2 using a good knot placement KNOTS ABS ERROR L.S. ERROR (2500 data points) 10 X 10 1.38E-2 1.36E-2 1.21E-2 1.18E-2 1.15E-2 (-.265,.551 ) (-.265,.632 ) (-.510,-.510) (-.204,.469 ) (-.204,.551 ) 3.52E-3 (225 data points) 10 X 10 1.37E-2 1.34E-2 1.21E-2 1.21E-2 1.19E-2 (-.265,.551 ) (-.265,.632 ) (-.204,.469 ) (-.510,-.510) (-.204,.551 ) 3.53E-3 Table IV l i s t s some error results from function f4. Notice that the optimal knot spacing over t h i s function i s an even mesh and the optimal weighting i s also an even d i s t r i b u t i o n . I t i s not possible for the spline to accurately f i t t h i s function without spending additional time on areas that do not need the attention. Table IV: Errors on function f3 KNOTS ABS ERROR 10 X 10 (even) 5.07E-3 ( .224,.612 ) 5.07E-3 ( .612,.224 ) 4.84E-3 ( .224,-.184) 4.84E-3 (-.184,.224 ) 4.72E-3 ( .469,.612 ) 30 X 30 (even) 1.95E-5 ( .142,.224 ) 1.91E-5 (-.673,-.591) 1.88E-5 ( .633,.551 ) 1.87E-5 ( .224,.142 ) 1.85E-5 (-.183,-.265) L.S. ERROR 1.84E-3 3.99E-6 Table V l i s t s some representative results from (Franke 79). The 42 results presented were obtained by testing the method on 100 scattered data points. Table V: Errors on function f2 by selected methods from (Franke 79) METHOD f l f2 f3 ABS. L.S. ABS. L.S. ABS. L.S. Methods based upon points: Mod. Quad. Shepard .0125 .00194 .0573 .0128 .0468 .00551 Shepard .273 .0620 Methods based upon triangles: Lawson .0565 .00359 .0951 .0124 .0280 .00448 Akima .0274 .00423 .0647 .0125 .0520 .00609 Methods based upon squares: Franke-TPS .0165 .00273 .0940 .0164 .0295 .00483 Hardy .00461 .00052 .0225 .00357 .0244 .00330 The modified quadratic Shepard's formula i s b a s i c a l l y equation (2.5) modified to have l o c a l support. Both Lawson's method and Akima's method use basis functions over a triangular g r i d . Franke's method uses l o c a l (Thin plate) functions over a rectangular mesh. Hardy's method i s somewhat d i f f e r e n t from the methods previously discussed. I t uses a l o c a l basis function ( l i k e the B-splines) but the support of the basis function i s a c i r c u l a r region. A large sparse system of equations ( l i k e the one which could result from using B-splines for scattered data) must be solved to generate the surface. 43 In Table VI we present some timing results from the generation and evaluation of our spline surface. Table VII presents the timing results given i n (Franke 79) for the methods presented i n Table V. The times are i n seconds. An Amdahl 470-V6-II was used to generate and evaluate the spline surfaces. An IBM 360-67 was used for a l l timings taken from (Franke 79). Table VI: Timings from generation and evaluation of spline SPLINE SURFACE PREPROCESSING EVALUATION (Generating the spline) (625 points) (2500 data points) (100 basis functions) 0.13 0.26 (900 basis functions) 0.19 0.26 (100 data points) (100 basis functions) 0.05 0.26 Evaluation speed could be increased considerably since the algorithm used to generate the B-splines to evaluate the spline surface i s for splines of a r b i t r a r y order and therefore quite slow. Table VII: Timings on selected methods from (Franke 79) METHOD PREPROCESSING EVALUATION (1089 data points) Mod. Quad. Shepard 8.6 15 Shepard 17 Lawson 1.8 1.7 Akima 2.2 0.8 Franke-TPS 2.7 6.5 Hardy 7.1 13 44 Even allowing for the difference i n speed between computers, the ef f i c i e n c y of the one dimensional spline approach to surface f i t t i n g i s c l e a r l y demonstrated. The fourth test function was included to show the adaptability of the spline using a good knot placement, and the a b i l i t y of the spline to handle non-rectangular domains. The domain of the f i r s t example of function f4 i s the square ([-1,1] X [-1,1]). The function i s given by: f4(X,Y) = 1/(1+25X2) COS (11/2 Y) This function i s one dimensional i n nature. The surface varies the most along the X axis near the o r i g i n . S U R F A C E F I T T I N G . 1 0 O 2 ' ? S 6 1 2 E T 0 1 . 1 2 U 7 6 S 3 4 9 E - 0 6 . 5 8 H 0 S 3 . 7 1 5 E + 01 I i ' !— Figure 12 Plot of function f4 45 S U R F A C E f I T T J N G . 1 0 0 2 7 S 8 1 2 E T 0 1 . I 2 4 7 6 3 3 U 9 E - 0 6 . 7 6 2 2 1 8 0 ' 9 ' 4 . R > 0 I I — i ; ! : ; • , , , Figure 13 Function f4 rotated 90 degrees The error results for a spline of even mesh over function f4 are given in Table V I I I . Table VIII: Errors from function f4 over an even mesh KNOTS ABS ERROR L.S. ERROR 13 X 13 (even) 1.33E-2 at ( .143,-.204) 3.71E-3 1.32E-2 ( .143,.612 ) 1.31E-2 ( .143,-.102) 1.29E-2 ( .143,.142 ) 1.28E-2 (-.204,-.204) A good knot sequence was chosen using routines based upon the one dimensional optimization technique presented by Dodson as described i n Chapter 4. 46 Table IX: Spline approx. of function f4 over an optimized mesh KNOTS ABS ERROR L.S. ERROR 13 X 13 (opt) 2.73E-3 ( .224,-.204) 7.74E-4 2.71E-3 ( .224,.612 ) 2.68E-3 ( .224,-.102) 2.65E-3 ( .224,.143 ) 2.62E-3 ( .224,-.184) The optimization routines w i l l not produce a better knot unless the f i t i s already quite good. In t h i s case, 13 knots were needed along the X axis before the routines would work. Table X presents the results of applying the techniques described in Chapter 4 to extend our spline method to non-rectangular domains. The domain of the function to be approximated i s defined by the user supplied Fortran function "RI" (Rl(x,y) i s true i f (x,y) i s i n the domain of f ( x , y ) , otherwise i t i s f a l s e ) . For t h i s example the domain of the function i s L-shaped with ([-1,0] X [-1,0]) as the "hole" in the rectangular domain. The surface varies very rapidly near (0,0). A more gentle surface would show much better results than t h i s example. Error estimates are given after the i n i t i a l patch and after each successive surface i s generated. The maximum errors tend to move out of the "hole" as the weights are increased. Weights are applied to scan l i n e s that "touch" the hole. As a result the weight on the data points i n the hole i s the product of the weight on the two l i n e s that intersect at that data point. Table X: S p l i n e approx. of f u n c t i o n f4 over an L shaped doma i KNOTS 13 X 13 (OPT) (PATCH WEIGHT 0.2) 13 X 13 (opt) (SMOOTH WEIGHT .2) 13 X 13 (opt) (SMOOTH WEIGHT .5) 13 X 13 (opt) (SMOOTH WEIGHT 1 ) ABS ERROR 1.16E-2 AT 1.14E-2 1.11E-2 1.06E-2 1.01E-2 5.53E-3 3.17E-3 3.16E-3 3.06E-3 3.01E-3 5.53E-3 3.04E-3 2.91E-3 2.86E-3 2.86E-3 5.53E-3 3.08E-3 2.85E-3 2.82E-3 2.78E-3 -.204,-.428) -.204,-.346) -.204,-.510) -.204,-.265) -.204,-.591) .959,-1 ) .224,-.265) .224,-.183) .224,-.346) .224,-.102) .959,-1 ) -.347, .142) .224,-.183) .224,-.102) .224,-.265) .959,-1 ) -.347, .142) .224,-.183) .224,-.102) .224,-.265) L.S. ERROR 1.67E-3 8.06E-4 7.69E-4 7.64E-4 48 Chapter 6 : Conclusions A method of generating a polynomial spline surface s(x,y) approximating a. function f(x,y) represented by a set of regularly spaced data points over a "nearly" rectangular g r i d has been presented. The surface, although "global" in nature, can be generated e f f i c i e n t l y for f a i r l y large sets of data points. The point evaluation cost i s low and can be reduced i f the exact degree of the spline surface i s known i n advance (the B-splines can be evaluated more e f f i c i e n t l y ) . The cost may be reduced s t i l l further at the expense of space i f the representation of the surface i s converted to a tensor product piecewise polynomial representation. As i t i s quite hard for the human eye to detect discontinuity i n derivatives higher than 2 i n a function, the surface approximant needs to be no more than C 2 to appear smooth. For lower degree f i t s (say k = 4 "cubic polynomials") the surface appears smooth and does not tend to o s c i l l a t e as much as surfaces based upon higher degree polynomials. When comparing t h i s method to other approaches of surface f i t t i n g , i t should be noted that other methods generally do not have a l l the r e s t r i c t i o n s t h i s method has but very few are able to approximate the data (least squares or otherwise) rather than interpolate i t . We have (perhaps unfairly) compared the method to methods for interpolating scattered data discussed i n (Franke 79). The better methods presented there f i t the test functions f l , f 2 , and f3 i n Chapter 5 well using one hundred data points. The spline f i t does almost as well (in terms of accuracy) as the better methods presented on t h i s amount of data with an even mesh and 100 basis functions (interpolation). Since a least-squares f i t i s used, i t i s 49 possible to increase the number of data points without increasing the number of basis functions. As the number of data points i s increased the accuracy of the spline f i t i s also increased and becomes comparable with the best results i n (Franke 79) while the computational cost i s s t i l l very reasonable. The surface f i t t i n g package was used to c a r e f u l l y choose a good knot sequence over the function f2. The accuracy of the surface using the c a r e f u l l y chosen knot sequence was considerably better than the f i t generated using an even knot sequence as seen from the the results i n Chapter 5. The cost of generating the spline surface and evaluation (even with a r e l a t i v e l y large number of points) i s nominal compared to the cost of generating and evaluating the more general surface f i t t i n g methods. The r e s t r i c t i o n s imposed on the data and the domain imply that the data can be stored i n « 1/3 the space i t would take to store the ( x i , y i , Z j ) t r i p l e t s . S i m i l a r l y , the B-spline representation of the surface approximant i s very compact; there are two vectors specifying the location of the surface patches (knots for the B-splines) and s l i g h t l y more than one c o e f f i c i e n t per patch i s needed to specify the spline surface over these patches. In terms of space and speed t h i s makes the spline f i t very a t t r a c t i v e r e l a t i v e to the more general surface f i t t i n g techniques. As discussed i n Chapter 4, there are ways of extending the method to non-rectangular domains. I f there are space or speed constraints on a problem or the data has errors i n i t then i t may be j u s t i f i e d to re-sample a problem defined with scattered data by using a more general surface f i t t i n g method and using the least squares spline f i t over the re-sampled data. The interactive environment supplied by the surface f i t t i n g 50 package provides a powerful tool to the user: A user can dynamically a l t e r the f i t of the spline surface u n t i l he i s s a t i s f i e d with the results (the graphics f a c i l i t i e s allowing him to v i s u a l i z e and therefore judge the accuracy of the f i t ) . There probably i s not much more work that can be done to improve the behavior of polynomial spline product surfaces. There may be some promise i n surfaces generated using rational splines or non-polynomial functions but we believe the current research into t h i s approach i s well developed. There i s considerable opportunity for research into more general two dimensional surface approximation techniques. In conclusion the graphics package provided i n t h i s thesis e f f e c t i v e l y solves two dimensional surface f i t t i n g problems provided that the constraints on the data are met. 51 BIBLIOGRAPHY B a r n h i l l , Robert E. and Riesenfeld R.E., Computer Aided Geometric Design. New York: Academic Press, 1974. B a r n h i l l , Robert E., "Representation and Approximation of Surfaces". Mathematical Software I I I , ed. Rice, J . R., London: Academic press, 1977, pp 71-120. de-Boor, C a r l , "Least Squares Cubic Spline Approximation II - Variable Knots". Purdue University: Technical Report CSD TR 21, 1968. de-Boor, C a r l , "Package for Calculating with B-splines".SIAM J . Numer. Anal.. 1977, pp 441-473. de-Boor, C a r l , A P r a c t i c a l Guide to Splines. New York: Springer Verlag, 1979. de-Boor, Carl and Rice, John R., "An Adaptive Algorithm for Multivariate Approximation Giving Optimal Convergence Rates". University of Wisconsin: MRC-TSR-1773, 1977. Coons, S. A., "Surfaces for Computer Aided Design". Design D i v i s i o n , Mech. Engin. Dept., MIT, 1964, revised 1967. Dodson, David S., "Optimal Order Approximation by Polynomial Spline Functions". Purdue University: Ph.D. thesis, 1972. Fowler, R. J . and L i t t l e J . L., "Automatic Extraction of Irregular Network D i g i t a l Terrain Models", ACM SIGGRAPH, 1979, pp 199-207. Franke, R., "A C r i t i c a l Comparison of Some Methods for Interpolation of Scattered Data". Naval Postgraduate School: NPS-53-79-003, 1979. Jupp, D. L. B., "Approximation to Data With Free Knots". SIAM J . Numer. Anal.. 1978, pp 328-343. Lawson, C. L., "Software for C 1 Surface Interpolation". Mathematical Software I I I , ed. Rice, J . R., London: Academic press, 1977, pp 161-193. McLain, D. H. "Drawing countours from a r b i t r a r y data points". The Computer Journal. 1974, pp 318-324. Mi t c h e l l A. R. and Waite R., The F i n i t e Element Method i n P a r t i a l D i f f e r e n t i a l Equations. New York: John Wiley & Sons, 1977. Shepard D., "A Two Dimensional Interpolation Function for Ir r e g u l a r l y Spaced Data". ACM P r o c . 23rd Nat. Conf., 1965, pp 517-523. Spath, H., Spline Algorithms for Curves and Surfaces. Winnipeg: U t i l i t a s Mathematica inc., 1974. 52 Strang, G. and G. F i x , An A n a l y s i s of the F i n i t e Element Method. New J e r s e y : P r e n t i c e - H a l l , 1973. 53 Appendix A : A Surface F i t t i n g Package The code for the package i s divided into two sections: The command interpreter and the support routines. a) The Command Interpreter The command interpreter i s designed to provide an interactive environment which a user can f i t a tensor product polynomial spline surface to his data. Its command structure i s designed to protect the user as much as i s p r a c t i c a l while trying not to "get i n the way". A description of the basic data structures and commands follow: The Data Structures There are b a s i c a l l y four d i f f e r e n t "data objects" the user may manipulate i n the package. These are: 1) Data: The data that s p e c i f i e s the surface consists of L - A matrix of dimension at least LM x LN which contain the values of F^ at the grid intersection points. XX - A vector of dimension at lease LM containing the x coordinates. XY - The corresponding vector i n the y di r e c t i o n . WX - A vector specifying the weights to apply to the g r i d l i n e s specified i n XX. WY - The corresponding vector i n the y di r e c t i o n . LM - the number of gr i d l i n e s perpendicular to the x axis. LN - the number of g r i d l i n e s perpendicular to the y axis. LMDIM - the f i r s t dimension of the matrix L. 54 2) Spline Surface: The components that specify the spline surface approximant. GAMMA - A matrix of dimension at least M x N containing the co e f f i c i e n t s of the B-spline surface approximant. TX - The knot sequence specifying the B-splines i n the x di r e c t i o n . TY - The corresponding vector i n the y di r e c t i o n . KX - The order of the B-splines i n the x di r e c t i o n . KY - " " y di r e c t i o n . M - The number of B-splines i n the x di r e c t i o n . N - " " y di r e c t i o n . MDIM - The f i r s t dimension of GAMMA. 3) Graphics: The internal 3-D graphical object (IMAGE) X - A vector specifying the position of the g r i d l i n e s perpendicular to the x axis. Y - The corresponding y vector. Z - The value of the surface at the points X 0) Y NX - The number of g r i d l i n e s i n x. NY - The number of gr i d l i n e s i n y. NZ - The f i r s t dimension of the matrix Z. 4) Graphics: The IG representation of the surface (DISPLAY) (data available i n surface package) THETA - The angle the surface i s viewed from (Azimuth). PHI - The height (Elevation) surface i s viewed from. 55 The system provides commands for generating, changing, i n i t i a l i z i n g , reading and writing these objects. I t also t r i e s to keep the system consistent (e.g. a user may not generate the coe f f i c i e n t s of a spline approximant before the data has been read i n and the knots have been specified). Any non-ambiguous subset of a command may be used to specify that command (e.g INIT), and i n most cases l i n e boundaries are ignored when reading the commands. The commands available are: INITIALIZE I n i t i a l i z e mode allows the user to i n i t i a l i z e the following variables: ORDER KX KY I n i t i a l i z e s the order of the splines i n the x and y di r e c t i o n . I t invalidates the current spline surface i f there i s one. KNOTS M N I n i t i a l i z e s the number of B-Splines to be used i n the system. The knots are i n i t i a l i z e d to an even mesh with M B-splines i n the x dir e c t i o n and N i n the y. The knots are r e - i n i t i a l i z e d every time the bounds may change ( i . e . i f new data i s read i n ) . DISPLAY NX NY I n i t i a l i z e the number of (evenly spaced) mesh l i n e s i n the display. The display l i n e s are also i n i t i a l i z e d every time the bounds may change. 56 REGION i Indicate which user supplied function i s to be used to specify the domain of the function to be approximated. I f the domain i s rectangular no function needs to be specified I f a region has been previously specified then invoking the command "REGION" with no region specified cancels the region (the knots specify a rectangular region i m p l i c i t l y ) . This mode 4 internal sections which allow the user to change the KNOTS, DATA, WEIGHTS and DISPLAY data structures i n t e r a c t i v e l y . ADD K Add K knots to the current knot sequence. The user chooses a knot by positioning the cross hairs (Tektronix) at the place desired on the knot sequence displayed on the screen and pressing anything but <ENTER>. DELETE K Delete K knots from the current knot sequence. OPTIMIZE This mode allows the user to obtain a good knot placement (for a one dimensional problem along a scan line) using the knot placement algorithm as discussed i n Chapter 4. END To return to the top l e v e l . CHANGE 57 PICK Interactively choose the scan l i n e the optimization technique i s to be applied to. TRY I SAVE K Apply the optimization technique to the current scan l i n e . Save the suggested mesh i f SAVE i s specified. ROTATE Rotate the viewing d i r e c t i o n 90 degrees (change the axis p a r a l l e l to the screen). Return to KNOT mode. ROTATE Rotate the viewing axis. END Return to CHANGE mode. This mode f a c i l i t a t e s the generation of data as described i n Chapter 4. These commands do not invoke any graphics a c t i v i t i e s . Use a Coons' patch to patch the holes i n the domain. A region must already have been specified i f t h i s command i s invoked. END DATA PATCH 58 SMOOTH fWEIGHT WX Generate a spline approximant over the domain with weight WT on the scan l i n e s that pass through holes i n the domain. F i l l i n the holes with data generated using the spline. Weight mode allows the user to i n t e r a c t i v e l y change the weights applied to the (lines of data) over the surface. WEIGHT WT Specify the value of the weight to be inserted. CHANGE K Interactively change the weight to WT on K data l i n e s . ROTATE Rotate to other axis. END Return to CHANGE mode. Interactively change the grid spacing on the internal graphical object (IMAGE). I t only allows the user to change the g r i d i f a spline surface i s to be displayed as the data i s discrete and cannot be sampled a r b i t r a r i l y . ADD K Interac t i v e l y add K l i n e s on the display. This does not immediately affect the graphical object, screen or image. END Return to CHANGE mode. WEIGHT DISPLAY DELETE K Interactively delete K l i n e s to the display. REFRESH Refresh the screen adding the new display l i n e s . ROTATE Rotate the display on the screen adding the new display l i n e s . END Return to CHANGE l e v e l . END Return to top l e v e l . REFRESH/DISPLAY Draw screen using current (IG) graphical object. ROTATE Rotate the display on the screen. ( COEFFICIENTS IMAGE K READ \ KNOTS ( FROM filename DATA Read the corresponding item from the f i l e filename. The formats for the f i l e s are contained i n the routines RDATA etc. the I/O section of the subroutine l i b r a r y . For "IMAGE" K must be i n (0 .. 3) where 0 means the current display (SCREEN) which i s used to generate the IG graphical object. 60 WRITE J j COEFFICIENTS IMAGE K KNOTS DATA SCREEN ONTO filename / DEVICE= 'CALC' V 'PRINT1 Write the corresponding item onto the f i l e filename. COEFFICIENTS GENERATE ERROR DISPLAY IMAGE FROM f SPLINE ^ SURFACE V ERROR Generate the corresponding item. MOVE(IMAGE KITO SCREEN IMAGE L SCREEN Move the corresponding item (K in (1 .. 3)). SUBTRACT|IMAGE KjFROM ( IMAGE L 1 1 SCREEN j SCREEN Subtract the corresponding units placing the r e s u l t i n the second entry. END Terminate the run. 61 b) The Support Routines The support routines are divided into f i v e main sections: GRAPHICS Routines that generate the graphical object, display i t on the screen, read positions for display l i n e s and knots o f f of the screen etc. I/O Routines that handle the input and output of DATA, KNOTS, COEFFICIENTS, & UNITS (Data that specifies enough information to generate a graphical object). These routines should be used to prepare the data and knots when using t h i s package to generate a spline surface. OPTIMIZATION Routines that handle the one dimensional optimization. PATCHING Routines find the rectangles on the rectangle D1 that enclose the regions not i n the domain of the function to be approximated and generate the data to f i l l these "holes" using a Coons' patch. GENERATION AND EVALUATION The routines that generate the spline surface, evaluate i t and provide some error estimates. 62 Appendix B : Generation and Evaluation These two routines use the spline package i n (de-Boor 77) to generate and evaluate a tensor product polynomial spline surface of arbit r a r y order. C C "SPLINE GENERATION AND EVALUATION" C C SUBROUTINE LSCOEF(GAMMA, TX, TY, KX, KY, M, N, NDIM, XX, XY, 1 WX, WY, L, LM, LN, LMDIM, A, X) C C CALCULATE THE COEFFICIENTS "GAMMA" FOR THE C B-SLINE SURFACE APPROXIMATION C c GAMMA - COEFFICIENT MATRIX RETURNED. c TX - KNOT SEQUENCE ALONG X AXIS c TY - KNOT SEQUENCE ALONG Y AXIS c KX - ORDER OF SPLINE ALONG X AXIS c KY - ORDER OF SPLINE ALONG Y AXIS c M - NUMBER OF B-SPLINES IN X DIRECTION c N - NUMBER OF B-SPLINES IN Y DIRECTION c NDIM - SIZE OF GAMMA, AND ONE SIDE OF X (NDIM > MAX(M, c XX - X VALUES OF DATA POINTS. c XY - Y VALUES OF DATA POINTS. c WX - WEIGHTS ALONG X DIRECTION. c WY - WEIGHTS ALONG Y DIRECTION. c L - DATA MATRIX. c LM - NO. OF DATA POINTS ALONG X AXIS c LN - NO. OF DATA POINTS ALONG Y AXIS c LMDIM - FIRST DIMENSION OF DATA MATRIX c A - WORK SPACE FOR GENERATING SPLINE c X - WORK SPACE FOR GENERATING COEFFICIENT MATRIX. c c REAL TX(1), TY(1), XX(1), XY(1), A ( l ) , L(LMDIM,LN) , 1 X(NDIM,LN), GAMMA(NDIM,NDIM), VNIKX(20), WX(1), WY(1) C C ZERO MATRICES C MNMAX = M IF ( N .GT. MNMAX) MNMAX = N KMAX = KX IF (KY .GT. KMAX) KMAX = KY IMAX = KMAX * MNMAX DO 10 1=1,IMAX A(I) = 0. 10 CONTINUE 63 DO 20 1=1 ,M DO 20 J = 1, LN X(I, J ) = 0. 20 CONTINUE C ILEFT = KX IMK = 0 DO 70 L I = 1, LM C C FIND INTERVAL "TX(ILEFT) <= XX(LI) < TX(ILEFT+1)" C 30 IF (ILEFT .EQ. M) GO TO 40 IF (XX(LI) .LT. TX(ILEFT + 1)) GO TO 40 ILEFT = ILEFT + 1 IMK = ILEFT - KX GO TO 30 C C PLACE VALUE OF ALL NON ZERO BSPLINES EVALUATED ON XX(LI) C IN VNIKX C 40 CALL BSPLVN(TX, KX, 1, XX(L1), ILEFT, VNIKX) C C T C GENERATE A L AND PLACE IN X C DO 60 J J = 1, KX I = IMK + J J WXBSPL = WX(L1) * VNIKX(JJ) DO 50 J = 1, LN X(I,J) = X( I , J ) + WXBSPL * L(L1,J) 50 CONTINUE C C T C GENERATE A A AND PLACE IN A C DO 60 Ml = J J , KX J = IMK + Ml U K = IJ(J,I,KX) A(IJK) = A(IJK) + WXBSPL * VNIKX (Ml) 60 CONTINUE 70 CONTINUE C C C T (J) T (J) C SOLVE A A X = A L TO OBTAIN THE MATRIX X C RATIO = 1.0E-10 EPS = 1.0E-77 NRHS = 1 ITER = 0. NSCALE = 0 DO 90 I = 1, LN CALL FBAND(A,X(1,I),M,KX,NRHS,RATIO,DET,JEXP,NSCALE, 1 ITER,EPS) 64 NRHS = 2 90 CONTINUE C C CLEAR MATRICES C DO 95 I=1,IMAX ' A(I) = 0. 95 CONTINUE DO 100 I = 1, MNMAX DO 100 J = 1, MNMAX GAMMA(I,J) = 0. 100 CONTINUE C ILEFT = KY IMK = 0 DO 150 L I = 1, LN C C FIND THE INTERVAL "TY (ILEFT) <= XY(L1) < TY (ILEFT+1) 1 1 C 110 IF (ILEFT .EQ. N) GO TO 120 IF (XY(L1) .LT. TY(ILEFT + 1 ) ) GO TO 120 ILEFT = ILEFT + 1 IMK = ILEFT - KY GO TO 110 C C OBTAIN VALUE OF ALL NON ZERO B SPLINES OF XY(L1) C 120 CALL BSPLVN(TY, KY, 1, XY(L1), ILEFT, VNIKX) C C T T C GENERATE B X AND PLACE IN GAMMA C DO 140 J J = 1, KY I = IMK + J J WYBSPL = WY(L1) * VNIKX(JJ) DO 130 J = 1, M GAMMA(I,J) = GAMMA(I,J) + WYBSPL * X(J,L1) 130 CONTINUE C C T C GENERATE B B AND PLACE IN A C DO 140 Ml = J J , KY J = IMK + Ml U K = IJ(J,I,KY) A(IJK) = A(IJK) + WYBSPL * VNIKX(Ml) 140 CONTINUE 150 CONTINUE C C T T(J) T T C SOLVE B B GAMMA = B X TO OBTAIN GAMMA AND PLACE IN GAMMA C NRHS = 1 RATIO = 1.0E-10 65 DO 170 I = 1, M CALL FBAND(A,GAMMA(1,1),N,KY,NRHS,RATIO,DET,JEXP,NSCALE, 1 ITER,EPS) NRHS = 2 170 CONTINUE CALL COPY(GAMMA,TX,TY,KX,KY,M,N,NDIM) RETURN END C REAL FUNCTION SPL2D(X, Y) C C THIS FUNCTION EVALUATES THE B-SPLINE FUNCTION DEFINED BY THE C COEFFICIENT MATRIX GAMMA AND THE KNOT SPQUENCES TX AND TY C COMMON /CA/ GAMMAT(30,30), TX(35), TY(35), M, NDIM, N, KX, KY REAL BCOEF(20) C C SPL2D(X,Y)=SUM (SUM GAMMA(I,J) B (Y)) B (X) C I J J,K,T I,H,S C CALL INTERV(TX, M + 1, X, LEFTX, MFLAG) SPL2D1 = 0 . IF (MFLAG .NE. 0) GO TO 20 C C INNER LOOP FIRST C DO 10 J = 1, KX BCOEF(J) = BVALUE(TY,GAMMAT(1,LEFTX - KX + J),N,KY,Y,0) 10 CONTINUE C SPL2D = BVALUE(TX(LEFTX - KX + 1),BCOEF,KX,KX,X,0) 20 CONTINUE RETURN END
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Interactive least squares surface fitting
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Interactive least squares surface fitting Samsom, Anthony Harm 1980
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Interactive least squares surface fitting |
Creator |
Samsom, Anthony Harm |
Publisher | University of British Columbia |
Date Issued | 1980 |
Description | This thesis is concerned with the design and implementation of a surface fitting package in an interactive graphics environment. Surface fitting techniques are used to generate a smooth looking, easy to evaluate, bivariate function given a set of data points on some domain in the plane, and are thus useful for a variety of applications. We consider the implementation of a surface fitting technique using weighted least squares with tensor products of B-splines on regular data grids (i.e. the position of the data points can be represented as the cross product of two vectors). While somewhat more restrictive than other surface fitting methods, this technique, when applicable, is extremely efficient. Knot placement and weight placement are discussed as methods of adapting the spline surface to rapidly varying regions on the domain. A disadvantage of the original method used to solve for the coefficients of the spline surface is that the domain of the function to be approximated must be rectangular. An algorithm to extend the surface fitting method to non-rectangular domains, thus removing this restriction, is presented. An interactive surface fitting package is provided, which allows a user to fit a spline surface to a set of data points on a regular grid. This provides a powerful tool which may be used to effectively modify the spline surface and indicate the accuracy of the approximation. |
Subject |
Computer graphics |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-03-22 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0051811 |
URI | http://hdl.handle.net/2429/22271 |
Degree |
Master of Science - MSc |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1980_A6_7 S35.pdf [ 3.22MB ]
- Metadata
- JSON: 831-1.0051811.json
- JSON-LD: 831-1.0051811-ld.json
- RDF/XML (Pretty): 831-1.0051811-rdf.xml
- RDF/JSON: 831-1.0051811-rdf.json
- Turtle: 831-1.0051811-turtle.txt
- N-Triples: 831-1.0051811-rdf-ntriples.txt
- Original Record: 831-1.0051811-source.json
- Full Text
- 831-1.0051811-fulltext.txt
- Citation
- 831-1.0051811.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0051811/manifest