- Email: [email protected]
Two-dimensional interpolation of irregularly spaced data using polynomial splines
180
Physics of the Earth andPlanetary Interiors, 12 (1976) 180—187 © Elsevier Scientific Publishing Company, Amsterdam — Printed in The Netherlands
TWO-DIMENSIONAL INTERPOLATION OF IRREGULARLY SPACED DATA USING POLYNOMIAL SPLINES J.C.DOOLEY Bureau of Mineral Resources, Geology and Geophysics, Canberra, A. C.T. (Australia) (Accepted for publication February 3, 1976) Dooley, J.C., 1976. Two-dimensional interpolation of irregularly spaced data using polynomial splines. Phys. Earth Planet. Inter., 12: 180—1 87. A method has been developed for interpolating field data specified at irregularly spaced points over an area, such as is often the case with geophysical observations. The area is divided into triangles, of which the data points form the vertices. The interpolating function consists of a polynomial within each triangle. Across the boundary of each adjacent pair of triangles the values of the polynomials and their derivatives are matched to a prescribed degree. Third-degree polynomials permit matching of second derivatives, and fifth-degree polynomials are required for matching second derivatives. Boundary conditions may be imposed on the derivatives, and equations may be introduced to improve the smoothness of the interpolating function. The analytical nature of the function permits certain manipulations, for example calculation of directions of rays and intersection points in refraction seismology. The function could also be used for automatic contouring.
I. Introduction Since about 1960, polynomial spline functions have been used successfully for interpolation of data specifled at points spaced regularly or irregularly along a line (see e.g., Ahlberg et al., 1967). The most commonly used form consists of a different cubic polynomial in each interval with continuous first and second derivatives at the data points or nodes; the third derivative may be discontinuous at the nodes, In extending splines to two dimensions, most attention has been devoted to data specified at points forming a rectangular grid (e.g., Ahlberg et al., 1967; Bhattacharya, 1969). Bicubic polynomials were used in each rectangle; first derivatives, and some but not all second derivatives, were continuous at the nodes and edges of the rectangles. Hessing et al. (1972) extended bicubic interpolation to irregularly spaced data by drawing lines through the data points to form quadrangles; extra points are needed at some intersections of these lines in order to complete the quadrangles, and values were calculated from neighbouring data for these points before beginning the interpolation. Briggs (1974) developed a method for contouring
irregularly spaced data, based on minimizing the curvature f(V2S~2of the spline function; this is equivalent to putting V4S = 0 everywhere except at the nodes. Any finite set of points in the plane can be connected by a mesh of which each element is a triangle; the points form the vertices of these triangles, and no additional intersecting points are introduced by crossing edges. This suggests interpolation by a spline function defmed as a polynomial in each triangle. Hall (1969) divided a rectangular mesh into triangles by diagonals, and showed that bicubic polynomials could not be used for second-order spline functions on this mesh. Splines based on a triangular mesh have been used for approximating the solutions of certain differential equations, e.g. by Birkhoff and Mansfield (1974), who give references to previous work. Some of these splines use certain rational functions as well as tricubic polynomials (i.e., a polynomial which is cubic on the edges of the triangle). In the following, we investigate the use of poiynomial spline functions on a triangular mesh to interpolate irregularly spaced data, giving globally continuous derivatives to first or second order.
181
2. Statement of the problem Let (x,y) be coordinates in a cartesian system in the plane. Let {x1,y1: I = 1, V}be data points, at each ...,
of which a number g is specified. We require a function S(x,y) with continuous derivatives to order l(i.e., S E C1) such that: S(x~,y 1) g1, = (1) ~,
...,
We seek a solution by joining pairs of the V points byE non-intersecting line intervals which form the edges ofT triangles. In each triangle T1 (I = 1, 7), we put: ...,
n restricted to polynomials of degree m
so that: N~= [~(n + lXn + 2) 3(n m)] T (4) 4S = 0 inside We also consider polynomials satisfying V each triangle; these must haven = 4 with the additional restriction: 4P V 1(x,y) = 24a40 + 8a22 + 24a04 = 0 For these, N~= 14T. —
—
2.1. Homogeneous coordinates S(x,y) = P1.(x,y) (2) where F1 is a polynomial in x andy of degree n. We bear in mind the following requirements: (1) The function S(x,y) should have at least as many parameters (i.e., coefficients of polynomials) as there are relations to be satisfied, in order to ensure continuity of derivatives to order 1. (2) Any excess of parameters over relations wifi enable boundary conditions to be imposed and/or smoothness conditions to be satisfied (i.e., minimization of certain second or higher deriv~tives,or of functions or combinations of these). (3) The degree of the polynomial should not be too high, as this may permit unwanted oscifiations between the data points. (4) S(x,y) should be of a uniform character, e.g. independent of azimuth or choice of coordinate axes as far as possible, and using the same degree of polynomial for eachP~. Most attention to date has been given (e.g., Birkhoff, 1969) to bicubic polynomials, i.e. polynomials of the form: Pj(x,y)= ~afkx1yk,
0~Ø,k)~3
It is convenient to use homogeneous triangular coordinates, z1, z2, z3, defined as the perpendicular distances of a point from the three edges of a triangle (Fig. 1). Let the triangle A 1A2A3 (Fig. 1) be numbered sequentially in a clockwise sense, with edges of length o~opposite vertex A1. Let q1 be the length of the perpendicular from A. to the opposite edge, and a~be the mternal angle at A1, taken as positive. Let R be the point’(T1 ,z2,z3). Any two of the z1 are sufficient to specify R; they satisfy:
E z1/q1
1=1
=
1
(5)
Ai (x1y)r(q,o,o)
y
a
(3)
These reduce to cubic polynomials along lines parallel to the x ory axis, but in general are sixth-degree polynomials in a variable u along a direction oblique to these axes. Thus there is an inherent difference between directions parallel to the coordinate axes and other directions. We consider here two-dimensional polynomials of degree n, i.e. of form (3) but with 0 ~j + k ~ n without other restriction (which have (n + lXn + 2)/2 parameters each), and two-dimensional polynomials of degree
82/
~-~--~
—~‘
~ 83
a~
a,
~,
~
A,
A2(x,,y,) =(oq,,o) 82
(x,,y,) (o,o,q,)
~x Fig. 1. Homogeneous triangular coordinates.
182
Let the edge A,A1 of the triangle make an angle 0k’ where k = 6 (i + /)~ clockwise from they axis. Then the triangular coordinates of a point (x,y) are derived from:
A,
—
zk=(x—xf)cosOk—(y—yI)sinOk
Ti
(6) Zi
where x.,y, are the. coordinates of the vertex A,.
zk
is
Z 3
positive for points inside the triangle ifA1A1 is in a clock-wise direction. twoWe sides are collinear coincident. assume that the or triangle is not degenerate, i.e. no The polynomial: ~3a~~xmyk,
A, A,
2
Z2
u
__________
1
B, /“~ z3
Ti’
n>~(m+k)>~0
a
A, A2’
~‘
A,’
can be transformed to a homogeneous form:
Fig. 2. Relations across the common boundary of two adjacent
triangles. P(z 1,z2,z3)=
E
l+m +kn
clmkzllz2mzak
(7) 3. Equations for matching derivatives
There are the same number of coefficients clmk as there were amk in the original polynomial. The degree of the polynomial will remain the same for any such transformation. This symmetrical form has the advantage that the polynomial can be evaluated at any vertex (z1 = z1 = zk = 0), or along any edge of the triangle (z1 = 0). Also the condition that a point lies inside the triangle is simplyz1>0,i 1,2,3. NowA.let(see B. be point in the sideand A/Ak vertex Fig.a 2, where i = 3), let opposite u 1 be the distance A1B, where / = i mod 3 + 1. Let v1 be the distance AkBI, so that Ui + Uf = Then for B1 we have z1 = 0, and hence:
F(q1 ,0,0) = c~00q~” = gj and hence: 7’ (11) cnOO = g1/q1 with similar expressions for C0~0,C00~. For continuity across the common edge of two tnangles, two polynomials Q, Q’ of the form (9) must be equated for all u = v’, v = u’ along the edge (Fig. 2).
zk
This + 1 coefficients, but11 twoatrelations existinvolves betweenmthem by virtue of eq. the end
=
u1 sin
(8)
u1 sin ak
I
For each triangle, three equations are derived from expression (1), e.g.:
points; hence m 1 additional equations are introduced for each common edge(where m is the degree of the edge polynomials: m —
By these transformations, e.g.,F(z1,z2,za)1z3=o can be expressed in the form:
3.2. First order
E bjkua’u3” k+1m
(9)
bik = clko(cosec a2~(cosec~
(10)
Q(u3,v3)
3.1. Zero order
where:
with similar expressions for the polynomials along the other sides of the triangle.
The gradient at vertex V1 can be represented by two components G1~= dS/dx, G1~= dS/dy, which are taken as unicnown. For each edge meeting at V,,, there is an equation of the form: dQ
=
u0 orp
cos 0 + G1~sin 0
(12)
183
For continuity, G,~and Gjy must have the same values for all edges meeting at V1. Since each edge
4. Triangulation of a set of points
joins two vertices, there wifi be 2E such equations. Elimination of the 2 V unknown quantities {G1~,G1~: / = 1, V} results in 2(E J/) relations between the coefficients of the polynomials ~Q and hence of {F}. As = Q’, derivatives of all orders with respect to u (or v) are also equal across common edges. If z3 and z’3 are the coordinates perpendicular to the edge (Fig. 2), then we need:
One of the main requirements is to select S(x,y) so that it has sufficient parameters to enable the required derivatives to be matched at vertices and edges; thus it is important to know the numbers of triangles and edges in a triangulation of a set of V points. Any network of polygons enclosed in and covering a bounded polygon in the plane satisfies Euler’s equation for a “sphere with one boundary” in the language of topology: F+V=E+l (16)
...,
—
},
Q
dFdI~
(13)
where F = the number of polygons (non-overlapping);
Now:
V = the number of vertices; and E = the number of edges. For a triangulated set of points, let the outer bound-
dFI dz~ z30
I
can be expressed as a polynomial in (u,v) of degree n 1. Thus n relations between the two sets of coefficients are needed, butatasthe theend-points, gradients have been matched by eq. 12 n 2already new relations are required for each common edge. Thus for matching first derivatives, we have: —
—
3T equations of type (11) (m— l)Ea equationsfromQ=Q’,whereE~isthe number of internal edges 2(E V) equations derived from eq. 12 (n 2)Ea equations derived from eq. 13
ary form a polygon with Eb sides and Vb = Eb vertices. Then we have in addition to eq. 16, replacing F by T (the number of triangles):2Ea + Eb (17) 3T = 2E Eb = E + Ea = where the suffix a is used to denote interior edges or vertices; and from eqs. 16 and 17 we can deduce further relations: —
—
T= V+ Va~~~2 E = 2V+ Va 3 —
(18) (19)
—
—
The total number is then: Dt3T2Va+(m+n_1)Ea
(14)
From eqs. 18 and 19 it is clear that for a given set of V points, once the number Vb of points in the boundary has been determined, then the total numbers of triangles andedgesareknown.
3.3. Second order
5. Parameters and equations
The second-order derivative or curvature is specified by three parameters:
The numbers of parameters N~determined from eq. 4 for various polynomials are given in Table I. Also listed are the numbers of equations D and D2have required for matching from eqs. 14 and 15.1 These been transformed by eqs. 16, 17, 18 and/or 19 into a form suitable for comparison with N~. The lowest-degree polynomial for which first den-
2S ddx2’
d2S dy2’
d2S dxdy
By arguments similar to those for the gradients, and noting that derivatives of the form d2P/dzdu are already equated across each side, we need 2E 3 V equations to match the second derivatives at the vertices, —
and n 3 equations across each internal edge. This gives a total of: —
D
4Va 2 = 3T
—
V+ (m + 2,z
—
2)Ea
(15)
vatives can be matched for any set of points is of degree 3 with no reduction in degree along the sides. This leaves V + Vb + 1 degrees of freedom which can be used for boundary or smoothing conditions, e.g., one boundary condition for each boundary edge, plus one relation at each vertex, plus one additional-general
184 TABLE I Parameters and equations for matching first- and second-order derivatives
n
m
32 3 4 4 4 4
32 2 4 3 2 V~P=0
5 5
5
N~_bi
N~
5
4 3
106T .7 15 12 9 14 21 18 15
6 6 6 6
6 5 4 3
28 25 22 19
5T_2+3Ea 5T_2+lEa
2 5 4 3
0
5
7 6 5
14 12 10 8
—2
19 16 13 10
relations. Other polynomials, such as n = 4, m = 2 for first derivatives, and n = 5, m = 4 for second derivatives, could be used for some but not all sets of points, as they depend on the relative numbers of internal and boundary points. These also have disadvantages with respect to the simple third- and fifth-degree polynomials as regards conditions (4) and (2) above. It is of interest however, that second derivatives could be matched by a cubic polynomial if: a
i.e.: Vb >(2V—
7T~V~4+3Ea
+4 5 +5 —3 3+2 0 1 —1 +3 4 +3 +3 9 +11 —9 7 +8 —6 5 -‘-5 —3
9 8 7 6
or
2
3 2Vb4lVa+l ~lVa~
relation; or two boundary conditions for each edge, plus one relation for each internal vertex, plus one additional general relation, For matching second derivatives the lowest degree polynomial is n = 5, m = 5; this permits four boundary conditions for each boundary edge, plus two relations for every vertex, or some equivalent combination of
Vb >2Va7
N~-D
7)13
2 6 5
4 6 9 8 7
—17 —14 —11 —8
12 11 10 9
lVb2Va+ 7 —1 —5 +10 3 —1 + 6 1—4+9 —1 7 +12 2 —3 +8 6 +2 + 3 4 —1 +6 2 —4 +9
10 8 6 4
+7 +4 +1 —2
—
+ + +
2 1 4 7
should be satisfactory e.g., for representation of a seismic refractor for two-dimensional ray-tracing or similar computations. We show how to reduce the problem from 1OT unknown parameters to 2 V parameters and the same number of equations. We have three equations of type (11), which give c 300, c030 and c003 directly in terms of g1 thus we can reduce the number of unknown parameters from lOT to 7T. From eq. 9, we get for the gradient at a vertex in the direction along an edge: dQ 2(3b =p 03 b1 2) (20) u0 b03 is related to c030 by eq. 1Q;hence eq. 12 expresses b12 in terms of G~,Gy and known quantities. For any three edges E1 (1 = 1, 2, 3) meeting at a common vertex V, let B, denote b12 or b21 according as u = 0 or v = 0 at V. Then eliminating G~and G~from the three equations of type (12) gives: 2B 2B (p, 1 —3g/p1)sin(03—02)+(,p2 2 3g/p2) sin(0j 03) + (p 2B 3 3 3g/p3) sin(02 0~)= 0 (21) —
—
—
2 points
This condition is satisfied by a square array of n ifn~<5. 6. Cubic polynomials It was decided in the first place to study third-degree polynomials, matching first derivatives; such a surface
—
—
—
-
For practical calculation, two edges most nearly at nght angles are selected as reference lines, and an equation of form (21) is written for these two edges together with each other edge in turn meeting at that .
185 vertex. We defme a vector B of dimension 2V, whose elements consist of the reference {B1}. From such equations at all vertices, the coefficients b12 and b21 in eq. 9 can be determined for every edge in terms the known we get 2Eof B2and V equations of quantities this nature.p. g and 0; thus —
The coefficients c012, c021, c102, c201, c120 and c210 for each triangle can then be determined by eq. 10 in terms of B. As the data values and gradients are matched at both ends of the edge for the polynomials F and F’ the edge triangles on eitherQside it follows thatinthe polynomials andofQ’the areedge, matched along theedge (Fig. 2). The transverse gradient at an edge is given by expressions such as:
3
z3=0
=
(c2o1
—
c210 cos a1
—
3c300 cos a2)z1 2
+(c111—2c12ocosa1—2c210cosa2)z1z2 2 +(c021 c030 cos ai c120 cos a2)z2 (22) —
—
As the gradients are matched at the ends of each 2 and z1 forF = 0;z2 = 0) zedge, 2 are(i.e., equal and F’.the Wecoefficients write Ba = bofz1 2 12 = b~1 andBb=bai b~2andC=c111thenmatchingthe coefficients of z1z2 and z~z’2gives: Csm a1 sin aa2 + C sm a1 sin a2 = 2Ba(cot 1 + cot a~)+ 2Bb(cot a2 + cot a~) (23) AS Ba and Bb can be expressed in terms of the vector B, the RHS of eq. 23 can be transformed to an expression involving elements of B as the only unknown parameters. Together with C1, i = 1, ..., T, we now have 2V + Tleaving unknown with Ea them, 2 Vbparameters, + Va + 1 degrees of relations freedom.between 6.1. Boundary conditions The surplus degrees of freedom may be2S)2 usedcould to minimize the curvature of S(x,y), i.e., f(v be minimized, where the integral is taken over the area covered by the triangles. This involves use of Lagrangian parameters or generalized inverses, and increases the size of the problem substantially. This approach will be discussed elsewhere. Another approach is to devise Va + 21’i, + 1 extra
equations which would tend to reduce but not necessarily minimize the curvature; there are then exactly enough equations to determine the required parameters. 2S = 0 Following Briggs,of it the is desirable have along the boundary polygon.toAs V2PV 1 is a linear function in the homogeneous coordinates (u,v) along any edge of a triangle, and hence along any boundary edge 2P of the polygon, we can achieve this by putting V 1 2Vb = 0 at both endsleaving of eacha boundary equations, further Va edge; + 1 tothis be gives found. N ow. V2P=h 1z1 +h2z2 +h3z3 (24) where: h1=6c3o0+2c120+2c,02—2c111cosa1 h
4ccosa4ccosa 2
=
6C~j30+ 2c210 + 2c0i2
2c111 cos a2
—
+ -,
— —
—
‘-“-003
4c102 cos
—
-~
~.C2Ol +
a2
—
—
4c021 C05 Cl1
25)
4c120 cos a3 A
~“—012
cos a1
2c1 ii cos a3
If the edge z, = 0 is a boundary edge, then the bothsdary conditions are: h2
=
h3
=
0
(26)
It is possible that a boundary point is a vertex of only one triangle. In this case, putting c~p1= 0 at this vertex gives only one equation instead of the two normally obtained for a boundary vertex. Toorcompensate 2f/du2 d2f/du2 for for this, it is proposed to equate d the two edges where they meet at this point; this means that the principal axes of curvature will bisect the angles formed by these edges. All the coefficients except c 1 ii can be expressed in terms of for B; hence either h2 = 0 or = 0 gives an expression c 111 in terms ofB; we also get an additional relation between elements ofB. The boundary edges give a total of Vb such relations, and reduce the number of unknowns by Vb. We can express the remaining C, in terms of B by the following procedure. We arrange the triangles in sets (Fig. 3) defined inductively as follows:
186
pliers or a similar technique would have to be used. On 2
2
2 2 2
I
the principle that, in minimizing a sum of squares of linear quantitiesthe first move is to make average equal to zero, the approach adopted here istheir to equate the average (or the sum, which comes to the same thing)
I
4
-
2
~ 5 3
6
5 6 6
5
of the values of V2P1 at each internal vertex to zero. 2F We still need one equation. Somewhat arbitrarily, it was decided to equate the sum of values of V 1 at all boundary vertices to zero. Since the values for boundary triangles (i.e., triangles of set 1) are already zero, this means the sum of the values for triangles with just one vertex on the boundary.
I
~-
-
2
2
4
6 6
3
2 4
4
5
4
•
5 4
3
32
I
2
7. Conclusions 2
2
Fig. 3. Ordering of triangles in sets to facilitate reduction of equations for coefficients C.
Set I consists of triangles with an edge forming part of the boundary polygon. For k> 1, set k consists of triangles which do not belong to any set m(m
—
It has been shown that a third-degree polynomial is capable of matching first derivatives, and that the problem may be reduced to solving 2 V equations for 2 V parameters. By similar procedures, a fifth-degree poly.
nomial may be used for matching second derivatives; the number of equations and parameters for this case is 5 V. An alternative approach is to minimize the integrated least-squares total curvature or differential curvature. Also, 2S/dz2 different conditions may be applied, = 0boundary across boundary edges instead of eq. 26. e.g. d It is intended to discuss these approaches elsewhere. The method is applicable where a smooth interpolant is required which can be calculated readily at intermediate points. The function S(x,y) could be used for drawing contours, but it has not yet been demonstrated whether it would have any advantage over other methods such as that of Briggs (1974).
—
—
Acknowledgement The permission of the Director, Bureau of Mineral Resources, Geology and Geophysics, Australia, to publish this paper is acknowledged.
References J.H., Nilson, E.N. and Walsh, J.L, 1967. The Theory of Splines and their Applications. Academic Press, New York, N.Y., 284 pp. Bhattacharya, B.K., 1969. Bicubic spline interpolation as a Ahlberg,
187 method for treatment ofpotential field data. Geophysics, 34: 402—423. Birkhoff, G., 1969. Piecewisebicubic interpolation and approximation in polygons. In: I.J. Schoenberg (Editor), Approximations with Special Emphasis on Spline Functions, Academic Press, New York, N.Y., pp. 185—221. Birkhoff, G. and Mansfield, L., 1974. Compatible triangular fmite elements. J. Math. Anal. Appl., 47: 531—553.
Briggs, I.C., 1974. Machine contouringusing minimum curvature. Geophysics, 39: 39—48. Hall, (IA., 1969. Bicubic interpolation over triangles. J. Math. Mech., 19: 1—11. Hessing, R.C., Lee, H.K., Pierce, A., and Powers, E.N., 1972. Automatic contouring using bicubic functions. Geophysics, 37: 669—674.
Physics of the Earth andPlanetary Interiors, 12 (1976) 180—187 © Elsevier Scientific Publishing Company, Amsterdam — Printed in The Netherlands
TWO-DIMENSIONAL INTERPOLATION OF IRREGULARLY SPACED DATA USING POLYNOMIAL SPLINES J.C.DOOLEY Bureau of Mineral Resources, Geology and Geophysics, Canberra, A. C.T. (Australia) (Accepted for publication February 3, 1976) Dooley, J.C., 1976. Two-dimensional interpolation of irregularly spaced data using polynomial splines. Phys. Earth Planet. Inter., 12: 180—1 87. A method has been developed for interpolating field data specified at irregularly spaced points over an area, such as is often the case with geophysical observations. The area is divided into triangles, of which the data points form the vertices. The interpolating function consists of a polynomial within each triangle. Across the boundary of each adjacent pair of triangles the values of the polynomials and their derivatives are matched to a prescribed degree. Third-degree polynomials permit matching of second derivatives, and fifth-degree polynomials are required for matching second derivatives. Boundary conditions may be imposed on the derivatives, and equations may be introduced to improve the smoothness of the interpolating function. The analytical nature of the function permits certain manipulations, for example calculation of directions of rays and intersection points in refraction seismology. The function could also be used for automatic contouring.
I. Introduction Since about 1960, polynomial spline functions have been used successfully for interpolation of data specifled at points spaced regularly or irregularly along a line (see e.g., Ahlberg et al., 1967). The most commonly used form consists of a different cubic polynomial in each interval with continuous first and second derivatives at the data points or nodes; the third derivative may be discontinuous at the nodes, In extending splines to two dimensions, most attention has been devoted to data specified at points forming a rectangular grid (e.g., Ahlberg et al., 1967; Bhattacharya, 1969). Bicubic polynomials were used in each rectangle; first derivatives, and some but not all second derivatives, were continuous at the nodes and edges of the rectangles. Hessing et al. (1972) extended bicubic interpolation to irregularly spaced data by drawing lines through the data points to form quadrangles; extra points are needed at some intersections of these lines in order to complete the quadrangles, and values were calculated from neighbouring data for these points before beginning the interpolation. Briggs (1974) developed a method for contouring
irregularly spaced data, based on minimizing the curvature f(V2S~2of the spline function; this is equivalent to putting V4S = 0 everywhere except at the nodes. Any finite set of points in the plane can be connected by a mesh of which each element is a triangle; the points form the vertices of these triangles, and no additional intersecting points are introduced by crossing edges. This suggests interpolation by a spline function defmed as a polynomial in each triangle. Hall (1969) divided a rectangular mesh into triangles by diagonals, and showed that bicubic polynomials could not be used for second-order spline functions on this mesh. Splines based on a triangular mesh have been used for approximating the solutions of certain differential equations, e.g. by Birkhoff and Mansfield (1974), who give references to previous work. Some of these splines use certain rational functions as well as tricubic polynomials (i.e., a polynomial which is cubic on the edges of the triangle). In the following, we investigate the use of poiynomial spline functions on a triangular mesh to interpolate irregularly spaced data, giving globally continuous derivatives to first or second order.
181
2. Statement of the problem Let (x,y) be coordinates in a cartesian system in the plane. Let {x1,y1: I = 1, V}be data points, at each ...,
of which a number g is specified. We require a function S(x,y) with continuous derivatives to order l(i.e., S E C1) such that: S(x~,y 1) g1, = (1) ~,
...,
We seek a solution by joining pairs of the V points byE non-intersecting line intervals which form the edges ofT triangles. In each triangle T1 (I = 1, 7), we put: ...,
n restricted to polynomials of degree m
so that: N~= [~(n + lXn + 2) 3(n m)] T (4) 4S = 0 inside We also consider polynomials satisfying V each triangle; these must haven = 4 with the additional restriction: 4P V 1(x,y) = 24a40 + 8a22 + 24a04 = 0 For these, N~= 14T. —
—
2.1. Homogeneous coordinates S(x,y) = P1.(x,y) (2) where F1 is a polynomial in x andy of degree n. We bear in mind the following requirements: (1) The function S(x,y) should have at least as many parameters (i.e., coefficients of polynomials) as there are relations to be satisfied, in order to ensure continuity of derivatives to order 1. (2) Any excess of parameters over relations wifi enable boundary conditions to be imposed and/or smoothness conditions to be satisfied (i.e., minimization of certain second or higher deriv~tives,or of functions or combinations of these). (3) The degree of the polynomial should not be too high, as this may permit unwanted oscifiations between the data points. (4) S(x,y) should be of a uniform character, e.g. independent of azimuth or choice of coordinate axes as far as possible, and using the same degree of polynomial for eachP~. Most attention to date has been given (e.g., Birkhoff, 1969) to bicubic polynomials, i.e. polynomials of the form: Pj(x,y)= ~afkx1yk,
0~Ø,k)~3
It is convenient to use homogeneous triangular coordinates, z1, z2, z3, defined as the perpendicular distances of a point from the three edges of a triangle (Fig. 1). Let the triangle A 1A2A3 (Fig. 1) be numbered sequentially in a clockwise sense, with edges of length o~opposite vertex A1. Let q1 be the length of the perpendicular from A. to the opposite edge, and a~be the mternal angle at A1, taken as positive. Let R be the point’(T1 ,z2,z3). Any two of the z1 are sufficient to specify R; they satisfy:
E z1/q1
1=1
=
1
(5)
Ai (x1y)r(q,o,o)
y
a
(3)
These reduce to cubic polynomials along lines parallel to the x ory axis, but in general are sixth-degree polynomials in a variable u along a direction oblique to these axes. Thus there is an inherent difference between directions parallel to the coordinate axes and other directions. We consider here two-dimensional polynomials of degree n, i.e. of form (3) but with 0 ~j + k ~ n without other restriction (which have (n + lXn + 2)/2 parameters each), and two-dimensional polynomials of degree
82/
~-~--~
—~‘
~ 83
a~
a,
~,
~
A,
A2(x,,y,) =(oq,,o) 82
(x,,y,) (o,o,q,)
~x Fig. 1. Homogeneous triangular coordinates.
182
Let the edge A,A1 of the triangle make an angle 0k’ where k = 6 (i + /)~ clockwise from they axis. Then the triangular coordinates of a point (x,y) are derived from:
A,
—
zk=(x—xf)cosOk—(y—yI)sinOk
Ti
(6) Zi
where x.,y, are the. coordinates of the vertex A,.
zk
is
Z 3
positive for points inside the triangle ifA1A1 is in a clock-wise direction. twoWe sides are collinear coincident. assume that the or triangle is not degenerate, i.e. no The polynomial: ~3a~~xmyk,
A, A,
2
Z2
u
__________
1
B, /“~ z3
Ti’
n>~(m+k)>~0
a
A, A2’
~‘
A,’
can be transformed to a homogeneous form:
Fig. 2. Relations across the common boundary of two adjacent
triangles. P(z 1,z2,z3)=
E
l+m +kn
clmkzllz2mzak
(7) 3. Equations for matching derivatives
There are the same number of coefficients clmk as there were amk in the original polynomial. The degree of the polynomial will remain the same for any such transformation. This symmetrical form has the advantage that the polynomial can be evaluated at any vertex (z1 = z1 = zk = 0), or along any edge of the triangle (z1 = 0). Also the condition that a point lies inside the triangle is simplyz1>0,i 1,2,3. NowA.let(see B. be point in the sideand A/Ak vertex Fig.a 2, where i = 3), let opposite u 1 be the distance A1B, where / = i mod 3 + 1. Let v1 be the distance AkBI, so that Ui + Uf = Then for B1 we have z1 = 0, and hence:
F(q1 ,0,0) = c~00q~” = gj and hence: 7’ (11) cnOO = g1/q1 with similar expressions for C0~0,C00~. For continuity across the common edge of two tnangles, two polynomials Q, Q’ of the form (9) must be equated for all u = v’, v = u’ along the edge (Fig. 2).
zk
This + 1 coefficients, but11 twoatrelations existinvolves betweenmthem by virtue of eq. the end
=
u1 sin
(8)
u1 sin ak
I
For each triangle, three equations are derived from expression (1), e.g.:
points; hence m 1 additional equations are introduced for each common edge(where m is the degree of the edge polynomials: m —
By these transformations, e.g.,F(z1,z2,za)1z3=o can be expressed in the form:
3.2. First order
E bjkua’u3” k+1m
(9)
bik = clko(cosec a2~(cosec~
(10)
Q(u3,v3)
3.1. Zero order
where:
with similar expressions for the polynomials along the other sides of the triangle.
The gradient at vertex V1 can be represented by two components G1~= dS/dx, G1~= dS/dy, which are taken as unicnown. For each edge meeting at V,,, there is an equation of the form: dQ
=
u0 orp
cos 0 + G1~sin 0
(12)
183
For continuity, G,~and Gjy must have the same values for all edges meeting at V1. Since each edge
4. Triangulation of a set of points
joins two vertices, there wifi be 2E such equations. Elimination of the 2 V unknown quantities {G1~,G1~: / = 1, V} results in 2(E J/) relations between the coefficients of the polynomials ~Q and hence of {F}. As = Q’, derivatives of all orders with respect to u (or v) are also equal across common edges. If z3 and z’3 are the coordinates perpendicular to the edge (Fig. 2), then we need:
One of the main requirements is to select S(x,y) so that it has sufficient parameters to enable the required derivatives to be matched at vertices and edges; thus it is important to know the numbers of triangles and edges in a triangulation of a set of V points. Any network of polygons enclosed in and covering a bounded polygon in the plane satisfies Euler’s equation for a “sphere with one boundary” in the language of topology: F+V=E+l (16)
...,
—
},
Q
dFdI~
(13)
where F = the number of polygons (non-overlapping);
Now:
V = the number of vertices; and E = the number of edges. For a triangulated set of points, let the outer bound-
dFI dz~ z30
I
can be expressed as a polynomial in (u,v) of degree n 1. Thus n relations between the two sets of coefficients are needed, butatasthe theend-points, gradients have been matched by eq. 12 n 2already new relations are required for each common edge. Thus for matching first derivatives, we have: —
—
3T equations of type (11) (m— l)Ea equationsfromQ=Q’,whereE~isthe number of internal edges 2(E V) equations derived from eq. 12 (n 2)Ea equations derived from eq. 13
ary form a polygon with Eb sides and Vb = Eb vertices. Then we have in addition to eq. 16, replacing F by T (the number of triangles):2Ea + Eb (17) 3T = 2E Eb = E + Ea = where the suffix a is used to denote interior edges or vertices; and from eqs. 16 and 17 we can deduce further relations: —
—
T= V+ Va~~~2 E = 2V+ Va 3 —
(18) (19)
—
—
The total number is then: Dt3T2Va+(m+n_1)Ea
(14)
From eqs. 18 and 19 it is clear that for a given set of V points, once the number Vb of points in the boundary has been determined, then the total numbers of triangles andedgesareknown.
3.3. Second order
5. Parameters and equations
The second-order derivative or curvature is specified by three parameters:
The numbers of parameters N~determined from eq. 4 for various polynomials are given in Table I. Also listed are the numbers of equations D and D2have required for matching from eqs. 14 and 15.1 These been transformed by eqs. 16, 17, 18 and/or 19 into a form suitable for comparison with N~. The lowest-degree polynomial for which first den-
2S ddx2’
d2S dy2’
d2S dxdy
By arguments similar to those for the gradients, and noting that derivatives of the form d2P/dzdu are already equated across each side, we need 2E 3 V equations to match the second derivatives at the vertices, —
and n 3 equations across each internal edge. This gives a total of: —
D
4Va 2 = 3T
—
V+ (m + 2,z
—
2)Ea
(15)
vatives can be matched for any set of points is of degree 3 with no reduction in degree along the sides. This leaves V + Vb + 1 degrees of freedom which can be used for boundary or smoothing conditions, e.g., one boundary condition for each boundary edge, plus one relation at each vertex, plus one additional-general
184 TABLE I Parameters and equations for matching first- and second-order derivatives
n
m
32 3 4 4 4 4
32 2 4 3 2 V~P=0
5 5
5
N~_bi
N~
5
4 3
106T .7 15 12 9 14 21 18 15
6 6 6 6
6 5 4 3
28 25 22 19
5T_2+3Ea 5T_2+lEa
2 5 4 3
0
5
7 6 5
14 12 10 8
—2
19 16 13 10
relations. Other polynomials, such as n = 4, m = 2 for first derivatives, and n = 5, m = 4 for second derivatives, could be used for some but not all sets of points, as they depend on the relative numbers of internal and boundary points. These also have disadvantages with respect to the simple third- and fifth-degree polynomials as regards conditions (4) and (2) above. It is of interest however, that second derivatives could be matched by a cubic polynomial if: a
i.e.: Vb >(2V—
7T~V~4+3Ea
+4 5 +5 —3 3+2 0 1 —1 +3 4 +3 +3 9 +11 —9 7 +8 —6 5 -‘-5 —3
9 8 7 6
or
2
3 2Vb4lVa+l ~lVa~
relation; or two boundary conditions for each edge, plus one relation for each internal vertex, plus one additional general relation, For matching second derivatives the lowest degree polynomial is n = 5, m = 5; this permits four boundary conditions for each boundary edge, plus two relations for every vertex, or some equivalent combination of
Vb >2Va7
N~-D
7)13
2 6 5
4 6 9 8 7
—17 —14 —11 —8
12 11 10 9
lVb2Va+ 7 —1 —5 +10 3 —1 + 6 1—4+9 —1 7 +12 2 —3 +8 6 +2 + 3 4 —1 +6 2 —4 +9
10 8 6 4
+7 +4 +1 —2
—
+ + +
2 1 4 7
should be satisfactory e.g., for representation of a seismic refractor for two-dimensional ray-tracing or similar computations. We show how to reduce the problem from 1OT unknown parameters to 2 V parameters and the same number of equations. We have three equations of type (11), which give c 300, c030 and c003 directly in terms of g1 thus we can reduce the number of unknown parameters from lOT to 7T. From eq. 9, we get for the gradient at a vertex in the direction along an edge: dQ 2(3b =p 03 b1 2) (20) u0 b03 is related to c030 by eq. 1Q;hence eq. 12 expresses b12 in terms of G~,Gy and known quantities. For any three edges E1 (1 = 1, 2, 3) meeting at a common vertex V, let B, denote b12 or b21 according as u = 0 or v = 0 at V. Then eliminating G~and G~from the three equations of type (12) gives: 2B 2B (p, 1 —3g/p1)sin(03—02)+(,p2 2 3g/p2) sin(0j 03) + (p 2B 3 3 3g/p3) sin(02 0~)= 0 (21) —
—
—
2 points
This condition is satisfied by a square array of n ifn~<5. 6. Cubic polynomials It was decided in the first place to study third-degree polynomials, matching first derivatives; such a surface
—
—
—
-
For practical calculation, two edges most nearly at nght angles are selected as reference lines, and an equation of form (21) is written for these two edges together with each other edge in turn meeting at that .
185 vertex. We defme a vector B of dimension 2V, whose elements consist of the reference {B1}. From such equations at all vertices, the coefficients b12 and b21 in eq. 9 can be determined for every edge in terms the known we get 2Eof B2and V equations of quantities this nature.p. g and 0; thus —
The coefficients c012, c021, c102, c201, c120 and c210 for each triangle can then be determined by eq. 10 in terms of B. As the data values and gradients are matched at both ends of the edge for the polynomials F and F’ the edge triangles on eitherQside it follows thatinthe polynomials andofQ’the areedge, matched along theedge (Fig. 2). The transverse gradient at an edge is given by expressions such as:
3
z3=0
=
(c2o1
—
c210 cos a1
—
3c300 cos a2)z1 2
+(c111—2c12ocosa1—2c210cosa2)z1z2 2 +(c021 c030 cos ai c120 cos a2)z2 (22) —
—
As the gradients are matched at the ends of each 2 and z1 forF = 0;z2 = 0) zedge, 2 are(i.e., equal and F’.the Wecoefficients write Ba = bofz1 2 12 = b~1 andBb=bai b~2andC=c111thenmatchingthe coefficients of z1z2 and z~z’2gives: Csm a1 sin aa2 + C sm a1 sin a2 = 2Ba(cot 1 + cot a~)+ 2Bb(cot a2 + cot a~) (23) AS Ba and Bb can be expressed in terms of the vector B, the RHS of eq. 23 can be transformed to an expression involving elements of B as the only unknown parameters. Together with C1, i = 1, ..., T, we now have 2V + Tleaving unknown with Ea them, 2 Vbparameters, + Va + 1 degrees of relations freedom.between 6.1. Boundary conditions The surplus degrees of freedom may be2S)2 usedcould to minimize the curvature of S(x,y), i.e., f(v be minimized, where the integral is taken over the area covered by the triangles. This involves use of Lagrangian parameters or generalized inverses, and increases the size of the problem substantially. This approach will be discussed elsewhere. Another approach is to devise Va + 21’i, + 1 extra
equations which would tend to reduce but not necessarily minimize the curvature; there are then exactly enough equations to determine the required parameters. 2S = 0 Following Briggs,of it the is desirable have along the boundary polygon.toAs V2PV 1 is a linear function in the homogeneous coordinates (u,v) along any edge of a triangle, and hence along any boundary edge 2P of the polygon, we can achieve this by putting V 1 2Vb = 0 at both endsleaving of eacha boundary equations, further Va edge; + 1 tothis be gives found. N ow. V2P=h 1z1 +h2z2 +h3z3 (24) where: h1=6c3o0+2c120+2c,02—2c111cosa1 h
4ccosa4ccosa 2
=
6C~j30+ 2c210 + 2c0i2
2c111 cos a2
—
+ -,
— —
—
‘-“-003
4c102 cos
—
-~
~.C2Ol +
a2
—
—
4c021 C05 Cl1
25)
4c120 cos a3 A
~“—012
cos a1
2c1 ii cos a3
If the edge z, = 0 is a boundary edge, then the bothsdary conditions are: h2
=
h3
=
0
(26)
It is possible that a boundary point is a vertex of only one triangle. In this case, putting c~p1= 0 at this vertex gives only one equation instead of the two normally obtained for a boundary vertex. Toorcompensate 2f/du2 d2f/du2 for for this, it is proposed to equate d the two edges where they meet at this point; this means that the principal axes of curvature will bisect the angles formed by these edges. All the coefficients except c 1 ii can be expressed in terms of for B; hence either h2 = 0 or = 0 gives an expression c 111 in terms ofB; we also get an additional relation between elements ofB. The boundary edges give a total of Vb such relations, and reduce the number of unknowns by Vb. We can express the remaining C, in terms of B by the following procedure. We arrange the triangles in sets (Fig. 3) defined inductively as follows:
186
pliers or a similar technique would have to be used. On 2
2
2 2 2
I
the principle that, in minimizing a sum of squares of linear quantitiesthe first move is to make average equal to zero, the approach adopted here istheir to equate the average (or the sum, which comes to the same thing)
I
4
-
2
~ 5 3
6
5 6 6
5
of the values of V2P1 at each internal vertex to zero. 2F We still need one equation. Somewhat arbitrarily, it was decided to equate the sum of values of V 1 at all boundary vertices to zero. Since the values for boundary triangles (i.e., triangles of set 1) are already zero, this means the sum of the values for triangles with just one vertex on the boundary.
I
~-
-
2
2
4
6 6
3
2 4
4
5
4
•
5 4
3
32
I
2
7. Conclusions 2
2
Fig. 3. Ordering of triangles in sets to facilitate reduction of equations for coefficients C.
Set I consists of triangles with an edge forming part of the boundary polygon. For k> 1, set k consists of triangles which do not belong to any set m(m
—
It has been shown that a third-degree polynomial is capable of matching first derivatives, and that the problem may be reduced to solving 2 V equations for 2 V parameters. By similar procedures, a fifth-degree poly.
nomial may be used for matching second derivatives; the number of equations and parameters for this case is 5 V. An alternative approach is to minimize the integrated least-squares total curvature or differential curvature. Also, 2S/dz2 different conditions may be applied, = 0boundary across boundary edges instead of eq. 26. e.g. d It is intended to discuss these approaches elsewhere. The method is applicable where a smooth interpolant is required which can be calculated readily at intermediate points. The function S(x,y) could be used for drawing contours, but it has not yet been demonstrated whether it would have any advantage over other methods such as that of Briggs (1974).
—
—
Acknowledgement The permission of the Director, Bureau of Mineral Resources, Geology and Geophysics, Australia, to publish this paper is acknowledged.
References J.H., Nilson, E.N. and Walsh, J.L, 1967. The Theory of Splines and their Applications. Academic Press, New York, N.Y., 284 pp. Bhattacharya, B.K., 1969. Bicubic spline interpolation as a Ahlberg,
187 method for treatment ofpotential field data. Geophysics, 34: 402—423. Birkhoff, G., 1969. Piecewisebicubic interpolation and approximation in polygons. In: I.J. Schoenberg (Editor), Approximations with Special Emphasis on Spline Functions, Academic Press, New York, N.Y., pp. 185—221. Birkhoff, G. and Mansfield, L., 1974. Compatible triangular fmite elements. J. Math. Anal. Appl., 47: 531—553.
Briggs, I.C., 1974. Machine contouringusing minimum curvature. Geophysics, 39: 39—48. Hall, (IA., 1969. Bicubic interpolation over triangles. J. Math. Mech., 19: 1—11. Hessing, R.C., Lee, H.K., Pierce, A., and Powers, E.N., 1972. Automatic contouring using bicubic functions. Geophysics, 37: 669—674.