Asymptotically optimal disposition of tangent points for approximation of smooth convex surfaces by polygonal functions

COMPUTER AIDED GEOMETRIC DESIGN ELSEVIER

Computer Aided Geometric Design 14 (1997) 533-546

Asymptotically optimal disposition of tangent points for approximation of smooth convex surfaces by polygonal functions A.A. Ligun", A.A. Shumeiko a, S.R Radzevitch a, E.D. Goodman b,, a Dneprodzerzhinsk State Technical University, Dneprodzerzhinsk, Ukraine b Case Center for Computer-Aided Engineering and Manufacturing, Michigan State University, East Lansing, M1 48824, USA Received October 1995; revised September 1996

Abstract

Many applications require that smooth convex surfaces be approximated by polygonal functions. While there exist methods for determining optimal tangent points for definition of exterior facetings of such surfaces, they are not practical to compute for many real-world problems. Asymptotically optimal tangent points are nearly as useful, and using the result presented here, are far more easily computed. © 1997 Elsevier Science B.V.

Keywords: Approximation; Surface polygonization; Exterior tangent facets

A ubiquitous problem in engineering geometry is the description of a part surface by means of a mathematical form appropriate for the uses to be made of it. In this article, we investigate one such problem--namely that of finding the best description (in an asymptotic sense) for a smooth convex surface in terms of a set of tangent planes (exterior facets). While there exist methods for determining optimal tangent points for definition of exterior facetings of such surfaces, they are not practical to compute for many realworld problems. Asymptotically optimal tangent points are nearly as useful, and using the result presented here, are far more easily computed (Charrot and Gregory, 1984; Nielson, 1980). Before the result is derived, it is necessary to define the sense in which it will be asserted to be an asymptotically optimal description of the surface. To that end, some necessary concepts are introduced. * Corresponding author. E-mail: [email protected]. 0167-8396/97/$17.00 © 1997 Elsevier Science B.V. All rights reserved. PII S01 6 7 - 8 3 9 6 ( 9 6 ) 0 0 0 4 4 - 1

A.A. Ligun et. al / ComputerAided GeometricDesign 14 (1997) 533-546

534

Given a convex smooth surface defined by the equation z=z(x,y)

((x,y) e D ) ,

(l)

we'll say that domain D is partitioned into D1, D 2 , . . . , Dn, if

0

Di = D,

(2)

i=l and mes(Vi A Vj) = 0 (i ~ j),

(3)

(here mes(D) is a measure on domain D). By a polygonal function (the two-dimensional analog of a piecewise linear curve) we'll mean a continuous function consisting of "pieces" of planes, or facets. An exact definition of a polygonal function follows: A function fn(X, y) defined on domain D is called polygonal if it is continuous on the set D and there exists a partition of domain D into parts D1, D2, • •., Dn such that for every i there exists a plane Pi such that all points (x, y, fn(x, y) ((x, y) E D) belong to the same plane Pi. The polygonal description of surfaces has wide application in many fields--note, for example, that in machine-building, such a problem appears in the forming of a part by a plane cutter, grinding disk, etc. The problem of optimal machining of a surface and many other problems are reducible to the problem of polygonal approximation of the surface. Under the approximation of convex surface by polygonal functions we can consider two problems----descriptions by an external (consisting of tangent planes) polygonal function (Fig. 1) and by an internal (consisting of interpolational planes) polygonal function (Fig. 2). The main problem is the determination of an optimal disposition of tangent points (in the first case) or interpolation points (in the second case). Now we formulate a definition for the quality of approximation of a surface and an asymptotically optimal disposition of nodes (the tangent points or interpolation points). For clarity of exposition, we'll suppose that the part surface is given in the form of (1) and consider the problem of external milling. Denote this surface by 69. For any point M(x, y, z), let

Be(M) = BE(x,y, z) denote a ball with radius ~ and cutter position described at point M. The union Te(69) of all balls with radius ~ and with centers at points on surface 69 will be called the e-layer of surface ~9. The set of points M(x, y, z) ((x, y) E D) such that M E Tc(69) and z >~ z(x, y) will be called the external c-layer of surface ~9 and will be denoted by T + (69), and those such that z <<.z(x, y) will be called the internal E-layer and will be denoted by TZ(69 ). We'll say that polygonal function Zn(X, y) is C- permissible for surface 69 if

z~(x,y) C T+(69)

V(x,y) e D

(i.e., function Zn(X, y) belongs to T+(~9)).

(4)

A.A. Ligun et. a l / Computer Aided Geometric Design 14 (1997) 533-546

Fig. 1. Approximation of convex surface by external tangent facets.

Fig. 2. Approximate of convex surface by internal interpolatory facets.

535

A.A. Ligun et. al / Computer Aided Geometric Design 14 (1997) 533-546

536

The polygonal function Znoo~(x,y) (where integer nopt = nopt(e)) will be called esurface 69 if it is e-permissible and for any e-permissible function zn(x, y), the inequality n >~ nopt is valid. Therefore the e-optimal polygonal function has the minimal number of patches among all e-permissible functions for surface 69. Now let e ~ 0. The sequence of e-permissible polygonal functions {Zn. (x, y) } (where integer n . = n . (e)) will be called asymptotically optimal if

optimal for

lim

n.(e)

nopt(e)

_ 1.

(5)

This relationship is denoted as n. = nopt(1 + o(1)).

(6)

Let Ni = (xi.y~) be arbitrary points of domain D. The equation of the plane tangent to surface 69 at point Ni is

z(z,y)

(y - (z - z(x , = o (7) ~Y (~,y~) If surface 69 is convex, then the lower tangent surface to all surfaces P / i s a polygonal function. Our problem is to find a distribution of points Ni such that this bounding surface is in the e-layer T + (69), where + 0 z ( x , Y____A)

e*

lim - - = 1,

n--+ ~o e

(8)

or, similarly, e* = e ( l + o(1)).

(9)

An exact optimal solution to such problems reduces to a difficult minimax problem and can be found only in easy cases--for instance, if the surface is a sphere, a paraboloid of revolution or a parabolic cylinder. For real-world problems, when surfaces are composite and the number of tangent points (or interpolation points) is sufficiently large, the exact solution of such problems is practically impossible. In such a case, with a requirement for a high-precision description of the surface by a polygonal function, asymptotically optimal methods are almost as effective as the optimal methods. Then the realization of asymptotically optimal methods via computer is relatively simple and convenient. That is a key justification for using asymptotic methods in the solution of engineering problems. This problem may be decomposed into three sub-problems: The first is to build a constructive algorithm for finding the tangent points, the second is to show that the lower bounding surface of the tangent planes at these points lies in the e-layer, and the third is to show that any other polygonal e-permissible function for the surface consists of no fewer planes--i.e., that there is no other polygonal function which can yield a better result. For solution of engineering problems, it is enough to solve the first of the problems stated. The second problem must be solved to prove the assumption of the first problem.

A.A. Ligun et. al / ComputerAided GeometricDesign 14 (1997) 533-546

537

The third (which is the most difficult from a mathematical viewpoint) is necessary only for user confidence in the fact that the given algorithm could be simplified, or adapted to the user's problems, but it cannot be improved upon in the sense of c-optimality. The results obtained in this article are essentially based on the results obtained in references (Ligun and Shumeiko, 1984, 1989), and are their one-dimensional analogs. Let F(t) = (x(t), y(t), z(t)) be a parametrically defined curve and A.~ = { 0 = to < tl < . - - < t,~ = T}

(10)

be an arbitrary partition of the domain of the parameter. Let 7(t, A,~) be a function that coincides over each interval [ti, ti+l] ( / = 0~ 1 , . . . , m 1) with the tangent line to F at some point ti E (ti,tt+l) . If 7(t, Am) C T~(F), then 7(t, Am) will be called e-permissible. The function 7*(t, Amo) (where integer mo = too(E)) will be called e-optimal for curve F if it is e-permissible and for any e-permissible function 7(t, Am) the inequality m >1 too(e) is valid. The sequence of E-permissible functions {'~*(t, Am.)} (where integer m . = m.(e)) will be called asymptotically optimal if lim m . ( e ) _ 1. ~ o no(e)

(11)

For parametric function F(t)= (x(t), y(t), z(t)) such that x, y, z are twice continuously differentiable on the interval [0, T] we shall introduce the following function:

4~(F, t) = ~(t) = ~ IF'(t) x F"(t)] [r'(t)l '

(12)

where ~ × b is the vector product of vectors ~ and b, and ]~[ is the length of vector ~. In coordinate form, the function ~5 can be written as:

~(~)

~ ( z ' u " - x"u') 2 + (u'z" - u"z') 2 + (z'z" - z"x')2 (x,)2 + (y,)2 + (z,)2

(13)

Theorem A. Let functions x(t), y(t) and z(t) determining the curve F(t) be three times continuously differentiable on [0, T] and let c~ E (0, 1/2). Then for any ~ > O, the number of asymptotically optimal nodes of function {~/*(t, Am)} is equal to

m=m(~)=

[i ~

1

e(t) dt + l ,

(14)

o

where [-] denotes the integer part and where nodes Am = {t i}i=0 m are defined from equations t~

f 0

T

f( 0

,m/

Then the piecewise linear curve passing through tangent points Fi+l/2:

538

A.A. Ligun et. al / Computer Aided Geometric Design 14 (1997) 533-546

f~ (~(t) + m-~) dt

.

.

.

. .

.

. .

.

. . . . .

. .

z, / 1/ i

.

. . . . . . .

i j

. . . .

0

tl

~]

,l :

t2

tin-1

tm= T

Fig. 3. Example of asymptoticallyoptimal node placementaccording to Eq. (15). -p(t, ~,.) =

(z'i+l/2(t - ti+l/2) + X~+l/2; y~+,/2(t - t,+,n) + y , + , n ; z~+~/~(t - t~+~n)

+

z~+ln)

will be asymptotically optimal f o r curve F (when ~ ~ 0), where we denote ti+l/2 -(ti+l + t i ) / 2 and

x,+~/2 = x(t~+~/2), d + , / 2 = ~'(ti+~n) and the same is true f o r y(t) and z(t).

The algorithm for asymptotically optimal choice of nodes (15) is illustrated in Fig. 3. Before we begin to describe the algorithm for building the optimal polygonal function, we'll develop some auxiliary constructions which we will require. At each point (xi,j, yi,j, z (x~,j, y~,j ) ), the two (orthogonal) principal directions of curvature are defined, which at our point lie along tangent lines to the curves defined by the Cauchy problem: ( E B - F A ) ( dx) 2 + ( E C - G A ) dx dy + ( F C - G B ) ( dy) 2 = O; Y(Xid) = Yi,j

(16)

where r A= v/l+pZ+q2; F= l+p2;

B-

E=pq;

s ~/l+p2+q2;

C-

1 v/l+p2+qZ;

G=l+q2;

and P=

Oz(x,y). Ox '

q-

Oz(x,y). i~y '

r--

O2z(z'Y)" i~2~ ,

1 -- O2z(x' y) 02 x ;

s--

i~2z(x' y) OxOy

A.A. Ligun et. al / ComputerAided Geometric"Design 14 (1997) 533-546

539

There are two equations of first order: t

X

Y -----f i ( , Y ) ,

"

(17)

y(~c~,j) = y~,~

or, equivalently, x ' = g i ( x , y);

(18)

x(Yi,j ) : Xi, j. where for i = l, 2

fi= 2(FC-

'

GB)

(GA-

+ (-I)iv/(EC

EC

AG): - 4 ( F C - G B ) ( E B - F A ) )

and, similarly,

9' = 2 ( E B 1_ FA) (GA

EC

+ ( - 1 ) ' v / ( J ~ C - - AG) 2---- 4(FC - G B ) ( E B - F A ) ) . In the same way, through every point (xi,j, Yi,j, Z(Xi,j, Yi,j)) there will pass two orthogonal curves, at every point of which the tangent lines are the principal directions of curvature. Denote them by x=t; y = qo(t);

(19)

z=z(t,~(t)) and x = ~,j(t);

y = t;

(20)

Although we noted that point (xo,o,Yo,o, z(x0,o, Yo,o)) was arbitrary, it is better to choose the first point such that (Xo,o, Y0,o) is near the center of gravity for domain D. Let x=t; Y = ~o,o(t);

(21)

z = z(t,

and

{

5 =

y = t; = z(¢o,o(t),

(22) t).

A.A. Ligun et. al/ Computer Aided Geometric Design 14 (1997) 533-546

540

Now, thanks to the algorithm of Theorem A, we construct the collection of points t~ using curve F(t), with error 2~; i.e., we define the number n=

~ ( F , t ) at

+l,

(23)

0

and the nodes from the condition: t~

T

(i = 0 , 1 , 2 , . . . , n ) .

n

o

(24)

o

And from curve Fo(t) we define the collection of points tj,o with error ~ V'3/2; the number of partitions will be:

ltO "-~- ~

1

0

~(Fo, t) dt

+ l,

(25)

and the nodes, from the condition, are: tj,o

To

f(~(F°'t)+n°~)dt=J--f no (~(F°'t)+n°~)dt 0

(i = 0, 1 , 2 , . . . , n 0 ) .

(26)

0

Then, in both cases, we shall choose points corresponding to the values of parameter t~ (i = 0, 1 , . . . , n) and tj,o (j = O, 1 , . . . , no), that lie in domain D (Fig. 4). By F~(t) we shall denote the direction of one of the principal curvatures and passing through point Fi = F(ti) orthogonally to F(t) (Fig. 5), and on each of the given curves define the set tj,i with error c v/-3/2 (from Theorem A). Supplement the set of points obtained, Fi,j = (xi,j, Y i j , z ( x i , j , Yi,j))

(27)

by points Fi+l/2,j+l/2 = (Xi+l/2,j+l/2, Yi+l/2,j+l/2, Z(Xi+l/2,j+l/2, Yi+l/2,jq-l/2)),

(28)

where

xi+l/z,j+l/2 = ¼(x~,j + X~+l,j + xi,j+l + xi+l,j+l),

(29)

Yi+l/2,j+l/2 : !4\(Yi ,9" ~- Yi+l,j "~- Yi,j+l "~ Yi+l,j+l).

(30)

Then the polygonal surface formed from the tangent planes to z = z(x, y) at the points F~,j and Fi+l/2,j+l/2, i.e., (Figs. 6--7)

y) v~,j (V - Y~,j) - (z - z(x~,j, Ui,j)) = 0 Oz(X,ox y) r, ,~ (x - xi,j) + Oz(x, 0-----~ and

Oz(x, y)

Ox

r,+,/2~+,/ (x _

Xi+l/2'j+l/2 ) ~- Oz(x, y)

0----7-

F,+,/:

j+,/2 (y

Yi+l/Z,j+l/2)

(31)

A.A. Ligun et. al / Computer A i d e d Geometric Design 14 (1997) 533-546 /

f

541

ro(t)

(z2,o, y2,o) i1

/

[ (x,,o, yJ,o)

v(t)

II

~ .~r~"i " l j

_.~

..........

j.J"~

,Y0)

"J'J

~ ....

(X0,2, y0,2)

(X0,1, Y03)

1

Fig. 4. Selection of points along two curves orthogonal at point (xo, Vo) in domain D, according to Theorem A.

J

4 : /"2,l

I~[

/"

ff 1"20'

::

i /

"_r~,l

" r ,0

F1,2

?

• ...................... • ...........................

• ......................... "i

~/'2,2

Fo,o

....., .................. iv0,,

i

:

ii

! 1~0,2

~F-I,2

2

Fig. 5. Application of Theorem A along additional curves Fi (t) orthogonal to F.

--

(Z -- Z(Xi_bll2,jq_l/2,

ffi+l/2,j+l/2))

~--- 0

(32)

will be asymptotically optimal from above. To prove this fact, we first show that the given polygonal function lies within the e-layer of surface 69. As regards the possibility

542

A.A. Ligun er al / Computer Aided Geometric Design 14 (1997) 533-546 /

./

/

"

:"

/

.: Fi+I,j+ I

// Fi,j+l

f

::

i ii

/ f

•

f • :Fi+l/Z,j+l/2

•

!

i Fi+2 j+l

Fi+3/2,j+l/2

J i~

~i

,

~

* Fi+2,j F i + l ,j

:!r~,j Fi+3/2,j_l/2

"

l'i+l/2,j_l/2

•

/

Fi,j-1

,

Fi+l,j-I

' Fi+2,j-- 1

Fig. 6. Supplementing with intermediate points xi+l/2d+l/2 from Eqs. (29), (30).

/

\ M1 ,o

/

x~ ~t~''l"

M-iI2,1/

M-l,O

M-II2,~12

k)~

Fig. 7. Facets produced by tangent planes at the points shown on Fig. 6.

A.A. Ligunet. al / ComputerAidedGeometricDesign14 (1997)533-546

543

that it lies in the upper part of the e-layer, then it would follow from the condition of convexity upward of surface 69 and the construction method of the polygonal function (by the definition of convexity, that the surface will be convex upward on the domain D, if for any point (x, y) E D the corresponding point on any tangent plane lies higher than the corresponding point on the surface). Let us take point xi, j : 0 , Yi,j = 0. Then the tangent plane to the surface z = z(x, y) at the point (0, 0, zi,j) (here zi,j = z(0, 0)) will have equation (see (7))

Oz(x,oxy) (o,o)x + Oz(x,o___if__ y) (o,o)y - (z - zi,j) = O.

(33)

Then the distance from point M(x, y, z(x, y)) to this plane will be equal to

d(M) = I~z(x' Y)/~xl(°'°)x + ~z(z, Y)/~YI(o,o)Y - (z - zio) I V/(Oz(x, y)/oxl(o o)) 2 + (Oz(x, y)/oyl(o,o))2 + 1

(34)

The points of the surface which are removed from the tangent plane by a distance of no more than e will be described by inequality

d(M) <~e.

(35)

By means of Taylor's formula for z = z(x, y) in the neighborhood of point (0,0)

Oz(x, y) Oz(~, ~')1 z(~:,y)-- z~j + D---T- (oo)X + 0---7-- (oo)y + ~. \

~x 2

Ox~y

(o,o)x 2 +

(o,o)xy + ~ 2 ~

+ O((x + y)3),

(36)

we have

d(M) =

zi,j

+

Oz(x, y) Ox (o,o)x

-

__ Zi,J

Oz(x, - - y) (0,0)y

OZ(X,~xy) (o,o)X

I (02z(z'Y) I

2! ~,

+

~x2

+ O((x + y)3)

x2+ (o,o)

Oy

Oz(X,~yY) (o,o)Y 202z(x,y)

i~xi~y

(o,o)XY

O:z(x,y) i)2Y

y: )

(o,o)

1

Vf(~z(x,y)/~xl(o,o)) 2 + (~z(~, ~)/~vl(o,o)) 2 + 1

_-1~ lAx 2 + 2Bxy + Cy2l <~e+ O((z + y)3).

(37)

Therefore the set of points (x, y) such that values z(x, y) are removed from the tangent plane at point (0,0) by a distance of no more than e is described by the inequality

21lAx 2 + 2Bxy + Cy21 <~e + O((x + y)3).

(38)

544

A.A. Ligun et. al / Computer Aided Geometric Design 14 (1997) 533-546

It is not difficult to see that this inequality defines the interior of an ellipse with semi-axes 1/~ and 1/2x/2-~z, where Ai, the eigenvalues of the quadratic form A x 2 + 2 B x y + C y 2, are the solutions of the equation A-A

B

B

C-A

= O,

(39)

and Ai = C + A 4 - v / 4 B 2 + ( C - A ) 2 (40) 2 One can see from the construction that the given polygonal function is a covering of surface z = z ( x , y) by hexagons (Fig. 7). To prove the first part of the theorem, it suffices to show that the top of each hexagon lies at least in some ellipse ½ [A(x~,j, y i , j ) ( x - x~,j) 2 + 2 B ( x i , j , y i , j ) ( x - x~,j)(y - yi,j) + C(xi,j,yi,j)(y

- yi,j)2[ <~ e + O(((x - xi,j) + (y - y~,j))3),

(41)

where M i , j = (x~,j, y~,j, z ( x i , j , yi,j)) is a tangent point of the plane. Consequently, we show that this point lies within the e (1 + o( l))-layer of the surface. For example consider the vertex that is formed by the intersection of tangent planes at points F~,y, Fi+l/2,j+l/2 and Fi,j+l. This point is defined by the system:

{

p~,j(x

- x~,j) + q,,j(v

- w,j)

+ z - z~,~) = 0;

P i + l / 2 , j + l / 2 ( x -- X i + l / 2 , j + l / 2 ) + q i + l / Z , j + l / 2 ( Y -- Yi+l/2,j+l/2)

(42)

+ Z -- Zi+l/Z,j+l/2) : O;

pi,j+l (x - xi,j+l ) + qi,j+l (Y - Yi,j+l) + z - zi,j+l) = O.

where Pi,j-

~z(x, y) I

I(x~,j,yi,~)'

(43)

~ z ( x , y) (~, ,~,u,,j ),

(44)

~

qi,j : ~

etc. Solving this system and using a decomposition by Taylor's formula in the neighborhood of point ['i,j, we can to verify that at this point l-lA(xi,J,2

y i , j ) ( x - xi,j) 2 + 2 B ( x i , j , y , , j ) ( x - x i , j ) ( y - Yi,j) + C(xi,j,yi,j)(y-

Yi,3/)21I = e(1 + o(1))

(45)

which demonstrates that the polygonal function in the neighborhood of this point lies in the e-layer of the surface. The proof that the given method yields asymptotically optimal results is rather difficult. The idea of this proof, without detailed analysis, is based on the following reasoning: first, we approximate the given surface by elliptical paraboloids on the domains for which

A.A. Ligun et. al / Computer Aided Geometric Design 14 (1997) 533-546

545

diameters are on the order of e 3/2. Then each of these paraboloids is approximated by a best polygonal function. In that case, the projections of hexagons of the polygonal function on the plane X O Y (which don't touch the boundaries) will be congruent, and the problem will be reduced to the packing of the plane by identical hexagons; the error of approximation on each of these patches will be e. Then it may be seen that their number asymptotically equals the number of hexagons yielded by the algorithm under consideration. It then follows that for surfaces consisting of pieces of elliptical paraboloids, there is no better polygonal function. It remains to be noticed that any sufficiently smooth convex upward surface may be approximated by a surface composed of elliptical paraboloids (using the same number) with an accuracy of e 3/2, and this surface may be used as an intermediate approximation. Schematically, note how we prove the third part of the theorem, i.e., that we cannot construct a polygonal function with a smaller number of planes. Divide the domain D into N parts Di with approximately equal areas, where N=

[--~e2]+l

(46)

and on every ith piece choose a point (xi, Yi) approximately near the center of gravity of this part and calculate Ai, Bi and Ci at these points and ai = I/2V/~l and bi = 1 / 2 v ~ 2 (where A,- are the eigenvalues of quadratic form Aix 2 + 2Bixy + Ciy 2) or the solutions of the characteristic equation

Ai - A

Bi = 0.

Bi

(47)

Ci -

Now, for this piece, substitute for z = z(x, y) an elliptic paraboloid z -- (x - Yi)~ + (y - Yi)~ (ai) 2 (bi) 2

(48)

or z-

-

(bi) 2

+

(ai) 2

'

(49)

where xi, Y'i is some point from Di. By doing this, the error of description will be

(')

O ~'~

= O(E 2)

(50)

and the proposed algorithm will yield an asymptotically optimal covering of this paraboloid by a polygonal function. It is possible to verify that for every Di the polygonal function will consist of identical hexagons in which the most distant points (all vertices of hexagons) are removed from the paraboloid by e(1 + o(1)). So, in that manner, Di will be densely packed by these hexagons, so Di cannot be filled by a smaller number of hexagons. As the description by paraboloids is accurate to O(e 2) and by the polygonal function to O(e), that proves the third part of the theorem--i.e., that the given algorithm is asymptotically optimal.

546

A.A. Ligun et. al / Computer Aided Geometric Design 14 (1997) 533-546

We have proved the analog of Theorem A for the case of using interpolated piecewise curves. By use of this theorem (analogously to the use of Theorem A), we can construct an algorithm for definition of an optimal polygonal function that consists of simplexes which approximate a convex-up surface from within, and a convex-down (concave) surface from outside. We will discuss this problem in detail in a future article.

References Charrot, E and Gregory, T. (1984), A pentagonal surface (tP) patch for computer-aided geometric design, Computer-Aided Geometric Design 1. De Boor, C. (1978), A Practical Guide to Splines, Springer, New York. Ligun, A. and Shumeiko A. (1984), On optimal choosing of nodes for approximation of functions by interpolation splines, J. of Calc. Math. and Math. Ph. (USSR) 9, 1283-1293 (in Russian). Ligun, A. and Shumeiko, A. (1989), Asymptotically optimal algorithm by piecewise linear interpolator in problem the tooling by milling cutter, Izv. VUZov. Machine-building (USSR) 7 (in Russian). Nielson, G. (1980), Minimum norm interpolation in triangles, SIAM J. Numer. Anal. 17, 46--64. Schumaker, L.L. (1981), Spline Functions, Basic theory, Wiley, New York.