UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Interactive least squares surface fitting Samsom, Anthony Harm 1980

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

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

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  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo 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

Comment

Related Items