__ b!!id
&-ii!23 B
EUEVIER
Computer
Aided Geometric
Design 13
( 1996) 81-88
Scattered data interpolation and approximation using bivariate C1 piecewise cubic polynomials Ming-Jun Lai * Department of Mathematics, The University of Georgia, Athens, GA 30602, USA Received December
1993; revised January
1995
Abstract We show that if the scattered data over a polygonal domain can be quadrangulated, then the space of bivariate C’ piecewise cubic polynomial functions on a triangulation obtained from the quadrangulation has the full approximation order. We point out that our method is more efficient than the Clough-Tocher Kqvwordst Bivariate
scheme.
splines; B-net; Full approximation
order; Quadrangulation:
Scattered
data interpolation
1. Introduction
For a scattered data piecewise polynomial Usually, we triangulate To be precise, let A = S;(A)
set over a polygonal domain D, we would like to have a smooth surface over D to interpolate or approximate the given data. the data and construct a spline function over the triangulation. {ti} be a triangulation of a polygonal domain D c I@ and let
:= {s E C’(D):
s(~, E Pd,‘v’ti E A)
be the spline function space of smoothness r and degree d, where ti denotes a triangle in A and IEDddenotes the space of polynomials of total degree < d. These spaces are very useful for applications, e.g., numerical solution of PDE’s by using finite element methods and computer aided geometric design. To measure how well these spaces can be used to approximate unknown functions, we need a notion of approximation order which
* Email:
[email protected].
Supported
by the National
Science Foundation
0167-8396/96/$15.00 @ 1996 Elsevier Science B.V. All rights reserved SS’D/Ol67-8396(95)00007-O
under grant DMS-9303 12 I.
82
M.-J. Lb/Computer
Aided Geometric Design 13 (1996) 81-88
is defined as follows: the largest integer n obtainable and the spline space S;(n) satisfies dist(f,
S;(n))
such that the distance
between
f
6 C,lnl”
for all sufficiently smooth functions f is called the approximation order of Si( A). It is now well-known that the approximation order of SL (A) is d + 1 if d > 3r f2 as shown in (de Boor and Hiillig, 1988). The approximation order can be achieved by using a linear interpolatory operator described in (Chui and Lai, 1990). In practice, spline functions of lower degree are preferred due to their computational efficiency. For example, bivariate C’ cubits are most interesting in applications. However, for any given triangulation A, it is not known if we can construct a C’ piecewise cubic polynomial function based on n to interpolate any given unknown function f on the vertices of D although several algorithms have been described to do this job without knowing the existence of such interpolant; see (Gmelig Meyling, 1987; Gmelig Meyling and Pfluger, 1990; Grandine, 1987). If n is refined by subdividing each triangle of n into three subtriangles at any interior point of the triangle, then we can always construct a C’ piecewise cubic polynomial function to interpolate the function values and its derivatives of first order at the vertices of A. Such construction of C’ piecewise cubic polynomial functions is usually called Clough-Tocher’s scheme, see (Ciarlet, 1974). However, it is easy to see that the Clough-Tocher scheme increases the computational complexity by tripling the number of the triangles in A. The approximation order of SA(A) is dependent on the triangulation since de Boor and Hollig ( 1983) showed that the approximation order of S:(n) is 3 when n is a three-directional mesh while the author in (Lai, 1994) showed that approximation order of Si( A) is 4 when n is a four-directional mesh. Those features show bivariate C’ piecewise cubic polynomial functions are difficult to study. Let us revisit this problem: We are given a set of scattered data over a polygonal domain, we are required to construct a C’ surface to interpolate the data, we are forced to use cubic polynomials in order to get computational efficiency, and we have to have the surface close to the given function if it is sufficiently smooth. What should we do? If we triangulate the data, we come back to a difficult problem mentioned above. How about if we quadrangulate the data? Here is the definition of quadrangulation. Definition 1. A set of convex quadrilaterals is called a quadrangulation 0 of a polygonal domain D if the union of all those quadrilaterals is D and the intersection of any two qi and qi from 0 is either empty or a their common vertex or a their common edge. A typical example of quadrangulation is rectangular grid of a rectangular domain. The following is another example of an arbitrary quadrangulation over the unit square with 46 scattered data. Of course, we may not be able to quadrangulate for some data, e.g., the data consists of the vertices of a non-convex quadrilateral. (If the data cannot be quadrangulated, we may mix quadrilaterals with triangles and then refine each quadrilateral and triangle to obtain a triangulation. See Remark 2 in the next section.) Let us suppose a quadrangulation of the given data is possible. Then we subdivide each quadrilateral of 0 into four triangles
1
M.-J. Lni/Computer
Aided Geometric Design 13 (1996) 81-88
Y
Y
x x
09
i
0.9
x3
x
x
x
x.
x
x
-
x
07
x
x x
41.
x
Y
x
06
x
0
x
x
4‘
x x
04
x
x
x
x
0
x x
02
x x
x
m
0
01
02
Fig.
x
x
x
Ol-
03
04
* 0.5
06
07
0.8
09
1
I. Scattered data and a quadrangulation.
by its two diagonals. This results a triangulation $. The advantages of quadrangulation are 1. We can always construct a C’ piecewise cubic polynomial function based on $ interpolating the function values and the derivatives of first order at the vertices of 0; 2. The spline space S: ($) has full approximation order; 3. L,ocally supported basis functions in S:(e) of the finite element type can easily constructed. If this space is used to solve a PDE, then the resulting system is sparse. If this space is used to find a C’ interpolant, no large linear system needs to be solved as in (Gmelig Meyling, 1987; Gmelig Meyling and Pfluger, 1990; Grandine, 1987) ;
84
M.-J. Lb/Computer
Aided Geomeiric Design 13 (1996) 81-88
4. Comparing with Clough-Tocher’s scheme, our method needs fewer locally supported basis functions (see Remark 1 in the next section) and fewer triangles (about 2/3 of the number of triangles of Clough-Tocher’s scheme). Thus, our method is more efficient than Clough-Tocher’s scheme; 5. Graphically displaying such a surface can be done by displaying the subsurface over each quadrilateral of 0. The number of quadrilaterals of $ is just half of the number of triangles of n if we triangulate the data. The main purpose of this note is to establish the above results. From the advantages mentioned above, we see that although quadrangulation is not a panacea, it gives the designer of interpolatory C’ cubic spline surface a useful tool.
2. Main results In this section, we use Btzier representations to express cubic polynomial pieces. (If the reader is not familiar with that representation, please refer to (Farin, 1986) or (de Boor, 1987).) We shall construct so-called vertex splines and edge spline to have a locally supported basis of Sj($). The definition of vertex splines and edge splines on a quadrangulation is as follows. Definition 2. A spline function is S;( 9) is called a vertex spline if it is supported on the union of all quadrilaterals which share a common vertex. A spline function in Si ($) is called an edge spline if it is supported on the union of two quadrilaterals which share a common edge. Next we introduce the notion of derivative relative to an edge of 0 as in (Chui and Lai, 1990). For each interior edge e = [u~,~,v,,J] of 0, there are two triangles in $ sharing e. Let us choose the one which has larger area if possible. Otherwise let us fix any one of them. Denote by (v,,i , v,,~, v,). If e is a boundary edge, then there is one triangle (v,,i , vc,2, v,) containing e. Then derivatives relative to e are defined by 0,” = (D,-,,,,)“‘(D,,*-,,,,)“‘,
a=
where the derivatives are taken within standard differential operator D”f(x,
y> := D;‘D~f(x,
(a1,a2)
E z:,
the triangle
(v,,i, v,,~,v,).
Denote
by Da the
y).
We are ready to construct locally supported vertex splines and edge splines. Let V = {v} be the collection of all vertices of 0, & = {e} be the collection of all edges of 0 and let Q = {Q} be the collection of all quadrilaterals of 0. For each v E V, we define three vertex splines Vu,r, y E { (0,O) , (0,l) , ( 1,O) } which satisfy the following conditions: D%,,(u) D;‘+&(v,,~
= &,u&,y. ) = 0,
a E {(O,O), (l,O), e E &;
(0, l)},
U. E v;
(1.1) (1.2)
M.-J. Lb/Computer
Aided Geometric Design 13 (1996) RI-RR
85
v4
Fig. 2. An illustration of B-coefficients on a quadrilateral.
v,.,,
(1.3)
E C’(R2).
Here and throughout, as usual, the symbols S,,, and IS,,, are the Kronecker For each e E I, we define one edge spline V, satisfying the following D”&(u)
=o,
D~.‘~“V,(O,,,)
aE = s,,,,
{(0,0),(1,0),(0,1)}, c = (v,,, ,v,.2) E E;
v> E C’(lR2).
uEl);
delta.
(11.1) (11.2) (11.3)
Now we show the existence of those vertex splines and edge splines. For fixed v and y E {(O,O), (l,O),(O, l)}, consider V”[Q for a quadrilateral Q with vertices VI,V~,V~, and VJ. The B-net of G/e is displayed as shown in Fig. 2. Without loss of generality, we may assume that ~2, us, v4 are not v. By (I. 1) and the C’ smoothness condition, the B-coefficients of %,,[a are zeros around those vertices. I3y (1.2) on edges (~2, ~3) and (4, ~4) and the C’ smoothness condition (1.3), the B-coefficients located inside the triangle (~2, ~9, v4) have to be zero. However, vi may be v. We use (I. 1) to determine (~1, a~, and (~3. By using (1.2) on edges (VI, ~4) and (vi. vz), p and y may be a non-zero value. Finally, the remaining four B-coefficients of Ka,vla, i.e., -- 13 1 i-h’
- Y 1 fh’
r+kP I+k’
r+kP (l+k)(l+h)
are obtained by using C’ smoothness condition, where k = Iv4 - Ol/(vz - 01 and h = Iv, - Ol/jO - v3J with 0 being the intersection of the two diagonals of Q. If vi is not v, then LYI , LYZ,a~3 are simply zero. Thus, p and y are zero by (1.2) and hence, all B-coefficients located inside the Q are all zero. Therefore, V,,, is a C’ piecewise
86
M.-J. hi/Computer
Aided
Geometric
Design
13 (1996) 81&W
cubic polynomial function supported on the union of all quadrilaterals which share a common vertex U. (See Fig. 1 for the support of a vertex spline at IJ.) Similarly, we can show that V, is a C’ piecewise cubic polynomial function supported on the union of two quadrilaterals which share a common edge e. (See Fig. 1 for the support of an edge spline around e.) It follows from the cardinal interpolation conditions (1.1)) (1.2), (II.l), and (11.2) that those Vy,y’s and VP’s are linearly independent. The total number of those VU,y’s and VP’s is 3V + E, where V and E denote the number of vertices and edges of 0. Thus, dim($($))
3 3V+
E.
On the other hand, let us define linear functionals A,.,f := DYf( v), y E { (O,O), ( 1, O), e E 1. For any s E S:(e) satisfying &,s = 0 (0, I)}, o E V and A,f := D,(‘s”f(v,,l), for all v E V and y E {(O,O), (l,O), (0, I)} and h,s = 0 for all e E &, we can show as in the proof of the existence of Vi and V, that s = 0. Thus, dim(Si ($)) < 3V + E. Therefore, we have established the following Theorem 1. The dimension of Si ($)
is 3V + E with V and E denoting the number of
vertices and edges of 0. There exist 3V vertex splines VU9,E S: ($) V, E S&(e)
Let us define a linear interpolatory L,t(x)
and E edge splines
which form a basis of A’:($).
operator Lf on S: (9)
by
c DYfWV,,y(x)+~D:'*"f(~,r)V,.(x) PEE VEVYE{(O,O).(I,O),(O,l)}
=c
for any function f E C2( D). Clearly, L,, = p for all p E IP’3since p E S:(e) and V,,,, and K, are a locally supported basis of S: ($). Therefore, we are able to prove the following main result in this note. Theorem 2. The approximation dist(f,$($))
is full, i.e.,
< Cf1014.
Here, Cf is a constant dependent neighboring
order of Si ($)
on
f
and the maximal ratio of the area of any two
triangles of + hnd (01 is the maximum among the diameter of quadrilaterals
of 0.
Proof. We need to use a result in (Bramble and Hilbert, 1971): if a linear functional F defined on Ckt’ (D) satisfies the following two properties ( 1) 1F(f) 1 6 C Cta h’lj D’f II where C is independent of f and h and (2) F(p) = 0 for all p E pk, then there exists a positive constant K independent of f and h such that
IF(f>l
< Khk+‘llDk+‘fll>
where, for any integer 1 3 0,
M.-J. Lai/Cotnpuier Aided Geometric Design 13 (1996) RI-RR
Fig. 3. Franke’s test function based on the quadrangulation in Fig.
87
I
Fix a point x E D. Consider a linear functional F defined by F(f) = f(x) - Lf(x). It is easy to see that our linear functional F satisfies the above two properties. Indeed, F satisfies (2) for k = 3 as mentioned above. Since the locally supported splines Vu,, and V,, are bounded by a constant dependent only on the maximal ratio of two neighboring triangles of $, we conclude that our F satisfies ( 1) for k = 3 with 101 being the maximum among the diameter of quadrilaterals of 0. Thus, we immediately obtain that
If(x) - L,(x)/ 6 K101411D4fll for any x E D, where K is a positive constant
independent
of f, x, and 0.
0
The author has programmed with MATLAB this linear interpolatory operator Lf. For the well-known Franke’s test function f and for the quadrangulation of 46 data in Fig. 1, our Lf gives the plot in Fig. 3. The maximum error of Lf - f is 0.07367. We end the paper with the following remarks. Remarks. 1. Note that any triangulation of the data with the same boundary vertices as that of 0 has the same number of triangles. A triangulation a of the data may be easily obtained by drawn one diagonal of each quadrilateral of 0. Let & be CloughTocher’s triangulation of A. The dimension of S:(H) is 3V + 8, where V and i? denote the number of vertices and edges of & respectively. Since the number of vertices of n is the same as 0, i? = E + NY with NY being the number of quadrilaterals of 0. Clearly, N, = N,/2, where N, denotes the number of triangles of A. Thus, comparing with Clough-Tocher’s scheme, our method reduces the dimension by N,/2, a half of the number of triangles of A. 2. lf a given data cannot be quadrangulated, we should use both triangles and quadrilaterals to partition the data. This results a rectilinear partition. (See (Schumaker, 1984).) Then we subdivide each quadrilateral by its two diagonals and subdivide each triangle into three subtriangles at the center of the triangle. Now we have a refined triangulation. Our method and Clough-Tocher’s method together make possible construct C’ piecewise cubic polynomial function interpolating the data and the interpolatory scheme
88
M.-J. Lai/Computer Aided Geometric Design 13 (1996) 81-88
will have full approximation order. We should use as many quadrilaterals as possible since, as in Remark 1, the dimension of the spline space and the number of patches of the surface can be reduced. 3. Since a four-directional mesh is a triangulation $ based on an obvious quadrangulation, the theorems in this paper generalizes the results in (Lai, 1994). Also, a type-2 partition (see (Schumaker, 1984)) is also a triangulation based on a rectangular grid, i.e., a quadrangulation. Our results can be applied to that type-2 triangulation. 4. Our method can also be compared with Si (d) where 6 is such a triangulation that the number of edges attached to each interior vertex is odd. (6 may be a refinement of A.) Then by (Alfeld et al., 1987)) dim (Si (6) ) = 3V+ ,!?+ Vh, where Vh and I? denote the numbers of boundary vertices and edges of d, respectively. Since .!? 2 i? = E + Nq as in Remark 1, the dimension of Si (A) is bigger than that of Sj ($). Thus, our method is more efficient than using .Si (d). 5. It is natural to study S;(e) for d = 5,6,7. The study of St(+) is easy to do. For .S: (A), Alfeld used Clough-Tocher’s scheme twice to construct interpolating spline functions (see ( Alfeld, 1984) ) . Thus, Sg ($) may be a difficult problem. Recently, the dimension of Si( A) with Clough-Tocher’s triangulation 8 of a was determined in (Gao, 1993). It is possible to study the interpolation and approximation properties of Si ( $) . This was done in (Lai and Schumaker, 1995).
References Alfeld, P. ( 1984). A bivariate C* Clough-Tocher scheme, Computer Aided Geometric Design 1, 257-267. Alfeld, P., Piper, B. and Schumaker, L.L. ( 1987), An explicit basis for C’ quartic bivariate splines, SIAM _I. Numer. Anal. 24, 891-91 I. Bramble, J.H. and Hilbert, S.R. (1971). Bounds for a class of linear functionals with applications to Hermite interpolation, Numer. Math. 16, 362-369. Chui, C.K. and Lai, M.-J. (1990), On bivariate super vertex splines, Con&. Approx. 6, 399-419. Ciarlet, P.G. ( 1974). Sur I’Clement de Clough et Tocher, RAIRO Anal. Number. R2, 19-27. de Boor, C. ( 1987). B-form basics, in: Farin, G., ed., Geometric Modeling, SIAM, Philadelphia, PA, 131-148. de Boor, C. and HBllig, K. ( 1983), Approximation order from bivariate C’-cubits: a counterexample, Proc. Amer. Math. Sot. 87, 649-655. de Boor, C. and Hiillig, K. ( 1988). Approximation power of smooth bivariate pp functions, Math. Z. 197, 343-363. Farin, G. ( 1986). Triangular Bernstein-Btzier patches, Computer Aided Geometric Design 3, 83-127. Gao, .I. ( 1993), A scheme of C2 interpolation over triangulations, manuscript. Gmelig Meyling, R.H.J. (1987). Approximation by cubic C’-splines on arbitrary triangulation, Numer. Math. 5 I. 65-85. Gmelig Meyling, R.H.J. and Pfluger, P.R. ( 1990). Smooth interpolation to scattered data by bivariate piecewise polynomials of odd degree, Computer Aided Geometric Design 7, 439-458. Grandine, T.A. ( 1987), An iterative method for computing multivariate C’ piecewise polynomial interpolants, Computer Aided Geometric Design 4, 307-320. Lai, M.J. ( 1994). Approximation order from bivariate C’-cubits on a four-directional mesh is full, Computer Aided Geometric Design 11, 215-223. Lai. M.J. and Schumaker, L.L. ( 1995), Scattered data interpolation using C* supersplines of degree six, SIAM Numer. Anal., to appear. Schumaker, L.L. (1984), Bounds on the dimension of spaces of multivariate piecewise polynomials, Rocky Mountain J. Math. 14. 251-264.