INTERACTIVE LEAST SQUARES SURFACE FITTING by ANTHONY HARM SAMSOM B.Sc. U n i v e r s i t y o f B r i t i s h C o l u m b i a , 1978 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES (Department o f Computer S c i e n c e ) We a c c e p t t h i s t h e s i s a s c o n f o r m i n g to t h e r e q u i r e d s t a n d a r d . THE UNIVERSITY OF BRITISH COLUMBIA May, 1980 (c) Anthony Harm Samsom, 1980 DE-6 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. I t i s 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.tP B P 75-51 1 E W«W 2, /°i*0 ii ABSTRACT This t h e s i s i s concerned with the design and implementation of a surface f i t t i n g package i n an i n t e r a c t i v e graphics environment. Surface f i t t i n g techniques are used t o generate a smooth l o o k i n g , easy to evaluate, b i v a r i a t e function given a s e t of data points on some domain i n the plane, and are thus u s e f u l f o r a v a r i e t y o f a p p l i c a t i o n s . We consider the implementation of a surface f i t t i n g technique using weighted l e a s t squares with tensor products of B-splines on regular data g r i d s ( i . e . the p o s i t i o n of the data points can be represented as the cross product of two v e c t o r s ) . 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 a p p l i c a b l e , i s extremely e f f i c i e n t . Knot placement and weight placement are discussed as methods o f adapting the s p l i n e surface t o r a p i d l y varying regions on the domain. A disadvantage of the o r i g i n a l method used t o solve f o r the c o e f f i c i e n t s of the s p l i n e surface i s that the domain o f the f u n c t i o n to be approximated must be rectangular. An algorithm to extend the surface f i t t i n g method t o non-rectangular domains, thus removing t h i s r e s t r i c t i o n , i s presented. An i n t e r a c t i v e surface f i t t i n g package i s provided, which allows a user t o f i t a s p l i n e surface t o a s e t o f data points on a regular g r i d . This provides a powerful t o o l which may be used t o e f f e c t i v e l y modify the s p l i n e surface and i n d i c a t e the accuracy o f the approximation. iii TABLE OF CONTENTS Chapter 1 Introduction 1 Chapter 2 S u r v e y o f methods f o r S u r f a c e a p p r o x i m a t i o n 7 Methods based upon p o i n t s Chapter 3 Chapter 4 7 Methods based upon t r i a n g l e s 10 Methods based upon s q u a r e s 15 B-Spline Surface Approximation 18 One d i m e n s i o n a l s p l i n e s 18 G e n e r a l i z a t i o n t o two d i m e n s i o n s 23 Adapting the s p l i n e surface t o q u i c k l y varying d a t a and Non-Rectangular Domains 27 Knot Placement 27 Weight Placement 31 R e l a x i n g 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 Bibliography Appendix A Appendix B 51 A S u r f a c e F i t t i n g Package 53 The Command I n t e r p r e t e r 53 The S u p p o r t R o u t i n e s 61 G e n e r a t i o n 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 u s i n g e v e n l y spaced knots 39 Table I I : E r r o r s on f u n c t i o n f 2 u s i n g e v e n l y spaced knots 40 Table Table Table Table Table I I I : E r r o r s on f u n c t i o n f 2 u s i n g a good knot placement IV: E r r o r s on f u n c t i o n f 3 41 V: E r r o r s on f u n c t i o n f 2 by methods from (Franke 79) V I : T i m i n g s from g e n e r a t i o n and e v a l u a t i o n o f s p l i n e V I I : T i m i n g s on s e l e c t e d methods from 41 (Franke 79) 42 43 43 T a b l e V I I I : E r r o r s from f u n c t i o n f 4 over an even mesh 45 Table 46 Table IX: S p l i n e approx. o f f u n c t i o n f 4 over an o p t i m i z e d mesh X: S p l i n e approx. o f f u n c t i o n f 4 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 Figure 2: Smooth spline 18 Figure 3: Banded m a t r i x 24 Figure 4: L o c a l t r o u b l e a r e a 27 Figure 5: G l o b a l t r o u b l e a r e a 28 Figure 6: N o n - r e c t a n g u l a r domain 32 Figure 7: F u n c t i o n f l (Saddle f u n c t i o n ) 36 Figure 8: F u n c t i o n f l r o t a t e d 90 degrees 36 Figure 9: F u n c t i o n f 2 ( E x p o n e n t i a l t e s t f u n c t i o n ) 37 F i g u r e 10: F u n c t i o n f 2 (Rotated 90 degrees) 37 F i g u r e 11: F u n c t i o n f 3 ( C l i f f t e s t f u n c t i o n ) 38 F i g u r e 12: P l o t o f f u n c t i o n f 4 44 F i g u r e 13: F u n c t i o n f 4 r o t a t e d 90 degrees 45 vi ACKNOWLEDGEMENT I would l i k e to sincerely thank D r . U r i A s c h e r f o r h i s a d v i c e and f i n a n c i a l s u p p o r t i n t h e development and w r i t i n g o f t h i s t h e s i s . 1 C h a p t e r 1_ : I n t r o d u c t i o n Many problems a r i s e where the v a l u e known o n l y at techniques use approximates some d i s c r e t e p o i n t s the representation evaluate, this that of instance, d i s p l a y e d . The as milling fitting techniques analytic may be be In this such t e c h n i q u e be of similarly 74). expensive to i s much can a cheaper to to fields; generate i n such Design we to generate surfaces applied thesis using used the that analytic a r e used i n a v a r i e t y o f Computer A i d e d (Spath one used i t but smooth r e p r e s e n t a t i o n s can fitting function i s known but these techniques Geography, P h y s i c s , implementation of function approximates techniques applications an Surface A l t e r n a t i v e l y , i f the techniques i n Graphics and i t s domain. generate original fitting compact, a c c u r a t e areas the on function. accurately evaluate. Surface for to unknown surface function data of a b i v a r i a t e function i s and to be diverse industrial consider l e a s t squares the fits with t e n s o r p r o d u c t s o f B - s p l i n e s on r e g u l a r d a t a g r i d s . In general, there generating surfaces: points build to Shepard's f o r m u l a are three Points, surfaces basic Triangles are ( B a r n h i l l 77). and usually The building of the data u s u a l l y generates a points poor around Squares. based on used Methods some variant variants Shepard's (Barnhill 77) improve on the a p p r o x i m a n t ' 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 upon p o i n t s preprocessing requirements and represent approximant; however, the are using of the domain as a weighted i t . Although f i t , several Methods based for i d e a b e h i n d t h e s e methods i s t o g e n e r a t e t h e v a l u e o f an a r b i t r a r y p o i n t on average blocks u s u a l l y simple need which well. t o code, have m i n i m a l relatively point formula little evaluation is space to generally 2 q u i t e expensive. good surface evaluation Methods based upon t r i a n g l e s u s u a l l y can approximations. cost and They reasonable require a considerable have relatively preprocessing time. amount of space to store the provide cheap point However they t r i a n g l e s and b a s i s functions and are u s u a l l y d i f f i c u l t to code. Methods based upon squares are the l e a s t f l e x i b l e of the three. They adapt poorly to non-rectangular to irregularly domains spaced and data. may also Although have lookup trouble may rectangular g r i d than on a t r i a n g u l a r one, faster on a the basis functions over rectangles are more complex than the corresponding so point evaluation may be adapting be more expensive. ones over t r i a n g l e s However, such methods do have d i s t i n c t advantages i n more r e s t r i c t e d a p p l i c a t i o n s . I f the domain i s rectangular i n shape and the data points given on g r i d i n t e r s e c t i o n p o i n t s , i t i s p o s s i b l e to b u i l d a product piecewise polynomial surface over the 79). The point B-spline representation of the surface evaluation surface fitting dimensional cost can spline is be relatively considered approximation to the two Using (de-Boor i s compact and inexpensive. as tensor domain. B - s p l i n e s , the s p l i n e surface can be constructed e f f i c i e n t l y are This method g e n e r a l i z a t i o n of dimensions using the of one tensor products. A one dimensional (degree < k) may polynomial be represented s p l i n e approximation of as a l i n e a r combination of independent b a s i s functions where each basis function i s a polynomial of degree < k on the domain of the order k linearly piecewise f u n c t i o n . B-splines form such a basis with each b a s i s function having a l o c a l support. I f the number of b a s i s functions i s the same as the number of data points 3 (i.e. i n t e r p o l a t i o n ) a l i n e a r system obtain the coefficients of the of equations can be solved to spline approximant (a detailed explanation i s given i n Chapter 3). Using the B-spline b a s i s r e s u l t s in a banded matrix which can be solved q u i c k l y . I f there are more data p o i n t s than b a s i s f u n c t i o n s a l e a s t squares f i t can be generated by s o l v i n g 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 s p i t e of conditioning problems. The matrix A has depends upon the placement of of the data a zero structure relative which to the knots s p e c i f y i n g the B-splines and the order of the s p l i n e . Thus the normal equation matrix A A has a s i n g l e band whose width i s dependent only on T the order of the s p l i n e . One dimensional i n t e r p o l a t i o n and approximation can be extended to two dimensions using tensor products. s(x,y) = I (1.2) g i j B i (x) Bj (y) where s ( x j , y j ) = ( ( f ( x i , y j ) , i = 1, .., M ) j = 1 , .. , N) This c o n s t r u c t i o n 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 p o s s i b l e to solve f o r the < c o e f f i c i e n t s by s o l v i n g a s e r i e s of one dimensional problems. This reduces the computational cost considerably (from 0((MN)3) to 0(M MN 2 + NM 2 + N) 3 3 + i f the matrices are f u l l ) . Both systems are banded i f the B-spline b a s i s i s used but the d i r e c t approach r e s u l t s i n a l a r g e 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 l e a s t squares f i t can be generated i f there are more data points than b a s i s functions. The c o e f f i c i e n t s can be found by s o l v i n g the system d i r e c t l y or (as i n the i n t e r p o l a t i o n case) s o l v i n g a s e r i e s of one dimensional problems. In both cases i t i s advantageous t o form the normal equations as A^A has a much n i c e r zero s t r u c t u r e 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 a t l e a s t one data p o i n t i n the support o f each b a s i s function) and weights may be applied t o each data point. However, i f the system i s t o be solved as a s e r i e s o f one dimensional problems, the data must be on a regular g r i d and the weights may only be applied t o these g r i d l i n e s instead of the i n d i v i d u a l data p o i n t s . The one dimensional approach imposes some r e s t r i c t i o n s on the a d a p t a b i l i t y o f the surface t o trouble areas ( r a p i d l y changing data i n some region) and i t s a p p l i c a b i l t y t o non-rectangular domains. Knot placement i s r e s t r i c t e d t o changing knot l i n e s (because of the way the b a s i s functions are constructed). As a r e s u l t the a d a p t a b i l i t y o f the s p l i n e t o trouble spots on the domain by changing the knots i s a l s o l i m i t e d . There are some reasonable knot placement algorithms f o r 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 interactively i s applied choose a t o two "scan dimensional knot s e l e c t i o n dimensions line" technique by having on the domain which i s applied the user the one t o . The way that weights can be applied t o the domain implies s i m i l a r problems i n using weights t o force the s p l i n e surface t o adapt t o trouble areas. No method i s proposed f o r the automatic a p p l i c a t i o n o f weights. Weights are probably most e f f e c t i v e when applied c a r e f u l l y by the user i n an 5 i n t e r a c t i v e environment where the r e s u l t s can be c a r e f u l l y In both knot placement and the a p p l i c a t i o n o f weights monitored. the use o f a graphics d i s p l a y i n an i n t e r a c t i v e environment g r e a t l y improves the e f f e c t i v e n e s s of the method. This approach t o surface f i t t i n g i s restricted t o rectangular domains. Since many problems are posed over non-rectangular domains, i t would be nice t o extend the method to include them. The approach suggested i s t o f i t a s p l i n e surface over a rectangular domain which encloses the domain of the f u n c t i o n t o 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 s p l i n e . The method proposed t o generate t h i s data c o n s i s t s o f 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 o f the surface over the holes t o 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 s p l i n e surface t o a s e t of r e g u l a r l y spaced data points i s presented. The graphics c a p a b i l i t i e s provide a powerful t o o l a user can use t o adapt the s p l i n e f i t and evaluate the r e s u l t s . The package i s d i v i d e d i n t o two s e c t i o n s : the command i n t e r p r e t e r and the support routines. support The command routines are interpreter i n Fortran. i s written This i n Pascal guarantees and the portability r e s t r i c t i o n s . However, as the package i s h i g h l y dependent upon UBC graphics r o u t i n e s , i t i s unreasonable t o expect p o r t a b i l i t y anyway. The t h e s i s i s d i v i d e d into f i v e chapters and three appendices, 6 Chapter 1 being the i n t r o d u c t i o n . Chapter 2 contains a survey of some of the popular surface generalization presented fitting techniques. o f one dimensional i n Chapter 3. Chapter splines A development of the t o two dimensions 4 contains two short Adapting the s p l i n e by moving the knots and/or changing and adapting contains limited the s p l i n e some r e s u l t s comparisons conclusions t o non-rectangular from with are stated the surface methods presented i n Chapter 6. Chapter package i n (Franke Appendix A sections: the weights domains. fitting is 5 and some 79). The contains a d e s c r i p t i o n s o f the surface f i t t i n g package and Appendix B contains the source code f o r the two routines that generate and evaluate the s p l i n e surfaces. 7 Chapter 2 : Survey of methods f o r surface approximation The problem of surface approximation p r a c t i c e . Given a set of data p o i n t s some domain D, s(x,y) which one may arises = fCx^y^) frequently i n , i = 1 .. N i n wish to construct a smooth-loo king surface approximates the function f on D. The approximation s(x,y) should be easy to s t o r e , evaluate, d i f f e r e n t i a t e , integrate etc. i f f is Even i f the function f i s known everywhere i n D, complicated one may want to replace i t by the n i c e r f u n c t i o n s. There are three b a s i c b u i l d i n g blocks f o r 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 b u i l d i n g blocks with some d i s c u s s i o n of the p r o p e r t i e s , advantages and disadvantages of each method. a) Methods based upon points. The b a s i c idea behind Shepard's formula and i t s generalizations ( B a r n h i l l 77) i s to l e t the value of the approximating f u n c t i o n at an a r 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 p o i n t s c l o s e to (x,y) influence s(x,y) more than points f a r t h e r away. Shepard's formula can be w r i t t e n as: N 2 i=l N (2.1) SF = (x,y) j- (x yi) i f 2 l (x,y) = (Xi,yi) 8 Where = the distance from the i ' t h data point to (x,y). (2.2) di = ^((x-Xi) 2 + (y-yi) ) 2 yu = constant > 0 This can be r e w r i t t e n as N (2.3) 2 1 SF = Wi F i N I 1 Wi where Wi = f l j=l The slope, smoothness and c o n t i n u i t y of the surface depends upon u. I f 0 < u < 1 then there are cusps around the data p o i n t s . 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_ d x = 6_f_ = 0 d y B a r n h i l l i n ( B a r n h i l l 77) a s t h e t i c a l l y pleasing One surface suggests u=2 as a good compromise (for an surface). possible g e n e r a l i z a t i o n of Shepard's formula i s to force to i n t e r p o l a t e the first p a r t i a l d e r i v a t i v e s at point as w e l l as the value at the data points: the each data 9 N y Li F a? i=l (2.5) S]F Z 1/di 1 • \ F i (x,y) = This function provides a Since (X-irYi) (x,y) N every data (Xi,yi) continuous f i t to the data. point contributes to the value of s(x,y) the cost of one point evaluation of s(x,y) i s quite high r e l a t i v e to other methods Little for surface (Barnhill support. The (2.6) approximation 77) have (for proposed a reasonable simple scheme consists of replacing 1/dj N). Franke scheme to i n (2.1) or and localize (2.5) by (Ri Ridi where Rj = radius of support around ( Ri - d i (x^yj) (Ri - di) > 0 (Ri " <3i)+ = ( 0 otherwise The cost per point evaluation i s s t i l l high (for a reasonable f i t ) but i f the number of data points i n the domain i s large, t h i s considerably reduces the number of operations per evaluation. SF and polynomial for SF, S^F are rational functions. They do not however f i t functions very w e l l : since d f = d f = 0 d x d y i t cannot model a l i n e a r function exactly. Using argument i t can be shown that S^F cannot model a quadratic precisely. Various researchers have modified Shepard's a similar polynomial formula to 10 increase function precision relative to polynomials. (Maclean 74) (Barnhill 77). In (Barnhill 77) B a r n h i l l shows how to use the following theorem to increase function p r e c i s i o n . 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 interpolants i n h i s paper to construct two hybrid surface (Barnhill 77). P i s replaced by Shepard's formula and Q i s replaced by either a least squares f i t or a 9 degree triangular interpolant. In both cases, restrictions on the precision t h i s way domain. The the cost form of Q did not of i s high, however, and increasing the choice the place function of Q may place r e s t r i c t i o n s on the shape of the domain or on the p o s i t i o n of the data points or both. Shepard's formula and most r e s t r i c t i o n s on the points methods are and the shape of approximant evaluation of the i t s variants domain or generally i s expensive as (higher d e r i v a t i v e s e t c . ) . Information and i t s boundaries i s not implicit easy do not placement to have of the implement. i s increasing any data However, i t s precision about the shape of the domain to the method and must be kept separately. b) Methods based upon t r i a n g l e s . Using methods based on approximation is probably triangular elements to build more complex to implement a than surface other methods. Therefore, much more programming e f f o r t i s required to build 11 a system using this usually divided two Preprocessing 1) the triangles triangle are triangle over triangulation generation evaluation which that based upon triangles are stages: where and Methods of of the the domain function optimization of defined over each point (x,y) the handled. Function 2) 1) into approach. where contains triangle is the for point evaluated an is to arbitrary found obtain and the function defined s(x,y). Preprocessing There are triangular grid the following a) several (Lawson algorithm Choose a algorithms point 77)(Fowler (for p* D on available and Little to construct Lawson 79). the suggests convex). the boundary of D having the smallest X coordinate. b) Sort from p*. •• • c) the the ' The Let first sequence to is not d) Build second be form = P3 edge pj then an is not the and in D denote the first in increasing the remaining starting at i t pi of point pj with as pj the and 3 l i s t along C and and is P2 third shift Euclidian distance points P2, (pj_, as the the half next connected point of of P3, with to then angle up by and P3, the between from C and triangle p^ P2, in p^ coordinates. p3); ray point this p^, angular P2, as The P2. a l l others their i s measured through p^ consisting with triangle p^ and If p boundary of centroid a connecting triangle. relabel initial of drawn collinear occurrence coordinates ray p* points Pn• P2 C remaining to 1. a Let angular the p^. half The 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 is assigned the angular coordinate 2 fl e) For each point p^, k = 4, ..., N i ) Determine the angular coordinate of p^. ii) Find a p a i r of boundary points whose angular coordinates bracket t h i s angle. i i i ) Attach p^ to these p o i n t s and place edge opposite pt of new t r i a n g l e on a stack. iv) I f the stack i s not empty unstack one local o p t i m i z a t i o n procedure to edge and i t . I f the edge apply the i s to be swapped, do so and place the two edges opposite p^ on the stack. I t e r a t e u n t i l the stack i s empty. v) t r y to connect pj< to another neighboring boundary p o i n t . I f p o s s i b l e 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 s e v e r a l l o c a l o p t i m i z a t i o n techniques purpose of these becoming i s to prevent s m a l l . I f some angles any angles are allowed i n the and Waite). methods i s given in A comparison of (Lawson 77). The The triangular grid to become small this causes accuracy problems i n the d e r i v a t i v e s : See (Mitchell available. (Strang and Fix) or several l o c a l optimization max-min angle criterion is discussed here. The b a s i c 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 grid. This 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 . method 13 Consider the f o l l o w i n g f i g u r e : Figure 1 Max-min t r i a n g l e o p t i m i z a t i o n I f the smallest angle i s increased by replacing the edge between and T2, the swap i s made c r e a t i n g two new t r i a n g l e s as shown above. Once the t r i a n g u l a r g r i d generation has been b u i l t on the data points, of the functions defined over the t r i a n g l e can proceed. Since, f o r i n t e r p o l a t i o n , the b a s i s functions are u s u a l l y chosen so that they are independent over l o c a l ) no system of equations each t r i a n g l e ( i . e . the b a s i s i s need be solved to obtain the surface. As a r e s u l t , i f the functions are simple ( l i n e a r etc.) the generation of the surface may be postponed u n t i l the value o f s(x,y) i s asked f o r in that particular triangle. The f o l l o w i n g i s an example of the c o n s t r u c t i o n o f a piecewise l i n e a r surface over the t r i a n g u l a r g r i d . The surface can be considered to be constructed o f N piecewise l i n e a r b a s i s functions N (2.7) s(x,y) = ? F i B-i(x y) f j=l where f 1 Bj(Xi, y i ) = . v0 Any i = j otherwise t r i a n g l e i n the g r i d w i l l defined have three non-zero b a s i s on i t . These basis functions are defined t r i a n g l e as f o l l o w s : functions on the canonical 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 ((x,y) i s i n Since any + F3 v with values (Fl,-,F2,-,F3,-) at the v e r t i c e s ) t r i a n g l e i n the domain can be mapped onto the canonical t r i a n g l e using the f o l l o w i n g transformation: (2.10) u = (x - X!) ( y - yi) - (y - y i ) ( x - xj) det(J) 3 v = (y - y ) (xp 1 3 - (x - xj) ( y - y i ) det(J) X]) 2 where /(X2 " *i) \(Y2 - yi) (x - x) 3 x \ (y3 - y i ) / The value of any point i n the domain can be c a l c u l a t e d by: a) Finding the t r i a n g l e which contains i t . b) Mapping the point onto the canonical t r i a n g l e and using (2.3), where F]_, F2, F3 are the values of the f u n c t i o n at the v e r t i c e s of the t r i a n g l e . A number of functions s u i t a b l e f o r surface approximation t r i a n g l e s are given i n ( B a r n h i l l 77). To get over continuity, f i r s t d e r i v a t i v e s must be approximated. Since f o r higher degree CO approximation the c o e f f i c i e n t s are expensive to generate, are u s u a l l y c a l c u l a t e d and stored during the preprocessing The cost per f u n c t i o n evaluation depends on two f a c t o r s : how and they stage. fast 15 the t r i a n g l e which contains the point t o be evaluated can be found and time the complexity o f the basis functions. A reasonable (Lawson 77) i s 0(log2 N ) domain. The piecewise 2 linear f o r an a r b i t r a r y surface point involves search i n the around 15 m u l t i p l i c a t i o n s per f u n c t i o n evaluation plus search time. The cost per f u n c t i o n evaluation i s g e n e r a l l y considerably b e t t e r here than most v a r i a t i o n s of Shepard's formula (for reasonable accuracy). I f not a l l the data points are used t o generate approach can be made t o adapt t o l o c a l domain. In (Fowler and L i t t l e the surface, t h i s "trouble spots" 79) t h i s was done i n the by first generating a t r i a n g u l a r 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 t r i a n g l e d i d not meet i t s accuracy requirements a point was chosen i n t e r i o r t o t h i s t r i a n g l e and the t r i a n g l e was replaced by three new t r i a n g l e s . This process i s repeated u n t i l each t r i a n g l e meets i t s accuracy requirements. When compared t o 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 r e s u l t s were obtained using t h i s method i n terms o f data compaction. I f the user has a package a v a i l a b l e or can a f f o r d t o b u i l d one t h i s i s the method o f choice when working with large numbers o f data points on a r b i t r a r y domains and/or the data points are not r e g u l a r l y spaced. c) Methods based upon squares, (rectangles) Rectangles are the l e a s t f l e x i b l e of the three methods. They have domain restrictions and are not e a s i l y adaptable to i r r e g u l a r l y spaced data. However i f the domain i s s u i t a b l e and the 16 data points are r e g u l a r l y spaced using another method), t h i s (or can be moved to g r i d method has d e f i n i t e lines computational advantages. Most obtained i n t e r p o l a t i o n methods defined over from "transfinite" Coons' patches rectangles can be 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 u n i t square i s given by: (2.11) s(u,v) = [l-u,u] F(0,v) + [F(u,0),F(u l)] f F(l,v) v - [l-u,u] TF(0,0) F(1,0) where 1-v F(0,l)Ul-v F(l,l) (0 < u,v < 1) This patch i n t e r p o l a t e s to the data a l l around the u n i t 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 bilinear other interpolant interpolants results. can be S i m i l a r l y quadratic, c u b i c , and 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° i n t e r p o l a n t s however. Patches which when d i s c r e t i z e d give cl or better surfaces e x i s t . References may be found i n ( B a r n h i l l 77). I f the domain of f(x,y) i s rectangular i n shape and the data points e x i s t on g r i d build a tensor i n t e r s e c t i o n points then i t i s p o s s i b l e to product construction and s o l u t i o n polynomial of a l i n e a r surface. system Although the o f equations i s necessary, the l o c a l support offered by an appropriate basis such as B-splines r e s u l t s i n a system of equations which can be stored 17 and solved obtaining C* c o n t i n u i t y or better i s cheap as i s point evaluation of s(x,y) and storage of efficiently. the surface Using this representation. approach, 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 l e a s t squares f i t to data which i s u s u a l l y d i f f i c u l t to accomplish surface f i t t i n g . / using other approaches to 18 Chapter 3: B-Spline surface approximation. This approach t o surface approximation can be considered as the generalization dimensions of one dimensional spline approximation v i a tensor products. One dimensional to polynomial two splines are q u i t e popular as they are both computationally e f f i c i e n t and q u i t e flexible f o r most applications. They are a l s o well understood mathematically. a) One Dimensional splines: Let f l be a p a r t i t i o n (polynomial) spline a = t 0 < t]_ < ... < t function o f order k (degree M = b, then a k-1) with N - 2 i n t e r n a l knots i s a f u n c t i o n s(x) where s(x) i s a polynomial o f degree < k on each i n t e r v a l ( ( t j ^ t i + j ) i=0 , smooth s p l i n e i f s(x) i s i n C dimensional linear k - 2 [a,b]. M - l ) . s(x) i s s a i d to be a Let P j ir <f denote the f i n i t e space o f piecewise polynomials defined on [a,b] having breakpoints fl and having degree < k. Using the d e f i n i t i o n o f a smooth s p l i n e , i t i s f a i r l y easy t o construct a s p l i n e i n t e r p o l a n t to a given set o f data. Example: Smooth cubic s p l i n e i n t e r p o l a t i o n . 4^ ~t e ti - - - - - — — "t M Figure 2 Smooth s p l i n e s(x) i s a piecewise cubic polynomial which means: (3.1) s(x) = s i ( x ) = a i ( x - t i ) 3 + bi(x-ti) 2 + ci(x-ti) + di 19 since s(x) i s i n C [ a , b ] : 2 (3.2) 1) S i ( t i + 1 ) = s (t i + 1 i + 1 ) = f(t i + 1 3 ) i => a-jhi + b^hj + Cjh£ + 2) s ' i ( t i + 1 ) = s' (t i + 1 i + 1 = dj+ji ) => 3 a j h i + 2bjhj + c j = c j 3) s»i(t ) = s " i+1 (t i + 1 i + 1 ) => 6 a i h j + 2bj = 2bj where h j = ( t j + j - t j ) This implies that we need to solve a system of equations with N+2 unknowns i n N equations t o obtain a l l the c o e f f i c i e n t s . In p r a c t i c e , the two a d d i t i o n a l unknowns are u s u a l l y provided s"(a) = s"(b) = 0 (natural spline) or s e t t i n g s ( a ) 1 by e i t h e r s e t t i n g and s (b) 1 to some approximation of the slope a t the end points of the i n t e r v a l . Note that the above equations form a t r i d i a g o n a l system which i s t o be solved f o r the c o e f f i c i e n t s CJ, i=l,M. This approach, although simple f o r t h i s example i s not e a s i l y applied t o l e a s t squares problems or generalized t o two dimensions. A l t e r n a t i v e l y the problem can be re-stated as (3.3) M s ( X i ) = T a^B^x-i) j=l where each i s a piecewise i = 1,...,M polynomial i n P,g,7rand B^ i s i n C [ a , b ] . 2 Solving the r e s u l t i n g system of equations f o r a j 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 b a s i s 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 b a s i s has, along with the added computation involved in s o l v i n g a f u l l matrix, c o n d i t i o n i n g problems as w e l l . One b a s i s 20 with the desired l o c a l support i s the B-spline b a s i s . (3.4) Bi(x) = [ t i t f i + k ] ( t - x)+' where ( t - x ) = {max (t-x,0) } + and [ t i , . . . ,ti+|<] f denotes t h e k - t h d i v i d e d d i f f e r e n c e o f f u n c t i o n f d e f i n e d r e c u r s i v e l y here as (3.5) [ti,t i + 1 ]f = f(ti) - f(t ti " t i + i + 1 ) 1 [ t i , . . . ,ti ]J f = [tj,...,ti k-i1f - [tj-n,... ,tj k] f t i ~ ti+k + + + I t i s e a s y t o show t h a t B i has l o c a l s u p p o r t , if x > t i + k t h e n ( t - x ) = 0 => B i ( x ) = 0 + i f x < t i t h e n ( t - x ) ^ ! i s a p o l y n o m i a l o f degree k - 1 - [ti, and s i n c e t i + ] f = f k ( f ) where 2, i s some p o i n t on t h e l i n e k B i ( x ) = 0. U s i n g L e i b n i z ' f o r m u l a w h i c h 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) [ti,...,ti + k ]fg = T r=i [ti,...,t ]f r [t ,...,ti r + k ]g we c a n d e f i n e t h e r e c u r r e n c e r e l a t i o n (3.7) [ti,...,ti + k ]f = x - tj ti+k ~ i [ t i , — , t i + k ] f fc + i+k ~ Cti,—/ti+k]g ti+k ~ i fc x fc which c a n be used t o e f f i c i e n t l y g e n e r a t e 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 o v e r an even mesh have some v e r y 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 o f e q u a t i o n s t o s o l v e (de-Boor 79). 21 For (k=4) the corresponding B-spline b a s i s , B j ( x ) , on an even mesh i s given by S((x - Xj)/h) (h i s the distance between mesh p o i n t s ) . x < -2 / 0 (3.8) S(x) = (2 - x ) 24 3 - (1 - x ) 6 3 - x^ + (1 + x ) 4 6 (2 - x ) 24 3 - (1 - x ) 3 - (2 - x ) 3 - (1 - x ) 5 (2 - x ) 3 3 -2 < x <-l -1 < x < 0 0 < x < 1 3 1 < x < 2 ~2T x > 2= 2 For an uneven mesh and a given k (3.7) can be expanded to form the b a s i s functions e x p l i c i t l y . However each basis f u n c t i o n cannot be generated by s c a l i n g a "canonical basis f u n c t i o n " as i s done here f o r an even mesh. For an a r b i t r a r y k, equation (3.7) i s used implicitly to generate the b a s i s functions. There are some problems d e f i n i n g the B-splines at the end of the interval i n question using the current partition. points For the remainder of the t h e s i s 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 d e f i c i e n t . I f two knots coalesce at x j the s p l i n e approximation loses another degree of smoothness at that knot, i f k knots coalesce at X j there i s a jump d i s c o n t i n u i t y at that p o i n t . The knot sequence above applied to N data points s p e c i f i e s a NxN banded system of equations be solved. Solving the system to 22 (3.10) B i (xi) B (xi) Bi(x ) BM ( M) = X M f o r the c o e f f i c i e n t s (from 3.4) 'Fl — M provides the i n t e r p o l a t i n g s p l i n e s(x) This approach can be modified q u i t e r e a d i l y so that s(x) approximates the data instead of i n t e r p o l a t i n g 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 f o r the l e a s t squares approximation, the system (3.11) define (A A)c = A y where c = [ a i , a 2 , . . . , a ^ ] T the T normal equations which can be T solved quite efficiently using a d i r e c t method f o r s o l v i n g banded symmetric p o s i t i v e 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 c o n d i t i o n i n g problems f o r uneven meshes as K(A^A) = K(A)2 B) and K(A) (where K(B) denotes the c o n d i t i o n number of a matrix can get q u i t e large i f an abruptly varying mesh i s used. A knot sequence that r e s u l t s i n a high c o n d i t i o n number can be e a s i l y s p e c i f i e d . 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 w i t h i n the support of that B - s p l i n e ) , then there w i l l be l a r g e and very small diagonal elements i n A A. T However, the zero s t r u c t u r e of the matrix A depends on the placement of the data and s o l v i n g and (for s t o r i n g i t as e f f i c i e n t l y as most problems) makes the quite a t t r a c t i v e . For (A A) B-splines the is difficult. T s o l u t i o n using the practical knots, so This normal equations experience indicates that an uneven mesh w i l l not cause accuracy problems. I t appears that K(A) i s a poor i n d i c a t i o n of what i s happening f o r t h i s system, b) G e n e r a l i z a t i o n to two dimensions G e n e r a l i z a t i o n to two next step in the process. dimensions using A two tensor dimensional basis constructed q u i t e e a s i l y from B-splines using tensor (3.13) ( products i s the function products. V ( x , y ) = Bi(x) Bj(y) i j 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 s p l i n e i n the x and y d i r e c t i o n s as k ky is r e s p e c t i v e l y . A s p l i n e surface approximating the x data could and be characterized using t h i s basis as: (3.14) MxN s(x,y) = 7 gijBiMBjty) if j = l Note that when i n t e r p o l a t i o n i s used, t h i s c o n s t r u c t i o n implies that the data l i e s on a "wire frame" g r i d over a rectangular domain. Obtaining the coefficients for q u i t e straightforward. There are MN s(x,y) using equations i n MN interpolation i s 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: 24 2(k - 1)-1 bands Figure 3_ Banded matrix k and x ky are assumed t o be the same (=k) to s i m p l i f y the d e s c r i p t i o n but because there i s more than one band, (the number of bands depend on k and ky, and the distance between the bands depends x 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 p o s s i b l e , however, to solve t h i s system of eguations by s o l v i n g a s e r i e s of corresponding one dimensional problems i n each d i r e c t i o n . Let L be G M X N the matrix be the matrix of c o e f f i c i e n t s (Gjj = g j j ) and l e t containing the values of the f u n c t i o n 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 a x i s (i.e. Ajj = B i j = B i ( y j ) ) . From (3.14) (3.15) s(x ,yx) k M 2 N 2 9ij i j = kl B i=i j=i N M = 2(2 i=l j=l B L 9 i j Bi(x )) Bj(y ) k N = 2 j=l (3.16) (A< ) . G^ )) BiCy-.) row col R 1 s ( x , y i ) = (A G B ) T k k l Solving A G B^ = L f o r G i s done i n two steps: 2 BJ(XJ) and 25 (3.17) 1) A X(J) = L(J) j = 1,...N => X = G B T =>XT = B GT 2) B G 0 ) = X (J) T j = 1,...,M T where x ( J ) i s the j ' t h column o f X e t c . Since A and B are banded (1 band) matrices o f s i z e M x M and N x N, the storage and computation requirements are reduced considerably over d i r e c t l y s o l v i n g the system. De-Boor i n (de-Boor 79) estimates the complexity of the d i r e c t approach as 0((MN) ) and the separated system 3 as 0 ( M + MN + NM + N ) for a general basis ( f u l l matrices) but the 3 2 2 3 d i f f e r e n c e (for small k x and ky) i s even l a r g e r using a l o c a l b a s i s , 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 p o i n t s than approximating f u n c t i o n s . I f the d i r e c t approach i s used, an MN x MN system o f equations again has t o be solved t o o b t a i n the c o e f f i c i e n t s when the normal equations are used. We no longer have the strict c o n s t r a i n t s on the placement o f the data, c o n s t r a i n t that there must be a data point i n the support basis function (where the f u n c t i o n i s non-zero). only the o f each I t i s again advantagous t o form the normal equations, A A x = A b , i n s p i t e o f T possible c o n d i t i o n i n g problems as the zero T structure f o r A A i s T dependent only on M, N, and the support o f the b a s i s functions which makes the matrix e a s i e r t o store and manipulate. The matrix i s banded with bands separated by e i t h e r a distance o f M or N depending on the ordering. I f the data p o i n t s are on a rectangular g r i d , the problem can again be reduced t o s o l v i n g a s e r i e s o f one dimensional problems as i n 26 equation (3.17). (3.18) A?A G B B = ATL B T w h i c h c a n be s o l v e d (as f o r i n t e r p o l a t i o n ) i n 2 s t e p s : 1) (A A) xO) 2) (BT ) G (3) T B which T = (A L)(J) T j = l,...L (BX )(J) T = j =1,...,M i s e q u i v a l e n t t o (3.16) g i v e n and B B a r e banded m a t r i c e s T used to advantage in n t h e same d a t a . From (3.16) A A T (1 band each) a p r o p e r t y w h i c h c a n a l s o be constructing A A T efficiently. Although t h e o r e t i c a l l y t h e r e c o u l d be problems because o f t h e l a r g e c o n d i t i o n number t h a t t h i s system o f e q u a t i o n s c o u l d 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 . W e i g h t s c a n be added t o t h i s p r o c e s s a s i n t h e l a r g e r system A W A G B W B (3.19) where W T X x y T = A W L B W T X y B T and Wy a r e d i a g o n a l m a t r i c e s 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. T h i s i s c o n s i d e r a b l y d i f f e r e n t t h a n = A Wb i n t h e f u l l T point. system where each w e i g h t a p p l i e s t o a s i n g l e A WAX T data 27 Chapter 4: Adapting the Spline Surface to Quickly Varying Data and Non-Rectangular Domains Least squares method described evaluate. spline However, they and generated using i n Chapter 3 are quite inexpensive have some l i m i t a t i o n s adapt to "trouble" areas abruptly) surfaces on the separation to construct and their ability (where the surface to be approximated changes t h e i r a b i l i t y to adapt to non-rectangular domains. To change the shape of the spline surface either weights may the knot spacing may to be added or be changed. We discuss the knot placement f i r s t , a) Knot Placement Changing the knot sequence i n each d i r e c t i o n can dramatically however, improve since only the surface grid lines approximation. can be In i n some cases other moved, changing the cases, knot sequence may be i n e f f e c t i v e . I f we have a surface to be approximated which contains a single l o c a l trouble spot, knots can be added u n t i l the desired accuracy i s reached i n 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 i s not l o c a l to one area and not p a r a l l e l to one of efficiently. the axes, i t may be impossible to adapt the mesh 28 X Figure 5_: Global trouble area Only increasing the number of knots everywhere will improve the accuracy i n t h i s example. This time the whole surface gets a great deal more a t t e n t i o n , i n c l u d i n g 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 c a r e f u l placement o f knots. Since only knot " l i n e s " can be moved, i t i s p o s s i b l e to use one dimensional knot o p t i m i z a t i o n methods along these l i n e s f o r c e r t a i n problems. There are two b a s i c 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 r e s u l t s 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 f i n d i n g a good knot placement i s o u t l i n e d below, f o r more d e t a i l s see (Dodson 72). The distance from f ( x ) to s(x) i n the L2 norm i s given by: b (4.1) I If - s| | = ( J l f - s|P d x ) V p p a 29 Dodson presents the f o l l o w i n g approximation to t h i s d i s t a n c e : b (4.2) |! f - slip < C $ | f ( t ) f dt k a where cv = ^[k + l ] ^ 1 E q u a l i z i n g the distance from f(x) to s(x) across each s u b i n t e r v a l so that ti+1 (4.3) b $|f | k ff dt = 1 n - k ti $|f k f dt a should provide a good knot placement. Dodson suggests the f o l l o w i n g as an estimate f o r f ( t ) as i t i s k unavailable f o r most a p p l i c a t i o n s : a) Generate a s p l i n e approximation preliminary knot placement with s(x) to the data with some the d e s i r e d ' number of knots (usually provided by the user). b) D i f f e r e n t i a t e s(x) k-1 times t o obtain the piecewise constant function s l ( x ) . k _ c) J o i n the midpoints o f each adjacent piecewise s ~ ^ ( x ) to form a new piecewise l i n e a r f u n c t i o n f k k - 1 constant i n (x). d) D i f f e r e n t i a t e f l ( x ) to o b t a i n the piecewise constant f u n c t i o n k - f (x). k This approximation to f ( x ) i s used i n (4.3) to generate the d e s i r e d k mesh. I f the o r i g i n a l s p l i n e f i t s(x) approximates the f u n c t i o n f ( x ) reasonably w e l l then using the new knot placement obtained above w i l l 30 g e n e r a l l y improve the f i t of the s p l i n e function. However i f s(x) does not f i t the data w e l l to begin w i t h , the new knot sequence may not be an improvement over the o l d and the e r r o r I|f - s l i p may grow. Applying this any other dimensional one i f not impossible to do e f f e c t i v e l y . I f the problem has a then automatic selection difficult (Figure 4) i n an knot to t r o u b l e spot splines dimensional technique local two or manner i f the scan l i n e s f o r the dimensional o p t i m i z a t i o n are chosen c o r r e c t l y , the knots suggested the o p t i m i z a t i o n technique (along these lines) will improve is one by the o v e r a l l surface approximation. I f there i s more than one t r o u b l e spot i t i s not so obvious how to apply the o p t i m i z a t i o n 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 surface to choose good knot fitting problem with sequences f o r a two spline products. It dimensional 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 p o s i t i o n of the knots and deciding when the f i t i s acceptable (or when i t s time to g i v e up). In the package described i n Appendix A, a user can opt to apply the one dimensional o p t i m i z a t i o n techniques sequences to be used f o r h i s problem. provide a tool to suggest i n t e r a c t i v e l y by the user or knot to suggest A l l that i s intended sequences (more l i k e l y ) most powerful t o o l a v a i l a b l e to generate which can be knot i s to altered completely ignored. The a good s p l i n e 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 s p l i n e and to see the r e s u l t s . 31 b) Weight Placement Weights can a l s o be applied to the surface to a l t e r the shape of the s p l i n e approximation. They can be used i n two ways, e i t h e r to force a c l o s e r f i t i n one area by increasing the weight i n that area or to ignore the effect of some "impossible" area accuracy there to increase accuracy on the (sacrificing rest of the surface). Weights must be applied c o n s e r v a t i v e l y , however. Improperly they can cause severe conditioning problems u n l i k e the large c o n d i t i o n number caused i n the system applied which, by an i r r e g u l a r mesh, may have a severe e f f e c t on the accuracy of the computed s ( x , y ) . As stated i n Chapter 3, weights cannot be applied to i n d i v i d u a l data points i f equation (3.19) i s used to o b t a i n 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 s i 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 t r o u b l e spot without applying i t along the g r i d l i n e s outside the t r o u b l e area (a nonlocal weighting r e s u l t s ) . I f the problem i s not l o c a l Figure not 5, then the addition of weights will help as i n the f i t . Providing a system that a u t o m a t i c a l l y 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, d e l e t e and change weights, however, i s a very e f f e c t i v e approach to the problem. c) Relaxing Domain R e s t r i c t i o n s 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 efficiently 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 i s 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 method before complex, our approach may be taken. could also add considerable surface fitting The transformation, i f point evaluation and preprocessing cost. If instead, a rectangle D' i s superimposed over D, i t i s possible to generate a surface over D' providing that data i s given i n D' D and the data l i e s along mesh lines. The cost per point evaluation i s increased however as checking whether or not a point (x,y) i s i n D must be done for each point evaluation. D' = D U R Figure (5: Non-rectangular domain The problem now i s how to generate data points i n R which w i l l not adversely affect the spline f i t on D. I f i t were possible to weigh each data point individually, then R could be ignored but because of the way s(x,y) i s 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 many d i f f e r e n t methods of attempting t h i s goal we have chosen the rectangle be the following a) For each region R in D find smallest which encloses R. b) Use a v a r i a n t of a b i l i n e a r Coons' patch over t h i s rectangle to generate data i n s i d e 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 iteration does not improve the f i t or the f i t i s s a t i s f a c t o r y . I f the region R i s not rectangular, a weighted surface to reduce the e f f e c t of generated data on D may approximation be undesirable since i t w i l l a l s o a f f e c t the f i t along the boundary between D and R i n s i d e D. The rectangles enclosing R and the d i s c r e t i z e d Coons' patch are obtained i n the f o l l o w i n g 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 p o s i t i o n i s recorded i n two weight vectors at WXj and WYj. i i ) For each region the f o l l o w i n g 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 s t r a i g h t l i n e i n t e r p o l a t i o n between the two end points on that side i n D j u s t outside R. b) I f two adjacent sides are on the boundary of the rectangle, 34 t h e v e r t e x ( c o r n e r o f D') of v a l u e i s f i r s t computed as t h e average t h e two n e a r e s t c o r n e r s o f t h e r e c t a n g l e around R and (a) i s a p p l i e d t o b o t h edges. c) The Coons' p a t c h i s t h e n d i s c r e t i z e d by l e t t i n g t h e f u n c t i o n v a l u e s on t h e edge o f t h e u n i t square be t h e d a t a v a l u e s on rectangle e n c l o s i n g R. The discretized Coons' patch is the then e v a l u a t e d a t each p o i n t i n s i d e R. Generating data this way has varied results. Since for a non-rectangular R there are r e a l f u n c t i o n v a l u e s 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 t h e r e may be a jump d i s c o n t i n u i t y a t t h e boundary between R and D i f t h e r e a l d a t a i n t h e a r e a does n o t agree w i t h t h e Coons' p a t c h a p p r o x i m a t i o n . I f t h e s u r f a c e around R i s smooth t h e s p l i n e f i t may n o t be a d v e r s l y a f f e c t e d s i n c e t h e d i s c o n t i n u i t y s h o u l d be However i f t h e s u r f a c e moves q u i c k l y around R t h e n t h e w i l l p r o b a b l y be l a r g e and t h e o s c i l l a t i o n s caused to damp o u t . small. discontinuity by i t w i l l be hard 35 Chapter 5_ : Results and Evaluation To illustrate the e f f e c t i v e n e s s o f the i n t e r a c t i v e spline approximation, four t e s t functions were chosen. The f i r s t three were taken from Franke who used these functions t o t e s t a l a r g e number o f methods f o r i n t e r p o l a t i n g scattered data. Some accuracy and timing r e s u l t s from h i s paper have been included. However, since the methods discussed i n (Franke 79) are f o r scattered data, d i f f e r e n t data points had t o be used f o r the generation consequently, comparisons o f our s p l i n e surface and should be made c a u t i o u s l y as the methods have very d i f f e r e n t areas o f a p p l i c a t i o n . The fourth t e s t f u n c t i o n i s included t o demonstrate the knot s e l e c t i o n algorithm and the data i n s e r t i o n technique. The domain o f the three t e s t f u n c t i o n s from (Franke 79) i s the u n i t square. To improve the s c a l i n g o f the graphics, the sampled data was t r a n s l a t e d and scaled along the x and y axes to the square [-1,1] X [-1,1]. The three t e s t functions from Franke are: Saddle function: f l ( X , Y ) = (1.25 + cos(5.4Y))/(6(1 + ( 3 x - l ) ) 2 Exponential f u n c t i o n : f2(X,Y) = .75 EXP ( - ( 9 X - 2 ) + (9y-2) )/4 ) 2 2 + .75 EXP ( -(9X+l) /49 - (9Y+1) /10 ) 2 + .5 EXP ( -((9X+7) 2 2 + (9Y-3) )/4 ) 2 - .2 EXP ( - ( 9 X - 4 ) - (9Y-7) ) 2 C l i f f function: 2 f3(X,Y) = 1/9 (TANH(9Y - 9X) + 1 ) The t e s t functions are p l o t t e d below. In Figures 7, 9 and 11 the X a x i s i s h o r i z o n t a l , p a r a l l e l t o the page and the y a x i s i s pointed away from the reader. In Figures 8 and 10 the p i c t u r e has been 36 F i g u r e 7: F u n c t i o n f l (Saddle function) .37U88lJ8yUF.T0O .83327S297E-02 Figure 8: Function fl rotated 90 degrees The second test function contains two h i l l s and a depression: SURFACE FITTING . lOUt97S93Ef01 - . 1 8 5 0 0 4 7 1 IE + 00 .939S97893E+D; ( 1 1 1 ! i i i 1 1 i Figure 9 Function f2 (Exponential test function) SURFACE F I T T I N G . 1 0 U 1 9 7 5 9 8 E + 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 domain: F i g u r e 11 F u n c t i o n f 3 ( C l i f f t e s t f u n c t i o n ) the 39 The accuracy o f the s p l i n e f i t f o r these examples i s measured using 625 p o i n t s on the domain o f the function. The f i v e worst e r r o r s ( | f ( x i , y i ) - s ) ( x i , y j ) | ) and t h e i r p o s i t i o n s are given as w e l l as an approximation t o the l e a s t squares e r r o r . The s p l i n e surface f o r each example i s generated using 2500 data p o i n t s evenly spaced domain (unless otherwise s p e c i f i e d ) . Table I : E r r o r s on function f l using evenly spaced knots KNOTS ABS 10 X 10 1.08E-3 at (-.265,-1 ) 1.06E-3 (-.265,-.837) 1.05E-3 (-.265,-.918) .935E-3 (-.265,-.755) •852E-3 (-.265,-.673) 2.60E-4 30 X 30 1.43E-5 1.19E-5 1.19E-5 .894E-5 .835E-5 3.46E-7 ERROR L.S. ERROR (-.265,-1 ) (-.265,-.918) (-.265,-.837) (-.347,-1 ) (-.265,.959 ) 40 The second t e s t function i s more d i f f i c u l t t o f i t accurately and was used by Franke as h i s primary t e s t f u n c t i o n . This function i s examined more e x t e n s i v e l y here as seen i n the f o l l o w i n g r e s u l t s : The number o f knots, the knot spacing and the number o f data points are a l l varied. Table I I : E r r o r s on function f2 using evenly spaced knots KNOTS ABS 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 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 ERROR L.S. ERROR (225 data points) 10 X 10 (even) (100 data points =0> i n t e r p o l a t i o n ) 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 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 (2500 data points) 30 X 30 (even) 41 With a s e t o f good knots, which were chosen using the package, the improved r e s u l t given i n Table I I I are obtained. 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 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 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 (225 data points) 10 X 10 Table IV l i s t s some e r r o r r e s u l t s from function f 4 . Notice that the optimal knot spacing over t h i s f u n c t i o n i s an even mesh and the optimal weighting i s a l s o an even d i s t r i b u t i o n . I t i s not p o s s i b l e for the s p l i n e to accurately f i t this function without spending a d d i t i o n a l time on areas that do not need the a t t e n t i o n . Table IV: E r r o r s on f u n c t i o n f 3 KNOTS ABS 10 X 10 (even) 5.07E-3 5.07E-3 4.84E-3 4.84E-3 4.72E-3 ( .224,.612 ) ( .612,.224 ) ( .224,-.184) (-.184,.224 ) ( .469,.612 ) 1.84E-3 30 X 30 (even) 1.95E-5 1.91E-5 1.88E-5 1.87E-5 1.85E-5 ( .142,.224 ) (-.673,-.591) ( .633,.551 ) ( .224,.142 ) (-.183,-.265) 3.99E-6 Table V lists ERROR L.S. ERROR some representative results from (Franke 79). The 42 r e s u l t s presented were obtained by t e s t i n g the method on 100 scattered data p o i n t s . Table V: E r r o r s on function f2 by selected methods from (Franke 79) METHOD f l ABS. f2 f3 L.S. ABS. L.S. ABS. L.S. .00194 .0573 .0128 .0468 .00551 .273 .0620 Methods based upon p o i n t s : Mod. Quad. Shepard .0125 Shepard Methods based upon t r i a n g l e s : 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 (2.5) modified quadratic Shepard's formula to have local support. i s basically Both equation Lawson's method and Akima's method use b a s i s functions over a t r i a n g u l a r g r i d . Franke's method uses l o c a l Hardy's method discussed. the (Thin plate) functions over a rectangular i s somewhat d i f f e r e n t from the methods mesh. previously I t uses a l o c a l basis f u n c t i o n ( l i k e the B-splines) but support o f the b a s i s sparse system of equations function i s a c i r c u l a r region. A large ( l i k e the one which could r e s u l t from using B-splines f o r scattered data) must be solved to generate the surface. 43 In Table VI we present and some timing r e s u l t s from the generation evaluation of our s p l i n e surface. Table V I I presents r e s u l t s given i n (Franke 79) f o r the methods presented the timing i n Table V. The times are i n seconds. An Amdahl 470-V6-II was used to generate and evaluate the s p l i n e surfaces. An IBM 360-67 was used for a l l timings taken from (Franke 79). Table VI: Timings from generation and evaluation o f s p l i n e SPLINE SURFACE PREPROCESSING EVALUATION (Generating the s p l i n e ) (625 points) (100 basis functions) 0.13 0.26 (900 basis functions) 0.19 0.26 0.05 0.26 (2500 data points) (100 data points) (100 basis functions) Evaluation speed could be increased considerably since the algorithm used to generate the B-splines t o evaluate the s p l i n e surface i s f o r s p l i n e s of a r b i t r a r y order and therefore quite slow. Table V I I : Timings on selected methods from (Franke 79) METHOD PREPROCESSING EVALUATION (1089 data points) Mod. Quad. Shepard 8.6 Shepard 15 17 Lawson 1.8 1.7 Akima 2.2 0.8 Franke-TPS 2.7 6.5 Hardy 7.1 13 44 Even allowing f o r the d i f f e r e n c e i n speed between computers, the e f f i c i e n c y o f the one dimensional s p l i n e approach t o surface fitting i s c l e a r l y demonstrated. The fourth t e s t function was included t o show the a d a p t a b i l i t y o f the s p l i n e using a good knot placement, and the a b i l i t y o f the s p l i n e to handle non-rectangular domains. The domain o f the f i r s t example o f function f4 i s the square ([-1,1] X [-1,1]). The function i s given by: f4(X,Y) = 1/(1+25X ) COS (11/2 Y) 2 This function i s one dimensional i n nature. The surface most along the X a x i s near the o r i g i n . SURFACE FITTING . 10O2'?S612ET01 .12U76S349E-06 . 5 8 H 0 S 3 . 7 1 5 E + 01 I i ' !— Figure 12 P l o t o f function f4 v a r i e s the 45 S U R F A C E f ITT J NG . 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 e r r o r r e s u l t s f o r a s p l i n e o f even mesh over function f4 are given i n Table V I I I . 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 KNOTS ABS 13 X 13 (even) 1.33E-2 a t ( .143,-.204) 1.32E-2 ( .143,.612 ) 1.31E-2 ( .143,-.102) 1.29E-2 ( .143,.142 ) 1.28E-2 (-.204,-.204) A good knot sequence ERROR L.S. ERROR 3.71E-3 was chosen using routines based upon the one dimensional o p t i m i z a t i o n technique presented by Dodson as described i n Chapter 4. 46 Table IX: S p l i n e approx. of f u n c t i o n f4 over an optimized mesh KNOTS 13 X 13 The ABS (opt) ERROR 2.73E-3 2.71E-3 2.68E-3 2.65E-3 2.62E-3 o p t i m i z a t i o n routines L.S. ERROR ( ( ( ( ( .224,-.204) .224,.612 ) .224,-.102) .224,.143 ) .224,-.184) w i l l not 7.74E-4 produce a b e t t e r knot unless the f i t i s already q u i t e good. In t h i s case, 13 knots were needed along the X a x i s before the routines would work. Table X presents the r e s u l t s of applying the techniques described in Chapter 4 to extend our s p l i n e method to non-rectangular domains. The domain of the f u n c t i o n to be approximated i s defined by the supplied Fortran f u n c t i o n "RI" domain of f ( x , y ) , otherwise (Rl(x,y) i s true i f (x,y) user i s i n the i t i s f a l s e ) . For t h i s example the domain of the f u n c t i o n i s L-shaped with ([-1,0] X [-1,0]) as the "hole" i n the rectangular domain. The surface surface would varies very rapidly show much b e t t e r near results (0,0). than this A more gentle example. Error estimates are given a f t e r the i n i t i a l patch and a f t e r each successive surface "hole" i s generated. The as l i n e s that maximum e r r o r s tend the weights are "touch" the to move out increased. Weights are hole. As a result applied the weight on of the to scan the data points i n the hole i s the product of the weight on the two l i n e s that i n t e r s e c t at that data p o i n t . T a b l e X: S p l i n e approx. o f f u n c t i o n f 4 over an L shaped doma i KNOTS ABS 13 X 13 (OPT) (PATCH WEIGHT 0.2) 1.16E-2 AT 1.14E-2 1.11E-2 1.06E-2 1.01E-2 13 X 13 (opt) (SMOOTH WEIGHT .2) ERROR L.S. ERROR -.204,-.428) -.204,-.346) -.204,-.510) -.204,-.265) -.204,-.591) 1.67E-3 5.53E-3 3.17E-3 3.16E-3 3.06E-3 3.01E-3 .959,-1 ) .224,-.265) .224,-.183) .224,-.346) .224,-.102) 8.06E-4 13 X 13 (opt) (SMOOTH WEIGHT .5) 5.53E-3 3.04E-3 2.91E-3 2.86E-3 2.86E-3 .959,-1 ) -.347, .142) .224,-.183) .224,-.102) .224,-.265) 7.69E-4 13 X 13 (opt) (SMOOTH WEIGHT 1 ) 5.53E-3 3.08E-3 2.85E-3 2.82E-3 2.78E-3 .959,-1 ) -.347, .142) .224,-.183) .224,-.102) .224,-.265) 7.64E-4 48 Chapter 6 : Conclusions A method of generating a polynomial approximating a. f u n c t i o n f(x,y) represented spaced points data over a "nearly" spline surface by a set of s(x,y) regularly rectangular g r i d has been presented. The surface, although " g l o b a l " i n nature, can be generated efficiently for fairly large sets of evaluation cost i s low and can be reduced data points. The point i f the exact degree of the s p l i n e 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 c o s t may be reduced still 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 q u i t e hard f o r the human eye to detect d i s c o n t i n u i t y i n d e r i v a t i v e s higher than 2 i n a f u n c t i o n , the surface approximant needs to be no more than C 2 to appear polynomials") smooth. For lower the appears surface degree fits smooth (say and k does = 4 not "cubic tend to o s c i l l a t e as much as surfaces based upon higher degree polynomials. When comparing this method to other approaches of surface f i t t i n g , i t should be noted that other methods g e n e r a l l y do not have all the restrictions approximate the this data i n t e r p o l a t e i t . We method (least have methods f o r i n t e r p o l a t i n g has squares (perhaps but or unfairly) very few are otherwise) able rather to than compared the method to scattered data discussed i n (Franke 79). The b e t t e r methods presented there f i t the t e s t f u n c t i o n s f l , f 2 , and f3 i n Chapter does almost presented 5 w e l l using one hundred data p o i n t s . The as w e l l ( i n terms of accuracy) as spline f i t the b e t t e r methods on t h i s amount of data with an even mesh and 100 basis functions ( i n t e r p o l a t i o n ) . Since a least-squares f i t i s used, i t is 49 p o s s i b l e t o increase the number of data points without increasing the number of b a s i s f u n c t i o n s . As the number of data points i s increased the accuracy comparable o f the s p l i n e with computational the cost best is still f i t i s also results very in increased (Franke reasonable. 79) and becomes while The surface the fitting package was used to c a r e f u l l y choose a good knot sequence over the f u n c t i o n f 2 . The accuracy o f 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 r e s u l t s i n Chapter 5. The cost of generating the s p l i n e surface and evaluation (even with a r e l a t i v e l y l a r g e number of points) i s nominal compared to the cost o f 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 (xi,yi,Zj) triplets. Similarly, the B-spline representation of the surface approximant i s very compact; there are two vectors s p e c i f y i n g the l o c a t i o n of the surface patches (knots f o r 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 s p e c i f y the s p l i n e surface over these patches. In terms o f space and speed t h i s makes the s p l i n e 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 c o n s t r a i n t s on a problem or the data has e r r o r s 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 l e a s t squares s p l i n e f i t over the re-sampled data. The interactive environment supplied by the surface fitting 50 package provides a powerful t o o l to the user: A user can dynamically a l t e r the f i t o f the s p l i n e surface u n t i l he i s s a t i s f i e d with the results (the graphics facilities 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 t o improve the behavior of polynomial s p l i n e product surfaces. There may be some promise i n surfaces generated using r a t i o n a l s p l i n e s or non-polynomial functions but we b e l i e v e the current research i n t o t h i s approach i s w e l l developed. There i s considerable opportunity f o r research more general two dimensional surface approximation In conclusion the graphics e f f e c t i v e l y solves two dimensional package into techniques. provided in this thesis surface f i t t i n g problems provided that the c o n s t r a i n t s on the data are met. 51 BIBLIOGRAPHY B a r n h i l l , Robert E. and Riesenfeld R.E., Computer Aided Design. New York: Academic Press, 1974. Geometric 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 I I - V a r i a b l e Knots". Purdue U n i v e r s i t y : Technical Report CSD TR 21, 1968. de-Boor, C a r l , "Package f o r C a l c u l a t i n g Numer. Anal.. 1977, pp 441-473. de-Boor, C a r l , V e r l a g , 1979. A Practical Guide with to S p l i n e s . B-splines".SIAM J . New York: Springer de-Boor, C a r l and Rice, John R., "An Adaptive Algorithm f o r Multivariate Approximation Giving Optimal Convergence Rates". U n i v e r s i t y of Wisconsin: MRC-TSR-1773, 1977. Coons, S. A., "Surfaces f o r 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 S p l i n e Functions". Purdue U n i v e r s i t y : Ph.D. t h e s i s , 1972. Fowler, R. J . and L i t t l e J . L., "Automatic E x t r a c t i o n of I r r e g u l a r Network D i g i t a l T e r r a i n Models", ACM SIGGRAPH, 1979, pp 199-207. Franke, R., "A C r i t i c a l Comparison of Some Methods f o r I n t e r p o l a t i o n of Scattered Data". Naval Postgraduate School: NPS-53-79-003, 1979. Jupp, D. L. B., "Approximation to Data Numer. Anal.. 1978, pp 328-343. With Free Knots". SIAM J . Lawson, C. L., "Software f o r C Surface I n t e r p o l a t i o n " . Mathematical Software I I I , ed. Rice, J . R., London: Academic press, 1977, pp 161-193. 1 McLain, D. H. "Drawing countours from Computer J o u r n a l . 1974, pp 318-324. a r b i t r a r y data p o i n t s " . The M i 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 I n t e r p o l a t i o n Function f o r I r 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., S p l i n e Algorithms f o r Curves U t i l i t a s Mathematica i n c . , 1974. and Surfaces. Winnipeg: 52 S t r a n g , G. and G. F i x , An A n a l y s i s New J e r s e y : P r e n t i c e - H a l l , 1973. of the Finite Element Method. 53 Appendix A : A Surface F i t t i n g Package The code f o r the package i s divided into two s e c t i o n s : The command i n t e r p r e t e r and the support routines. a) The Command I n t e r p r e t e r The command i n t e r p r e t e r i s designed t o provide environment which a user can f i t a tensor product an i n t e r a c t i v e polynomial s p l i n e surface t o h i s data. I t s command s t r u c t u r e i s designed to protect the user as much as i s p r a c t i c a l while t r y i n g not t o "get i n the way". A d e s c r i p t i o n o f the basic data s t r u c t u r e s and commands f o l l o w : 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 a r e : 1) Data: The data that s p e c i f i e s the surface c o n s i s t s of L - A matrix of dimension a t l e a s t LM x LN which contain the values of F^ a t the g r i d i n t e r s e c t i o n p o i n t s . XX - A vector of dimension a t lease LM containing the x coordinates. XY - The corresponding vector i n the y d i r e c t i o n . WX - A vector s p e c i f y i n g the weights t o apply t o the g r i d l i n e s s p e c i f i e d i n XX. WY - The corresponding vector i n the y d i r e c t i o n . LM - the number of g r i d l i n e s perpendicular t o the x a x i s . LN - the number of g r i d l i n e s perpendicular t o the y a x i s . LMDIM - the f i r s t dimension of the matrix L. 54 2) Spline Surface: The components that s p e c i f y the s p l i n e surface approximant. GAMMA - A matrix of dimension a t l e a s t M x N containing the c o e f f i c i e n t s o f the B-spline surface approximant. TX - The knot sequence s p e c i f y i n g the B-splines i n the x direction. TY - The corresponding vector i n the y d i r e c t i o n . KX - The order o f the B-splines i n the x d i r e c t i o n . KY - " M - The number o f B-splines i n the x d i r e c t i o n . N - " MDIM - The f i r s t dimension o f GAMMA. " y direction. " y direction. 3) Graphics: The i n t e r n a l 3-D g r a p h i c a l object (IMAGE) X - A vector s p e c i f y i n g the p o s i t i o n o f the g r i d l i n e s perpendicular t o the x a x i s . Y - The corresponding y vector. Z - The value o f the surface a t the points X 0) Y NX - The number o f g r i d l i n e s i n x. NY - The number o f g r i d l i n e s i n y. NZ - The f i r s t dimension o f the matrix Z. 4) Graphics: The IG representation o f the surface (DISPLAY) (data a v a i l a b l e i n surface package) THETA - The angle the surface i s viewed from PHI (Azimuth). - The height (Elevation) surface i s viewed from. 55 The system initializing, keep provides commands reading and w r i t i n g the system consistent for generating, these objects. (e.g. a user may I t also changing, t r i e s to not generate the c o e f f i c i e n t s o f a s p l i n e approximant before the data has been read i n and the knots have been s p e c i f i e d ) . Any non-ambiguous subset o f a command may be used t o s p e c i f y that command cases l i n e (e.g INIT), and i n most boundaries are ignored when reading the commands. The commands a v a i l a b l e are: INITIALIZE I n i t i a l i z e mode allows the user t o i n i t i a l i z e the f o l l o w i n g variables: ORDER KX KY I n i t i a l i z e s the order o f the s p l i n e s i n the x and y d i r e c t i o n . I t i n v a l i d a t e s the current s p l i n e 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 t o be used i n the system. The knots are i n i t i a l i z e d t o an even mesh with M B-splines i n the x d i r e c t i o n and N i n the y. The knots are re-initialized 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 o f (evenly spaced) mesh l i n e s i n the d i s p l a y . The d i s p l a y l i n e s are a l s o i n i t i a l i z e d every time the bounds may change. 56 REGION i Indicate which user supplied f u n c t i o n i s t o be used t o s p e c i f y the domain o f the f u n c t i o n t o be approximated. I f the domain i s rectangular no f u n c t i o n needs t o be s p e c i f i e d I f a region has been p r e v i o u s l y s p e c i f i e d then invoking the command "REGION" with no region s p e c i f i e d cancels the region (the knots s p e c i f y a rectangular region i m p l i c i t l y ) . END To return t o the top l e v e l . CHANGE This mode 4 i n t e r n a l sections which allow the user t o change the KNOTS, DATA, WEIGHTS and DISPLAY data s t r u c t u r e s i n t e r a c t i v e l y . ADD K Add K knots t o the current knot sequence. The user chooses a knot by p o s i t i o n i n g the cross h a i r s (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 t o obtain a good knot placement (for a one dimensional problem along a scan l i n e ) using the knot placement algorithm as discussed i n Chapter 4. 57 PICK I n t e r a c t i v e l y choose the scan l i n e the o p t i m i z a t i o n technique i s t o be applied t o . TRY I SAVE K Apply the o p t i m i z a t i o n technique t o the current scan l i n e . Save the suggested mesh i f SAVE i s s p e c i f i e d . ROTATE Rotate the viewing d i r e c t i o n 90 degrees (change the a x i s p a r a l l e l to the screen). END Return t o KNOT mode. ROTATE Rotate the viewing a x i s . END Return t o CHANGE mode. DATA This mode f a c i l i t a t e s the generation o f data as described i n Chapter 4. These commands do not invoke any graphics activities. PATCH Use a Coons' patch t o patch the holes i n the domain. A region must already have been s p e c i f i e d i f t h i s command i s invoked. 58 SMOOTH fWEIGHT WX Generate a s p l i n e 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 s p l i n e . END Return t o CHANGE mode. WEIGHT Weight mode allows the user t o i n t e r a c t i v e l y change the weights applied t o the ( l i n e s o f data) over the surface. WEIGHT WT S p e c i f y the value o f the weight t o be i n s e r t e d . CHANGE K I n t e r a c t i v e l y change the weight t o WT on K data l i n e s . ROTATE Rotate t o other a x i s . END Return t o CHANGE mode. DISPLAY I n t e r a c t i v e l y change the g r i d spacing on the i n t e r n a l g r a p h i c a l object (IMAGE). I t o n l y allows the user t o change the g r i d i f a s p l i n e surface i s t o be displayed as the data i s d i s c r e t e and cannot be sampled a r b i t r a r i l y . ADD K I n t e r a c t i v e l y add K l i n e s on the d i s p l a y . This does not immediately a f f e c t the g r a p h i c a l o b j e c t , screen or image. DELETE K I n t e r a c t i v e l y d e l e t e K l i n e s t o the d i s p l a y . REFRESH Refresh the screen adding the new d i s p l a y l i n e s . ROTATE Rotate the d i s p l a y on the screen adding the new d i s p l a y lines. END Return t o CHANGE l e v e l . END Return t o top l e v e l . REFRESH/DISPLAY Draw screen using current (IG) g r a p h i c a l o b j e c t . ROTATE Rotate the d i s p l a y on the screen. ( COEFFICIENTS IMAGE K READ \ KNOTS ( FROM filename DATA Read the corresponding item from the f i l e filename. The formats f o r the f i l e s are contained i n the routines RDATA e t c . the I/O s e c t i o n of the subroutine l i b r a r y . For "IMAGE" K must be i n (0 .. 3) where 0 means the current d i s p l a y (SCREEN) which i s used t o generate the IG g r a p h i c a l o b j e c t . 60 j COEFFICIENTS IMAGE K WRITE J KNOTS ONTO f i l e n a m e DATA SCREEN / DEVICE= 'CALC' 'PRINT V Write the corresponding 1 item onto the f i l e filename. COEFFICIENTS ERROR GENERATE DISPLAY IMAGE FROM f SPLINE ^ SURFACE V ERROR Generate t h e c o r r e s p o n d i n g item. MOVE(IMAGE KITO IMAGE L SCREEN SCREEN Move t h e c o r r e s p o n d i n g item SUBTRACT|IMAGE KjFROM ( IMAGE L 1 SCREEN (K i n (1 .. 3 ) ) . 1 j SCREEN S u b t r a c t t h e c o r r e s p o n d i n g u n i t s p l a c i n g the r e s u l t i n t h e second e n t r y . END Terminate t h e r u n . 61 b) The Support Routines The support routines are d i v i d e d into f i v e main s e c t i o n s : GRAPHICS Routines that generate the g r a p h i c a l o b j e c t , d i s p l a y i t on the screen, read p o s i t i o n s f o r d i s p l a y l i n e s and knots o f f of the screen e t c . I/O Routines that handle the input and output o f DATA, KNOTS, COEFFICIENTS, & UNITS (Data that s p e c i f i e s enough information to generate a g r a p h i c a l o b j e c t ) . These routines should be used to prepare the data and knots when using t h i s package to generate a s p l i n e surface. OPTIMIZATION Routines that handle the one dimensional o p t i m i z a t i o n . PATCHING Routines f i n d the rectangles on the rectangle D 1 that enclose the regions not i n the domain o f the function to be approximated and generate the data t o f i l l these "holes" using a Coons' patch. GENERATION AND EVALUATION The routines that generate the s p l i n e surface, evaluate i t and provide some e r r o r estimates. 62 Appendix B : Generation and Evaluation These two routines use the s p l i n e package i n (de-Boor 77) t o generate and evaluate a tensor product polynomial s p l i n e surface o f a r b i t r a r y order. C C C C C C C C "SPLINE GENERATION AND EVALUATION" 1 CALCULATE THE COEFFICIENTS "GAMMA" FOR THE B-SLINE SURFACE APPROXIMATION GAMMA TX TY KX KY M N NDIM XX XY WX WY L LM LN LMDIM A X c c c c c c c c c c c c c c c c c c c c C C C 10 SUBROUTINE LSCOEF(GAMMA, TX, TY, KX, KY, M, N, NDIM, XX, XY, WX, WY, L, LM, LN, LMDIM, A, X) 1 - COEFFICIENT MATRIX RETURNED. KNOT SEQUENCE ALONG X AXIS KNOT SEQUENCE ALONG Y AXIS ORDER OF SPLINE ALONG X AXIS ORDER OF SPLINE ALONG Y AXIS NUMBER OF B-SPLINES IN X DIRECTION NUMBER OF B-SPLINES IN Y DIRECTION SIZE OF GAMMA, AND ONE SIDE OF X (NDIM > MAX(M, X VALUES OF DATA POINTS. Y VALUES OF DATA POINTS. WEIGHTS ALONG X DIRECTION. WEIGHTS ALONG Y DIRECTION. DATA MATRIX. NO. OF DATA POINTS ALONG X AXIS NO. OF DATA POINTS ALONG Y AXIS FIRST DIMENSION OF DATA MATRIX WORK SPACE FOR GENERATING SPLINE WORK SPACE FOR GENERATING COEFFICIENT MATRIX. REAL TX(1), TY(1), XX(1), XY(1), A ( l ) , L(LMDIM,LN) , X(NDIM,LN), GAMMA(NDIM,NDIM), VNIKX(20), WX(1), WY(1) ZERO MATRICES 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. CONTINUE 63 20 C DO 20 1=1 ,M DO 20 J = 1, LN X ( I , J ) = 0. CONTINUE ILEFT = KX IMK = 0 DO 70 L I = 1, LM C C C 30 C C C C 40 C C C C 50 C C C C 60 70 C C C C C FIND INTERVAL "TX(ILEFT) <= X X ( L I ) < TX(ILEFT+1)" I F (ILEFT .EQ. M) GO TO 40 I F (XX(LI) .LT. TX(ILEFT + 1)) GO TO 40 ILEFT = ILEFT + 1 IMK = ILEFT - KX GO TO 30 PLACE VALUE OF ALL NON ZERO BSPLINES EVALUATED ON X X ( L I ) IN VNIKX CALL BSPLVN(TX, KX, 1, X X ( L 1 ) , ILEFT, VNIKX) GENERATE T A L AND PLACE IN X 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 ( L 1 , J ) CONTINUE GENERATE T A A AND PLACE IN A DO 60 M l = J J , KX J = IMK + M l U K = IJ(J,I,KX) A ( I J K ) = A ( I J K ) + WXBSPL * VNIKX (Ml) CONTINUE CONTINUE SOLVE T (J) T (J) A A X = A L TO OBTAIN THE MATRIX X 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 90 C C C 95 100 C NRHS = 2 CONTINUE CLEAR MATRICES DO 95 I=1,IMAX ' A ( I ) = 0. CONTINUE DO 100 I = 1, MNMAX DO 100 J = 1, MNMAX GAMMA(I,J) = 0. CONTINUE ILEFT = KY IMK = 0 DO 150 L I = 1, LN C C C 110 C C C 120 C C C C 130 C C C C 140 150 C C C C FIND THE INTERVAL "TY (ILEFT) <= XY(L1) < TY (ILEFT+1) 11 I F (ILEFT .EQ. N) GO TO 120 I F (XY(L1) .LT. TY(ILEFT + 1 ) ) GO TO 120 ILEFT = ILEFT + 1 IMK = ILEFT - KY GO TO 110 OBTAIN VALUE OF ALL NON ZERO B SPLINES OF XY(L1) CALL BSPLVN(TY, KY, 1, X Y ( L 1 ) , ILEFT, VNIKX) GENERATE T T B X AND PLACE IN GAMMA 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 , L 1 ) CONTINUE GENERATE T B B AND PLACE I N A DO 140 M l = J J , KY J = IMK + M l U K = IJ(J,I,KY) A ( I J K ) = A ( I J K ) + WYBSPL * VNIKX(Ml) CONTINUE CONTINUE SOLVE T T(J) T B B GAMMA = B X NRHS = 1 RATIO = 1.0E-10 T TO OBTAIN GAMMA AND PLACE IN GAMMA 65 170 DO 170 I = 1, M CALL FBAND(A,GAMMA(1,1),N,KY,NRHS,RATIO,DET,JEXP,NSCALE, 1 ITER,EPS) NRHS = 2 CONTINUE CALL COPY(GAMMA,TX,TY,KX,KY,M,N,NDIM) RETURN END C REAL FUNCTION SPL2D(X, Y) C C C C C C C C THIS FUNCTION EVALUATES THE B-SPLINE FUNCTION DEFINED BY THE COEFFICIENT MATRIX GAMMA AND THE KNOT SPQUENCES TX AND TY COMMON /CA/ GAMMAT(30,30), T X ( 3 5 ) , T Y ( 3 5 ) , M, NDIM, N, KX, KY REAL BCOEF(20) SPL2D(X,Y)=SUM (SUM GAMMA(I,J) B I J J,K,T (Y)) B (X) I,H,S CALL INTERV(TX, M + 1, X, LEFTX, MFLAG) SPL2D1 = 0 . I F (MFLAG .NE. 0) GO TO 20 C C C 10 C 20 INNER LOOP FIRST DO 10 J = 1, KX BCOEF(J) = BVALUE(TY,GAMMAT(1,LEFTX - KX + J),N,KY,Y,0) CONTINUE SPL2D = BVALUE(TX(LEFTX - KX + 1),BCOEF,KX,KX,X,0) 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
Page Metadata
Item Metadata
Title | Interactive least squares surface fitting |
Creator |
Samsom, Anthony Harm |
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 |
Aggregated Source Repository | 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}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0051811/manifest