Generalized barycentric coordinates and ... - Semantic Scholar

Received 16/09/13

Generalized barycentric coordinates and approximations of convex functions on arbitrary convex polytopes Allal Guessab Laboratoire de Math´ematiques et de leurs Applications, UMR CNRS 4152, Universit´e de Pau et des Pays de l’Adour, 64000 Pau, France [email protected]

Abstract In this paper, we study the error in the approximation of a convex function obtained via a one-parameter family of approximation schemes, which we refer to as barycentric approximation schemes. For a given finite set of pairwise distinct points X n := {xi }ni=0 in Rd , the barycentric approximation of a convex function f is of the form: n X B[f ](x) = λi (x)f (xi ), i=0

{λi }ni=0

is a set of barycentric coordinates with respect to the point where set X n . The main contents of this paper is two-fold. The first goal is to derive sharp upper and lower bounds on all barycentric coordinates over the convex polytope conv(X n ). The second objective of the paper is to exploit the convexity assumption heavily and establish a number of upper and lower pointwise bounds on the approximation error for approximating arbitrary convex functions. These bounds depend solely on computable quantities related to the data values of the function, the largest and smallest barycentric coordinates. For convex twice continuously differentiable functions, we derive an optimal error estimate. We show that the Delaunay triangulation gives access to efficient algorithms for computing optimal barycentric approximation. Finally, numerical examples are used to show the success the method. Keywords: Barycentric approximation schemes–Barycentric coordinates–Convex functions– Delaunay triangulation–Optimal approximation–Nonexpansive approximation– Lower bounds– Upper Bounds–Polytopes Preprint submitted to Computers & Mathematics with Applications

September 16, 2013

1. Introduction, notations and preliminary results The basic idea behind barycentric approximation schemes is the following: Consider a given finite set of pairwise distinct points X n = {xi }ni=0 in P ⊂ Rd , with P = conv(X n ) denoting the convex hull of the point set X n . We do not make any assumptions about the distribution of the points in P except that they must not lie entirely in the same straight line. We are interested in approximating an unknown scalar-valued convex continuous function f : P → R from given function values f (x0 ), . . . , f (xn ) sampled at X n . In order to obtain a simple and stable global approximation on P , we may consider a weighted average of the function values at data points of the following manner: n X B[f ](x) = λi (x)f (xi ), (1) i=0

or, equivalently, a convex combination of the data values f (x0 ), . . . , f (xn ). This means that the system of functions λ := {λi }ni=0 forms a partition of unity, that is: n X

λi (x) ≥ 0, ∀x ∈ P, i = 0, . . . , n,

(2)

λi (x) = 1, ∀x ∈ P.

(3)

i=0

In addition, we shall also require the set of functions λ to satisfy the firstorder consistency conditions: x =

n X

λi (x)xi , (∀x ∈ P ).

(4)

i=0

We will call any set of functions λi : P → R, i = 0, . . . , n, barycentric coordinates if they satisfy collectively, for all x ∈ P, the three properties (2), (3) and (4). In view of this property we shall refer to the approximation schemes B as barycentric approximation (schemes). The barycentric approximation has many nice properties and can be very useful in some situations. For instance, conditions (3) and (4) together guarantee that affine functions are exactly reproduced by the approximation scheme B. For second-order partial differential equations, approximants that possess constant (3) and linear precision 2

(4) are sufficient for convergence in a Galerkin method [37]. The positivity condition (2) together with the partition of unity property (3), allow us to interpret the barycentric approximation as a weighted average of the data values. This viewpoint is common in computer graphics and geometric modeling, e.g. in B´ezier and B-Spline techniques [8, 33]. A very useful property of (2) and (3) is that if x is given such that 0 ≤ λi (x) ≤ 1 for all i, then x is inside the polytope P . Equations (2) and (3) also ensure that approximation operator B is bounded between the data’s minimum and maximum values: mini f (xi ) ≤ B[f ](x) ≤ maxi f (xi ). This guarantees a stability property of the operator B and it is also a statement of the discrete maximum principle and a requirement for the numerical discretization of the diffusion equation. Note that B is in general non-interpolatory on the point set X n . However, at the vertices of the polytope P the function λi associated to the vertex v i has to equal 1 at v i and 0 at all other points in X n \ {v i } , i.e., λi (v j ) = δij ( δ is the Kronecker delta). This property is the foundation of using barycentric coordinates for interpolation purposes. Indeed, it immediately implies vertex interpolation for any continuous function: B[f ](v i ) = f (v i ). In this regard, we should mention that vertex interpolation is satisfied automatically by any positive linear operator, see [20]. Clearly, if f is nonnegative, then so is B[f ]. Hence, the weighted average approximation B preserves the nonnegativity property and then it belongs to the class of positive operators. We will also see that all functions λi associated to points in the interior of X n vanish at the boundary of P , such property facilitates the imposition of linear Dirichlet boundary conditions in a Galerkin method, [4, 29]. It may happen that we know beforehand that the function f to be approximated is convex. To our knowledge, however, subject to the usual convexity, no careful analysis has been designed to be able to provide the correct error estimates for these barycentric approximation schemes. There is some literature on estimation of functions restricted by convexity. For instance, by allowing some smoothness on f to be approximated, say C 2 (P ), it played a crucial role in our research on the determination of the ‘best’ (or ‘optimal’) cubature formulae, see [14, 16, 17, 18, 19, 21], where the latter problem has been extensively reviewed both from the theoretical study as well as the numerical point of view. It is to note that, in the context of one dimensional case, many works have been devoted to the obtention of upper bounds for the discrete Jensen’s inequality, see e.g. [6, 28, 36] and references therein. Therefore, in this paper, we shall exploit the convexity assumption heavily, and establish a number of sharp upper and lower pointwise bounds on the 3

error inherent in barycentric approximation for approximating arbitrary convex functions. The remainder of this paper, which we have divided into six parts, is organized as follows. In Section 2, we first fix notation and provide a very brief survey and some relevant refinements on barycentric coordinates on polytopes. We also give basic existence results. We then derive in Section 3 sharp upper and lower bounds on all barycentric coordinates over convex polytopes. In section 4, we present some of the essential properties of barycentric approximations schemes, and their corresponding approximation errors. In particular, it is shown that a weak Kronecker-delta property at the boundary holds. This section also provides a characterization of the nonnegativity property of the approximation error on the class of convex functions. We also explain in this section why the framework given in this paper was limited to the case of polytope domains. In Section 5, we establish a number of upper and lower pointwise bounds on the approximation error for approximating convex functions. Our main result is Theorem 5.2, which shows how to control, for any convex function, the pointwise error at each x ∈ P committed by the approximation B in terms of computable quantities related to the data, the largest and smallest barycentric coordinates of the point x. In this section, sharp bounds are also derived through the use of the bounds on all barycentric coordinates obtained in Section 3. In Section 6, for convex twice continuously differentiable functions, we derive an optimal error estimate. We show that the Delaunay triangulation gives access to efficient algorithms for computing optimal barycentric approximation. Some numerical results presented in Section 7 illustrate the success of piecewise linear interpolation using Delaunay triangulation. 2. Generalized barycentric coordinates on polytopes and existence results In this section we will give the notation we will use, and extend the notion of barycentric coordinates to arbitrary polytopes. We will also show their existence under very general conditions. Consider a set of pairwise distinct points X n = {xi }ni=0 in Rd . To simplify the remaining presentation, let us introduce some notations: Throughout this paper d and n denote some fixed positive integers. The convex hull, P ⊂ Rd is the smallest (inclusionwise) convex set that contains all the point 4

set X n . The convex hull of a finite set of points in some Rd is called a (convex) polytope. For simplicity, we drop the qualifier “convex” throughout, since most polytopes under consideration will be convex. A point x of a polytope P is called an extreme point if there are no two points y 6= x and z 6= x in P such that x is a convex combination of y and z. For simplices, barycentric coordinates are a very common tool in many computations. Basically, they are defined as follows: let X d = {v 0 , . . . , v d } be any linearly independent set of d + 1 points in Rd , the simplex T with the set of vertices X d is the convex hull of X d , (e.g., a triangle in 2D or a tetrahedron in 3D). Let Ai (x) be the signed volume (or area) of the simplex created with the vertex v i replaced by x. Then the barycentric coordinate functions {λi , i = 0, . . . , d} of the simplex T with respect to its vertices are uniquely defined by λi (x) =

Ai (x) , vol(T )

(5)

where vol(T ) will mean the volume measure of T . It is easily Pd seen that each point x of T has a (unique) representation, that is x = i = 0 λi (x)v i and the barycentric coordinates λ0 , . . . , λd are nonnegative affine functions on T . The uniqueness of this representation allows the weights λi (x) to be interpreted as an alternative set of coordinates for point x, the so-called barycentric coordinates. Note that a simplex is a special polytope given as the convex hull of d + 1 vertices, each pair of which is joined by an edge. For a convex polytope P ⊂ Rd we will use ‘generalized’ barycentric coordinates (They are often called generalized barycentric coordinates to distinguish them from the original barycentric coordinates, which were only defined with respect to simplices.) While barycentric coordinates are unique for simplices, there are many possible solutions for polygons with more sides. In recent years, a wide range of generalized barycentric coordinates has been suggested. Two important examples of generalized barycentric coordinate functions are Wachspress’ coordinates described in [40] and the mean value coordinates of [8]. Both families produce smooth functions on the polytope P . While being initially described for the 2D case, both have generalizations to 3D. In recent years, the research on barycentric coordinates has been intensified and led to a general theory and extensions to higher dimensions [25, 24, 26, 11, 12]. Barycentric coordinates are natural coordinates 5

for meshes and have many applications including parameterization [5, 32], free form deformations [34, 24], finite elements [38]. For additional references as well as more details on this topic, refer to the complete web page (http://www.inf.usi.ch/hormann/barycentric/). The first result on the existence of barycentric coordinates for more general types of polytopes was due to Kalman [27, Theorem 2] (1961). The next lemma is due essentially to Kalman [27]. Our statements are stronger than the ones provided in [27], but the proof proceeds along the same lines as the proof of theorem 2 in [27], so we omit it. Lemma 2.1. Let X n = {x0 , . . . , xn } be a set of finite pairwise distinct points of Rd and let the polytope P = conv(X n ). Let x∗ be a given point in P with n n X X ∗ ∗ x = λi xi and λ∗i = 1 , λ∗i ≥ 0. (6) i=0

i=0

Then there exist nonnegative real-valued continuous functions λ0 , λ1 , . . . , λn defined on P such that x =

n X

λi (x)xi

and

i=0

n X

λi (x) = 1

(7)

i=0

for each x ∈ P and λi (x∗ ) = λ∗i . Moreover, we can choose the barycentric coordinates {λ0 , . . . , λn } in such a way that one of them is convex or concave. In addition to the assumptions of the previous Lemma, we now assume that there are no other points in X n than the vertices of the polytope P . This situation is of interest when P is a small convex polytope obtained by subdivisions of the original domain. The most important fact of the next result is that the barycentric coordinates are linearly independent. In fact, we have the following: Lemma 2.2. Let P be a polytope in Rd , {v 0 , v 1 , . . . , v m } its vertices. Then there are linearly independent nonnegative real continuous functions on P , λ = {λ0 , . . . , λm } defined on P such that x =

m X

λi (x)v i

and

i=0

m X i=0

6

λi (x) = 1

(8)

for each x ∈ P. Proof. The existence of a set of continuous barycentric coordinates λ is assured by Lemma 2.1. So it remains to show that the set consisting of the functions λ is linearly independent. Linear precision (8) shows in particular that each x may be represented as a convex combination of v 0 , v 1 , . . . , v m . Since each v i is an extreme point of P , we conclude by substituting x = v i in (8) that λi (v j ) = δij . Hence, the functions λ0 , . . . , λm satisfy the delta function property. Now, it is easy to see that this property implies that the set of functions λ must be linearly independent. Barycentric coordinates provide a basis for linear finite elements on simplices, and generalized barycentric coordinates naturally produce a suitable basis for linear finite elements on general polytopes. Here, the underlying principle is that one triangulates the polytope into simplices and then use the standard barycentric coordinate functions of these simplices. Interpolation properties of this scheme are well known from the standard analysis of the finite element method over simplicial meshes, but, as we will see in the next section, this construction will serve as an important point of comparison with any set of barycentric coordinates. Indeed, a most common form of approximation in the standard continuous linear finite element method is piecewise linear (piecewise affine) spline interpolation, see, e.g., Ciarlet [3] among many standard reference books on finite element theory. To describe a general set up, let X n = {xi }ni=0 in Rd be a finite subset of P ⊂ Rd with the property P = conv(X n ). A triangulation T of P with respect to X n is a decomposition of P into d-dimensional simplices such that X n is the set of all their vertices and the intersection of any two simplices consists of a common lower dimensional simplex or is empty. Triangulations of compact convex polytopes exist. Indeed, given any finite set X n of points that do not all lie on a hyperplane, Chen and Xu [2, p. 301] describe a lifting-andprojection procedure which results in a triangulation of the convex hull of X n . Let S1 , . . . Sl be the simplices of T and let Ni be the set of all integers j such that xi is a vertex of Sj . If x ∈ Sj and j ∈ Ni , then we denote by λij (x) the barycentric coordinate of x with respect to xi for the simplex Sj . It is easily verified that if x ∈ Sj ∩ Sk , then λij (x) = λik (x) if j, k ∈ Ni and λij (x) = 0 if j ∈ Ni , k 6∈ Ni . Therefore, setting ( λij (x) if x ∈ Sj and j ∈ Ni lif e (x) := 0 otherwise 7

for i = 0, . . . , n, we trivially obtain a set of well-defined barycentric coordinates, which are a generalization of barycentric coordinates when P is a simplex. We list some basic properties of the functions {l0f e , . . . , lnf e }

(9)

of which the following are particularly relevant to us: (1) They are well-defined, piecewise linear and nonnegative real-valued continuous functions. (2) The function lif e has to equal 1 at xi and 0 at all other points in X n \ {xi } , i.e., lif e (v j ) = δij ( δ is the Kronecker delta). 3. Sharp upper and lower bounds on all barycentric coordinates Floater, Hormann and K´os [7] derived sharp upper and lower bounds on all barycentric coordinates on convex polygons in the plane. The main objective of this section is to show that the results of [7] can be extended in a more general way for arbitrary polytopes in higher dimensions. The important property that will be used later is that the polytope P = conv (X n = {x0 , . . . , xn }) can also be characterized by means of a finite number of linear halfspaces.  P = x ∈ Rd : hk (x) = hak , xi + bk ≥ 0, k = 0, . . . , m , (10) here ak ∈ Rd and bk ∈ R, see, e. g., [41]. Let (nv + 1) be the number of vertices of P . We will use the notation V (P ) := {v 0 , v 1 , . . . , v nv } to denote the set of all vertices of P. Throughout this section P1 is the class of all polynomials in d real variables of degree at most 1, also called the class of affine functions on Rd , and λ := {λ0 , λ1 , . . . , λn } an arbitrary but fixed set of barycentric coordinates with respect to X n . The next Lemma provides an upper bound for each barycentric coordinate associated to a point of X n , which is in the interior of P . Here, hj , j = 0, . . . , m, are the hyperplanes defining the polytope P, see (10). 8

Lemma 3.1. Let {λ0 , λ1 , . . . , λn } be any barycentric coordinates system with respect to X n = {x0 , . . . , xn }. Then, for any xi ∈ X n ∩int(P ) the associated barycentric coordinate λi satisfies the bounds: hj (x) ≤ 1, (∀x ∈ P ). 0≤j≤m hj (xi )

0 ≤ λi (x) ≤ min

(11)

Proof. We first establish that any set of barycentric coordinates by definition reproduce affine functions. Indeed, let l ∈ P1 , then there exist a vector a ∈ Rd and a real number b such that l can be expressed as follows: l(x) = ha, xi + b, where h., .i denotes the usual inner product in Rd . Then, assuming (3) and (4), it follows that n X

l(xi )λi (x) =

i=0

n X

(ha, xi i + b) λi (x)

(12)

i=0

* =

a,

n X

+ λi (x)xi

+b

i=0

n X

λi (x)

i=0

= ha, xi + b = l(x), (∀x ∈ P ). Therefore, for all l ∈ P1 , B[l] = l.

(13)

This shows that B reproduces exactly all affine functions. Now take any point xi in X n ∩ int(P ), and let λi its associated barycentric coordinate. We next show that λi satisfies the following bounds: λi ≤

hj , (j = 0, . . . , m). hj (xi )

Indeed, making use the reproduction property of affine functions and the

9

nonnegativity of

hj , hj (xi )

we conclude that, for any j = 0, . . . , m, λi (x) = λi (x) ≤

n X

hj (xi ) hj (xi )

λk (x)

k=0

hj (xk ) hj (xi )

hj ](x) hj (xi ) hj (x) = . hj (xi )

= B[

(14) Hence, any barycentric coordinate associated to a point in X n ∩ int(P ) must satisfy the required inequality (11), as claimed. For the barycentric coordinates associated with each vertex, we get more precise bounds. Given a vertex v i of the polytope P. First of all, let us denote by Pi1 and Pi2 the two polytopes defined by the convex hull of X n \ {v i } and P \int(Pi1 ) respectively. Then, it is easily seen that Pi1 and Pi2 are two smaller polytopes of P, and they form a partition of P. As was shown in the example of Figure 1, the set Pi2 , in general, need not be convex. Let us denote by lji , j = 0 . . . , mi the hyperplanes defining the common facets of Pi1 and Pi2 . These hyperplanes with the additional normalization condition lji (v i ) = 1 play a predominant role in determining lower bounds on any barycentric coordinates. Let us make now some comments on the hyperplanes lji . Remark 3.2. The hyperplanes lji enjoy many truly remarkable properties: (i) Since lji , j = 0 . . . , mi , are supporting hyperplanes of Pi1 , then the polytope Pi1 is entirely contained in one of the two closed halfspaces determined by lji . Hence, it follows, in particular, from the fact that lji (v i ) = 1, lji (xk ) ≤ 0, k = 0, . . . , n, k 6= i. (15) (ii) For each i = 0, . . . , nv , the function associated with vertex v i defined by  ˜li := max 0, li (16) j j=0...,mi

10

vi b

Pi2 b

b

xi

b b

b

b

b b

b b

b b b

b b

b

b b

b

b

b

b b b

b

b

b b

b

b b

b

b

b

b b

b

b b

b b

b b

b b

b b b

b b

b

b b

b b

b b

b

b

b b

Pi1 b

b

Figure 1: Partition of the polytope P into two polytopes.

is obviously a nonnegative convex function on P, that it takes the value 1 at the vertex v i . We note that it also follows from (15) that ˜li (xk ) = 0, k = 0, . . . , n, k 6= i.

(17)

Hence ˜li satisfies the delta-property, and then, in particular, ˜li vanishes at all vertices of Pi1 . Thus, if follows from nonnegativity and convexity of ˜li , that it vanishes entirely on the polytope conv (X n \ {v i }) := Pi1 . After these preparations, we can prove the following result: Lemma 3.3. Under the above notations, any barycentric coordinate λi associated with a vertex v i satisfies for all x ∈ P the bounds: 0 ≤ ˜li (x) = max

j=0...,mi



˜ i (x) = min hj (x) ≤ 1, 0, lji (x) ≤ λi (x) ≤ L j∈∆i hj (v i )

(18)

where ∆i is the set of all indices j such that the hyperplane hj does not vanish at the vertex v i . Proof. Fix a vertex v i of the polytope P, and let λi its associated barycentric coordinate. The upper bound in (18) may be proved in a similar way by following the proof of Lemma 3.1, so we omit it. To prove the lower bound, we first note, as mentioned in Remark 3.2, (ii), that the function ˜li vanishes on Pi1 , then we may assume that x ∈ Pi2 , since otherwise the lower bound is 11

automatically satisfied. Hence, let us pick x any point in Pi2 . Then, since lji is an affine function lji (x) =

n X

λk (x)lji (xk )

k=0

=

λi (x)lji (v i )

+

n X

λk (x)lji (xk )

k=0,k6=i

(19) Hence, in view of (15), lji (x) can be bounded by lji (x) ≤ λi (x).

(20)

This implies the required lower bound and completes the proof of the lemma. Lemma 3.3 extends [7, Proposition 1] to any dimension and as immediate consequences of the last two Lemmas, we have the following general results. Corollary 3.4. Any set of barycentric coordinates satisfies: (i) The barycentric coordinates associated to each point belonging to the interior of the polytope P vanish on facets of P . (ii) The barycentric coordinates associated to the vertices satisfy the delta property. Some of these properties in different forms are known in many contexts. For instance, as a consequence of Corollary 3.4 B´ezier curves pass through the end control points and are tangent to the end control-polygon edges, see [8, 15]. Fix v i a vertex of P, and for each j ∈ {0, . . . mi }, let Qji be the polytope
defined by conv Fij ∪ {v i } , where Fij is the facet defining hyperplane lij . With this notation it is clear that these polytopes contain v i as a vertex and form a partition of Pi2 . The polytope Qji can be triangulated using no new vertices, and such that any simplex in T has v i as a vertex, see [1], see Figure 2. Note that the uniqueness of the barycentric coordinates for simplices implies that the barycentric coordinates associated to the vertex vi with respect to any simplex of T is exactly the equation of the hyperplane lji . 12

vi b

Qji

Pi2 b

b

b b

b

b

xi

b b

b b

b b b

b

b b

b b

b b b

b

b

b b

b b

b b

b

b

b

b

b

b b

b

b

b

b

b b

b b

b b

b

b b

b b

b b b

b b

b b

b

b b

Pi1 b

b

Figure 2: An example of a two-dimensional triangulation of the polytope Pi2 into simplices.

˜ i can be obtained The next result shows that no better bounds than ˜li and L for any set of barycentric coordinates. In view of this fact, we call them sharp lower and upper bounds. Indeed, we will show more precisely that there is a set of barycentric coordinates for which these bounds are themselves ith barycentric coordinates. We say that a function λi : P → R is an i-th barycentric coordinate if it is the i-th member of some set of barycentric coordinates. With the help of Lemma 2.1, we can prove the following: ˜ i associated to any vertex v i are Proposition 3.5. The functions ˜li and L i-th barycentric coordinates. Proof. Given a vertex v i of P. The main ingredient in the proof is the key Lemma 2.1, which says that there exists a set of barycentric coordinates {β0 , β1 , . . . , βn } such that the barycentric coordinate βi associated to v i is a convex function (or a concave function). We will now show that under the convexity assumption of βi , the two functions βi and ˜li are identically equal on P . Corollary 3.4, (ii) implies that βi must satisfy the delta property, that is βi (v j ) = δij , (21) and consequently βi vanishes on the polytope conv (X n \ {v i }) := Pi1 . Hence, βi and ˜li are equal on Pi1 . Let now x be any point of Pi2 = P \ int(Pi1 ), then it belongs to at least one simplex Sji of the triangulation T of the polytope Qji . Since βi satisfies the delta property at the vertices of any simplex of Sji , 13

it follows immediately from Jensen inequality that βi (x) ≤ lji (x), here, we have taken into account the fact that the barycentric coordinate with respect to any simplex of T associated the vertex v i is the hyperplane lji . Now, this immediately implies that βi (x) ≤ ˜li (x). Finally, this together with the opposite inequality proved in Lemma 3.3, satisfied by any set of barycentric coordinates, shows that βi = ˜li on Pi2 . ˜ i = βi . We can show in a completely analogous way that if βi is concave then L Thus we have proved Proposition 3.5. The result of Proposition 3.5 is precisely the result of Floater, Hormann and K´os [7, Proposition 2.4] derived on convex polygons in the plane. The next result shows that the reverse implication of Lemma 3.3 holds. Indeed, applying the same argument as in the case of polygons in [7], we can show that all functions between the upper and lower bounds, obtained in Lemma 3.3, are i-th barycentric coordinates. ˜ i : P → R is an i-th barycentric coordinate Proposition 3.6. A function λ associated to a vertex if and only if the following bounds hold for all x in P . ˜li (x) ≤ λ ˜ i (x) ≤ L ˜ i (x).

(22)

Proof. Fix i ∈ {0, . . . , n} , and take the vertex v i of the polytope P . We have already shown, in Lemma 3.3, that any barycentric coordinate λi associated to v i must satisfy the bound given in (22). Thus, the bounds (22) are necessary for λi to be an i-th barycentric coordinate with respect to a vertex. To show that the bounds are also sufficient, note first that they immediately imply ˜ i (x) = (1 − αi (x))˜li (x) + αi (x)L ˜ i (x), (∀x ∈ P ), λ (23) for some function αi : P → [0, 1]. By Proposition 3.5 there exist two sets of barycentric coordinates {λ0 , λ1 , . . . , λn } and {β0 , β1 , . . . , βn } such that λi (x) = ˜li (x) ˜ i (x). βi (x) = L Now define the set of functions {γ0 , γ1 , . . . , γn } by γk (x) = (1 − αi (x))λk (x) + αi (x)βk (x), (∀x ∈ P ).

(24)

It is easy to check that {γ0 , γ1 , . . . , γn } form a set of barycentric coordinates ˜ i = γi , this means that λ ˜ i is indeed an i-th barycentric coordinate. and since λ 14

Additional comparisons of barycentric functions can be found in the survey papers of Cueto et al. [4], Sukumar and Tabarraei [39].

4. General properties of barycentric approximation schemes Barycentric approximation schemes and their corresponding approximation errors share a number of desirable properties. The next section enumerates these properties. Let us denote the set of all continuous functions from the polytope P to R by C(P ). The subset of C(P ) formed by all continuous convex functions on P is denoted by K(P ). We will also denote by Rn the subspace of C(P ) spanned by the set of barycentric coordinates λ = {λ0 , λ1 , . . . , λn } . As a consequence of (3) and (4), the linear space Rn contains P1 as a subspace. We now state without proof some essential properties of barycentric approximation schemes. Theorem 4.1. The barycentric approximation scheme B has the following properties. (i) It maps C(P ) into Rn . (ii) It is nonexpansive with respect to the uniform norm k.k∞ , in the sense that if f, g ∈ C(P ) then kB[f ] − B[g]k∞ ≤ kf − gk∞ .

(25)

(iii) It cannot have a second-order approximation. The following statements, which easily follow from a simple calculation, give us some important properties relevant to the incurred approximation error, between f and its barycentric approximation. In the following we will often use the notation: E[f ](x) = E[f, λ, X n ] :=

n X

λi (x)f (xi ) − f (x).

(26)

i=0

Theorem 4.2. The error approximation E has the following properties. (i) It vanishes over affine functions, this means that B is a first-order approximation scheme, i. e., it reproduce exactly all affine functions. 15

(ii) It is nonnegative in K(P ). This means that B approximates convex functions from above. (iii) For any real continuous function f on P, the error E[f ] vanishes at the vertices. Hence, B satisfies the vertex interpolation: B[f ](v i ) = f (v i ), for all vertices of P . Let Ω be a compact convex subset of Rd and let p := {pi }ni=0 be a partition of unity, that is a collection of continuous functions from Ω to the unit interval [0, 1] and whose values sum the unity for all x ∈ Ω. Let {x0 , . . . , xn } be a set of distinct points in Ω. In a partition of unity method, an approximation of a function f defined on Ω is given by pu

A [f ](x) =

n X

pi (x)f (xi ).

(27)

i=0

In the case that p forms a set of barycentric coordinates, Theorem 4.2 shows that the nonnegativity property (ii) is satisfied by the approximation error Apu [f ] − f on K(P ). The next result shows that, in a certain sense, the reverse implication also holds. The main point of what follows it that the nonnegativity property on K(P ) of the approximation error E can be characterized in term of linear precision of the barycentric coordinates. Here, and in what follows, A(P ) will denote the space of all affine real functions on P . Theorem 4.3. Let Ω be a compact convex subset of Rd with positive measure. Assume that there are n + 1 points x0 , . . . , xn ∈ Ω, and a partition of unity {p0 , . . . , pn } on Ω. Then, the following statements are equivalent: (i) {p0 , . . . , pn } is a barycentric coordinates system. (ii) The error E is nonnegative on K(P ). (iii) The error E vanishes on A(Ω). Proof. We have already shown the implication (i) implies (ii) in Theorem 4.2, (ii). For the implication (ii) implies (iii), we use the fact that affine functions and their opposites belong to K(Ω). Then (ii) implies that for any affine function l we have E[l] ≥ 0 and E[−l] ≥ 0. Hence, by virtue of linearity of the operator E we deduce the vanishing property of E on A(P ). Finally, to prove that (iii) implies (i) it suffices to consider the usual projection functions e1 , . . . , ed , ei : x = (x1 , . . . , xd ) → xi . 16

5. Pointwise error estimations In this section, we will discuss the problem of finding upper and lower bounds for the approximation error between any continuous convex function and its barycentric approximation scheme defined by (1). Here, the error induced by the approximation process will be estimated by using a pointwise error measure. For a given convex function f , effective upper and lower bounds for E[f ] at a point x can be given in terms the data values of f , the largest and smallest barycentric coordinates of x. These bounds are presented in equations (34) and (50) of Theorems 5.2 and 5.4, respectively. More precisely, assume we are given a set λ = {λ0 , λ1 , . . . , λn } of barycentric coordinates with respect to a fixed set of (n + 1) pairwise distinct points X n = {xi }ni=0 of Rd , we are interested in estimating, for each x ∈ P, the quantity: E[f ](x) =

n X

λi (x)f (xi ) − f (x), (∀ f ∈ K(P )).

(28)

i=0

Before stating any error estimates, we fix some notation. Throughout, let V (P ) := {v 0 , v 1 , . . . , v nv } denotes the set of vertices of P . If there is no confusion, we will continue to simply denote E[f, λ, X n ] by E[f ]. The centroid w[X n ] of the polytope P = conv(X n ) with respect to X n is defined as the average of the points defining the polytope. Therefore, we have: n 1 X xi , (29) w[X n ] = n + 1 i=0 which we often write as w when X n is clear from context. We note that w is always located inside P . In what follows, we will let e[f ](w) denote n

1 X f (xi ) − f (w). e[f ](w) := n + 1 i=0

(30)

This fundamental quantity will play a significant role in estimating the error E. The reader may verify that under the conditions λi (w) =

1 , i = 0, . . . , n, n+1

(31)

e[f ](w) coincides with E[f ](w). Recall that Lemma 2.1 says that the set of barycentric coordinates λ can be chosen such that the equalities given in 17

equation (31) hold. If f ∈ K(P ), then Theorem 4.2, (ii), asserts that 0 ≤ E[f ]. The lower bound zero is of global nature since it holds for any x ∈ P, f ∈ K(P ) and any barycentric coordinates system. We observe, under the assumption that f is convex, that the following upper bound holds trivially. E[f ](x) ≤ max f (y) − f (x), (∀ f ∈ K(P )). y ∈P

(32)

We have to work a little more in order to improve on the bound given by (32). Indeed, in the subsequent analysis, we shall first prove a remarkable fact that the error E[f ] made in approximating f ∈ K(P ) at any point x ∈ P can be controlled by the error e[f ] at the centroid of X n . Indeed, the next result gives an upper global bound (depending solely on e[f ] at w) and it does not depend on the set of barycentric coordinates used by the approximation operator B: Theorem 5.1. With the above notation, if f is a continuous convex function on P , then we have: 0 ≤ E[f ](x) ≤ (n + 1)e[f ](w), (∀x ∈ P ).

(33)

Proof. To show this we first note that the centroid w can be written as n 1 X xi w = n + 1 i=0 ! n X 1 xi − x + x = n + 1 i=0 ! n n X X 1 = xi − λi (x)xi + x n + 1 i=0 i=0 =

n X 1 − λi (x) i=0

n+1

xi +

1 x. n+1

Which is obviously a convex combination of points of P . Now, Jensen’s inequality immediately implies n X (n + 1)f (w) ≤ (1 − λi (x)) f (xi ) + f (x) i=0

=

n X

f (xi ) − E[f ](x),

i=0

18

which is obviously equivalent to the desired inequality. As the reader may verify, the bound given in Theorem 5.1 is better than that given in (32). The next result gives more precise lower and upper bounds for the error E. To do so, it will be convenient to use two important functions. For any point x ∈ P with barycentric coordinates λ0 (x), . . . , λn (x), we define the largest and smallest barycentric coordinates as follows: λmax (x) := λmin (x) :=

max λi (x)

i=0,...,n

min λi (x).

i=0,...,n

For any given x ∈ P and f ∈ K(P ), Theorem 5.2 below shows how to control the error E[f ](x) in terms of the largest and smallest barycentric coordinates of x and the error e[f ](w). Theorem 5.2. If f is a continuous convex function on P , then, for all x in P , we have: (n + 1)λmin (x)e[f ](w) ≤ E[f ](x) ≤ (n + 1)λmax (x)e[f ](w).

(34)

Proof. Pick any x in P. We distinguish three cases: 1 1 . Then, λi (x) ≤ n+1 , i = 0, . . . , n, therefore • Case I: λmax (x) ≤ n+1
P n 1 1 λi (x) = n+1 , i = 0, . . . , n, since we have i = 0 n+1 − λi (x) = 0. This means that x = w, so inequality (34) obviously holds. 1 • Case II: n+1 ≤ λmin (x). Then, the same argument shows that λi (x) = 1 , i = 0, . . . , n, x = w, and again inequality (34) obviously holds. n+1 1 < λmax (x). • Case III: λmmin (x) < n+1

To simplify notation, let us denote l = (n + 1)λmin (x). Since 0 ≤ l < 1, we can write x as a convex combination of x0 , . . . , xn , namely, x = lw + (1 − l)y, (35) where y=

n X λi (x) − λmin (x) i=0

1−l 19

xi .

(36)

It is not hard to check also that y is a convex combination of the points x0 , . . . , xn . Hence, by using two times Jensen’s inequality we get f (x) ≤ lf (w) +

n X

(λi (x) − λmin (x)) f (xi )

i=0 n l X = lf (w) + f (xi ), λi (x)f (xi ) − n + 1 i=0 i=0 n X

which is equivalent to the left inequality in (34). To show the right-hand side inequality in (34), once again, we can express the centroid w in the form: 1 1 (37) w = x + (1 − )z, r r where, for ease of notation, we set z=

n X λmax (x) − λi (x) i=0

r−1

xi ,

(38)

and r = (n + 1)λmax (x). 1 n+1

Note that, since < λmax (x) then 0 < 1r < 1 and therefore z and w are written as convex combinations of points of P. Hence, by Jensen’s inequality applies to (37), we get after straightforward calculations. f (z) ≥

1 r f (w) − f (x). r−1 r−1

(39)

On the other hand, we have n n X r X f (xi ) = λmax (x)f (xi ) n + 1 i=0 i=0

= ≥ =

n X i=0 n X i=0 n X

λi (x)f (xi ) + (r − 1)

n X λmax (x) − λi (x) i=0

λi (x)f (xi ) + (r − 1)f

λi (x)f (xi ) + (r − 1)f (z). 20

f (xi )

n X λmax (x) − λi (x) i=0

i=0

r−1 r−1

! xi

Now, making use of (39), we rewrite the above inequality as follows: n n X r X f (xi ) ≥ λi (x)f (xi )f (xi ) + rf (w) − f (x). n + 1 i=0 i=0

(40)

Finally, short calculations based on the definition of r establish the right hand side of the inequality in (34). This validates the two bounds given in (34). There are two points that distinguish this result from what has been reported before. First one observes, the simple form of the upper bound that appears at the right side of the inequality (34) is better than those given by (33). Second, Theorem 5.2 yields a significantly better lower bound than the bound 0 provided by Theorem 5.1. It is generally not necessary to use all the data points to estimate the error E. Instead, we can use the data values at the vertices. Let us describe this technique by assuming, for instance, that we only know the values at the vertices of the convex function to be approximated. In this setting, we introduce a particularly useful notation which will be used throughout: nv

1 X f (v i ) − f (wv ), e [f ](w ) := v n + 1 i=0 v

v

(41)

where nv + 1 denotes the number of vertices of the polytope P. As before, we shall also use the notation v

n 1 X v w = v vi, n + 1 i=0

(42)

to denote the centroid of the polytope P with respect to its vertices. From now on the symbol E v will be reserved exclusively to denote the approximation error when the set of points X n is constituted only by vertices of the polytope P . Thus, we have v

v

E [f ](x)

:=

n X

λvi (x)f (v i ) − f (x)

i=0

here {λv0 , λv1 , . . . , λvnv } is a given set of barycentric coordinates of x with respect to the set of vertices {v 0 , . . . , v nv } . With these conventions in mind, 21

the next Lemma tells us how the approximation error E[f ], at any x ∈ P, can be controlled by E v [f ] := E v [f, αv , V (P )], where αv is a set of barycentric coordinates with respect to V (P ). Indeed we have the following result. Lemma 5.3. Fix λ := {λ0 , λ1 , . . . , λn } a set of barycentric coordinates with respect to X n = {x0 , . . . , xn } . Then, there exists a set of barycentric coordinates αv := {α0v , α1v , . . . , αnv v } with respect to {v 0 , . . . , v nv } , such that, for all continuous convex function f, and x in P , we have: E[f, λ, X n ](x) ≤ E v [f, α, V (P )](x).

(43)

Proof. Pick any set of barycentric coordinates {λv0 , λv1 , . . . , λvnv } with respect to V (P ). Clearly, the points x0 , . . . , xn are a convex combination of the vertices of P , that is, v

n X

xi =

λvj (xi )v j , i = 0, . . . , n.

j=0

Hence, by the convexity of f , n X

λi (x)f (xi )

=

i=0

n X

v

λi (x)f

i=0

≤

n X

n X

! λvj (xi )v j

j=0

λi (x)

i=0

nv X

! λvj (xi )f (v j )

(44)

j=0

v

=

n X

αjv (x)f (v j ),

j=0

where αjv (x)

:=

n X

λi (x)λvj (xi ), j = 0, . . . , nv .

i=0

It is easy to check that α0v , α1v , . . . , αnv v

(45)

form a set of barycentric coordinates with respect to {v 0 , . . . , v nv } . Finally, we get the desired result by subtracting f (x) from both sides of inequality (44). 22

Combining the previous Lemma and Theorem 5.2, we have: Theorem 5.4. With the above notation, there exists a set of barycentric coordinates αv := {α0v , α1v , . . . , αnv v } with respect to {v 0 , . . . , v nv } , such that, for all continuous convex function f, and x in P , we have: (nv + 1)λmin (x)ev [f ](wv ) ≤ E[f ](x) ≤ (nv + 1)αmax (x)ev [f ](wv ),

(46)

where αmax (x) := max v αjv (x). j=0,...,n

(47)

Proof. Lemma 5.3 says that E[f ](x) ≤ E v [f ](x).

(48)

Then, we readily see that the right hand side inequality now follows directly from Theorem 5.2, by taking the set of barycentric coordinates given in (45). To establish the left-hand side inequality of (50), it follows from Theorem 5.2 that it is enough to show that (nv + 1) v e [f ](wv ) ≤ e[f ](w). (n + 1)

(49)

To see this, we write the centroid w as follows: n

1 X w := xi n + 1 i=0 v n−n X

=

i=1

yi nv + 1 v + w , n+1 n+1

where y i , i = 1, . . . , n−nv are elements of X n located in int(X n ). Obviously, this representation is a convex combination of elements of P. Therefore, as a consequence of the convexity of f , v

n−n 1 X nv + 1 f (w) ≤ f (y i ) + f (wv ), n + 1 i=1 n+1 ! n nv X X 1 nv + 1 = f (xi ) − f (v i ) + f (wv ). n + 1 i=0 n + 1 i=0

Finally, a simple transformation shows that the last inequality is exactly (49) below. 23

We now use sharp upper and lower bounds, on all barycentric coordinates obtained in Lemma 3.3, to derive new bounds on the approximation error ˜ i be defined as in Lemma 3.3. For each x in P E v . To this end, let ˜li and L define ˜ max (x) := L ˜lmin (x) :=

˜ i (x) max L

i=0,...,nv

min ˜li (x).

i=0,...,nv

As an immediate consequence of Theorem 5.4 and Lemma 3.3, we obtain: Theorem 5.5. If f is a continuous convex function on P , then, for all x in P , we have: ˜ max (x)ev [f ](wv ). (nv + 1)˜lmin (x)ev [f ](wv ) ≤ E v [f ](x) ≤ (nv + 1)L

(50)

6. The case of a convex function f ∈ C 2 (P ) Clearly the lower and upper bound obtained in Theorems 5.1 and 5.2 depend on a choice of the set of barycentric coordinates λ := {λ0 , λ1 , . . . , λn } used to define the barycentric approximation B. This raises the question of what would be a best choice of λ. In this section, we will restrict our study to convex functions f ∈ C 2 (P ). For twice differentiable functions f : P → R, we denote by   ∂f (x) H[f ](x) := ∂xi ∂xj i,j=1,...,d the Hessian matrix of f at x and introduce 2 T D f := sup y H[f ](x)y . sup x∈P y ∈Rd ,ky k=1

(51)

All vectors are understood as column vectors, thus xT , or the transpose of x, is a row vector. The following result is a simple consequence of the classical Taylor’s formula for multivariate real valued function.

24

Proposition 6.1. Suppose f is a C 2 -convex function on P. Then, for all x in P, we have |D2 f | E[k.k2 , λ, X n ](x) 2 n |D2 f | X = λi (x) kx − xi k2 . 2 i=0

0 ≤ E[f, λ, X n ](x) ≤

(52)

Equality is attained for all functions of the form f (x) = a(x) + c, where c ∈ R and a is any affine function. In view of the error estimate (52), it seems desirablePto select the set of barycentric coordinates λ for which E[k.k2 , λ, X n ] := ni=0 λi kx − xi k2 is small. Hence, we want to solve the following constrained optimization problem (Px ): For a fixed but arbitrary x ∈ P minimize λ

E[k.k2 , λ, X n ](x),

subject to the constraints

λi ≥ 0, i = 0, . . . , n,

n X

λi = 1,

i=0

x =

n X

λi xi .

i=0

Since the parametric optimization problem (Px ) depends on x, its solution is also a function of x. We shall call a set of barycentric coordinates λ (global) optimal if for each x ∈ P , λ0 (x), λ1 (x), . . . , λn (x) is a solution of (Px ). Given a triangulation T of the point set X n , then the point x belongs to at least one simplex Sx of the triangulation T. The set of constraints of (Px ) provides (d + 1) constraints, hence, n − d weights can be chosen freely. One possible choice of weights would be to take the barycentric coordinates of the point x with respect to the simplex Sx and set the remaining (n − d) weights to zero. Of course, the resulting choice is not necessarily optimal. We will show that there exists a triangulation for which this choice is optimal. A very natural triangulation DT (P ) of P , that uses only the points of X n as triangulation vertices and such that no point in X n is inside the hypercircumsphere of any simplex of the triangulation. Such a triangulation exists 25

and it is called a Delaunay triangulation of P with respect to X n and it can be obtained as the geometric dual of the Voronoi diagram of X n , see, e. g., [9]. We recall from Section 2 that every Delaunay triangulation generates a set of (linear finite element) barycentric coordinates. The following theorem is an immediate consequence of a known extremal property of Delaunay triangulation. Theorem 6.2. Every set of barycentric coordinates, which is generated by a Delaunay triangulation, is optimal. Proof. Let the function g be defined by g(x) = min E[k.k2 , λ, X n ](x), (∀x ∈ P ). λ By [31, Lemma 10], formulated in our notation and for our situation, at g(x), for a fixed point x, the only nonzero values of λi occur for the vertices of the Delaunay simplex containing the point x. Thus, a necessary condition for minimizing the last expression is that λ is a set of barycentric coordinates, which is generated by a Delaunay triangulation. This completes the proof. Note that if P is a simplex in Rd , then n = d and the problem (Px ) yields for each x a unique solution for the λi . For n > d + 1, which is the case of interest in this paper, the problem (Px ) forms an under-determined system. If points x0 , . . . , xn are in general position, no d + 1 points of X n lie on a common hyperplane in Rd and no d + 2 points lie on a common hypersphere, then Delaunay triangulation is unique. An inspection of the proof of Theorem 6.2 shows that in this case, the global unique solution is described by the set of (linear finite element) barycentric coordinates generated by the unique Delaunay triangulation. When the points are not in general position, the optimization problem above does not have a unique optimal solution. Note that in this case, the problem has an infinity of solutions, indeed, a convex combination of global solutions of (Px ) is also a solution. So the optimal solutions to (Px ) form a convex set. It will be interesting to have a way of comparing optimal solutions and selecting favorable ones. This will be discussed in a forthcoming paper. It should be noted that not every set of optimal barycentric coordinates is generated by a Delaunay triangulation. In fact, consider a two-dimensional square S, Figure 3 shows two different Delaunay triangulations are possible 26

b

b

b

b

b

b

b

b

Figure 3: shows two different Delaunay triangulations of a two-dimensional square.

Figure 4: Delaunay triangulation for a 2-D set of points.

for S. Now every convex combination of these solutions provides a set of optimal barycentric coordinates which are not generated by a triangulation. We should mention that for numerical aspects there exist efficient algorithms for computing a Delaunay triangulation with respect to a finite point set X n whose convex hull is P , see [9]. Figure 4, (a) shows 21 equally spaced nodes on each edge of the boundary of square D := [0, 1]×[0, 1] and 361 scattered points randomly selected in D. Figure (b) shows Delaunay triangulation based on the data points given in (a). Once we have such a triangulation, we can obtain the set of optimal barycentric coordinates. Finally, for the corresponding L1 problem we refer to [16]. 7. Numerical experiments Given a set of scattered data {(xi , yi , fi }N i=1 , which are assumed to be 2 sampled from a convex function f : Ω ⊂ R → R. Taking the N scattered 27

Figure 5: (a) Delaunay triangulated for the domain of f , (b) the graph of f and (c) the graph for the linear interpolation of the data generated from f .

points as nodes, an optimal triangulation mesh, T , is constructed in domain Ω using Delaunay triangulation method. In this section, we present some numerical tests to illustrate the practical performance of the corresponding linear interpolating scheme developed in Section 6. With this aim, we approximate two test convex functions using randomly scattered data points in the indicated domain. Although our final findings and conclusions hold for many more examples we only present here two typical examples. Example 7.1. We take the following convex function from [30]: f (x, y) := x4 + y 4 . The data points are generated from the above function with the restriction of domain D := [0, 1] × [0, 1] . Figure 5, (a) is Delaunay triangulation based on 21 equally spaced nodes on each edge of the boundary of square D and 216 nodes in the square. The nodes in the domain are placed randomly selected from D while the nodes on the boundary is equally spaced. Figure 5, (b) presents the graph of f . Figure 5, (c) describes the graph for the linear interpolation of scattered data generated from the function f. 28

Figure 6: (a) Delaunay triangulated for the domain of g, (b) the graph of g and (c) the graph for the linear interpolation of the data generated from g.

Example 7.2. In the second example the data points is generated from the following test function taken from [23]: g(x, y) = x3 + 5(y 2 − 0.6)2 + 1. As it is mentioned in [23], the data points in this example belong to a surface that models part of a car. The data points are generated from the above function with the restriction of domain D := [0, 1] × [0, 1] , with data values taken from function g. Figure 6, (a) shows the Delaunay triangulation based on 21 equally spaced nodes on each edge of the boundary of square D and 216 nodes randomly selected in the square. Figure 5, (b) is the graph of g while Figure 5, (c) shows the linear interpolation of scattered data generated from the function g. From Figures 5 and 6 it is clear that the convexity of f and g has been preserved and there are no visual differences between the test functions and their linear interpolants.

29

Acknowledgment The author would like to thank the two referees for their useful proposals. This paper is dedicated to Professor Jean-Louis Gout on the occasion of his retirement. References [1] M. Beck, S. Robins, Computing the Continuous Discretely: Integerpoint Enumeration in Polyhedra, Springer, 2007. [2] L. Chen and J. Xu, Optimal Delaunay triangulations, J. Comput. Math. 22 (2004), 299–308. [3] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North Holland, 1979. [4] Cueto, E., Sukumar, N., Calvo, B., Martnez, M.A., Cegoino, J., Doblar, M.: Overview and recent advances in natural neighbor Galerkin methods. Arch. Comput. Methods Eng. 10(4), 307-384 (2003) [5] Desbrun M., Meyer M., Alliez P., Intrinsic parameterizations of surface meshes. Computer Graphics Forum 21 (2002), 209-218. [6] S. Dragomir, A converse of the Jensen inequality for convex mappings of several variables and applications, Acta mathematica vietnamica, vol. 29(1), 77-88, 2004. [7] M. S. Floater, K. Hormann, G. Kos, A general construction of barycentric coordinates over convex polygons, Adv. Comput. Math., 24(1), 311331 (2006). [8] M. S. Floater. Mean value coordinates. Comp. Aided Geom. Design, 20, 19-27, 2003. [9] George P.-L., Borouchaki H., Delaunay Triangulation and Meshing, ´ Application to Finite Elements (Translation from French), Editions Herm`es, Paris, 1998. [10] C. Gout, C k surface approximation from surface patches, Computers Mathematics with Applications, Volume 44, Issues 34, 2002, 389-406 30

[11] J. L. Gout, Construction of a Hermite rational Wachspress type finite element, Computers Mathematics with Applications, Volume 5, 337347, 1979. [12] J.L. Gout, Interpolation error estimates on hermite rational Wachspress type third degree finite element, Computers Mathematics with Applications, Volume 5, Issue 4, 1979, 349-357 [13] Guessab A., Direct and converse results for generalized multivariate Jensen-type inequalities, Journal of Nonlinear and Convex Analysis, V. 13, 4, 777-797, 2013 [14] Guessab A., Nouisser O., Peˇcari´c J., A multivariate extension of an inequality of Brenner-Alzer, Archiv der Mathematik, Volume 98, Number 3, 277–287, 2012. [15] Guessab A., Nouisser O., Schmeisser G., Multivariate approximation by a combination of modified Taylor polynomials, Journal of Computational and Applied Mathematics 196 (2006) 162-179 [16] Guessab A., Schmeisser G., Negative definite cubature formulae, extremality and Delaunay triangulation, Constr. Approx., Volume 31, 95113, 2010. [17] Guessab A., Schmeisser G., Construction of positive definite cubature formulae and approximation of functions via Voronoi tessellations, Adv Comput Math, 32, 25-41, 2010. [18] Guessab A., Nouisser O., Schmeisser G., A definiteness theory for cubature formulae of order two, Constr. Approx., 24, 263-288, 2006. [19] Guessab A., Schmeisser G., Convexity results and sharp error estimates in approximate multivariate integration, Math. Comp. 73 (2004), 1365– 1384. [20] Guessab A., Nouisser O., Schmeisser G., Enhancement of the algebraic precision of a linear operator and consequences under positivity, Positivity, 2009, Volume 13, Issue 4, 693-707. [21] Guessab A., Schmeisser G., Sharp Error Estimates for Interpolatory Approximation on Convex Polytopes, SIAM J. Numer. Anal., 43(3), 909-923, 2005 31

[22] R.L. Hardy. Multiquadric equations of topography and other irregular surfaces. J. Geophys. Res, 76(8), 1905-1915, 1971. [23] A. Li, Convexity preserving interpolation, Computer Aided Geometric Design 16 (1999), 127-147. [24] Ju T., Schaefer S., Warren J., Mean value coordinates for closed triangular meshes. ACM Trans. Graph. 24, 3 (2005), 561-566. [25] Ju T., Schaefer S., Warren J., Desbrun M., A geometric construction of coordinates for convex polyhedra using polar duals. In Proceedings of the Symposium on Geometry Processing (2005), 181-186. [26] Ju T., Warren J., General Constructions of Barycentric Coordinates in a Convex Triangular Polyhedron. Tech. rep., Washington University in St. Louis, Nov. 2005. [27] Kalman J. A., Continuity and convexity of projections and barycentric coordinates in convex polyhedra, Pacific J. of Math., II, 1017-1022 (1961) [28] M. Klaricic Bakula, J. Pecaric, J. Peric, On the converse Jensen inequality, Applied Mathematics and Computation 218 (2012) 6566-6575 [29] J. M. Pe˜ na, B-splines and optimal stability, Mathematics of Computation, 66 (1997), 1555- 1560. [30] R.J. Renka, Algorithm 833: CSRFPACK-Interpolation of Scattered Data with a C 1 Convexity-Preserving Surface, ACM Transactions on Mathematical Software 30(2) (2004), 200-211. [31] V. T. Rajan, Optimality of the Delaunay Triangulation in Rd , Discrete Comput Geom 12, 189-202 (1994) [32] Schreiner J., Asirvatham A., Praun E., Hoppe H., Inter-surface mapping. ACM Trans. Graph. 23, 3 (2004), 870-877. [33] Schrijver A., Theory of Linear and Integer Programming. Wiley, (1987). [34] Sederberg T. W., Parry S. R., Free-form deformation of solid geometric models. In Computer Graphics, Proceedings of ACM SIGGRAPH (1986), 151-160. 32

[35] D. Shepard. A two-dimensional interpolation function for irregularlyspaced data. In Proceedings of the 1968 23rd ACM national conference, pages 517-524. ACM New York, NY, USA, 1968. [36] S. Simic, Best possible bound for Jensens inequality, Applied mathematics and computation, vol. 215, 2224-2228, 2009. [37] G. Strang and G. Fix. An Analysis of the Finite Element Method. Prentice-Hall, Englewood Cliffs, N.J., 1973. [38] Sukumar N., Malsch E. A., Recent advances in the construction of polygonal interpolants. Arch. Comput. Meth. Engng. 13, 1 (2006), 129-163. [39] Sukumar, N., Tabarraei, A., Conforming polygonal finite elements. Int. J. Numer. Methods Eng. 61(12), 2045-2066 (2004) [40] Wachspress E. L., A Rational Finite Element Basis, vol. 114. Academic Press, New York, 1975. [41] Ziegler G. M., Lectures on Polytopes, Springer Verlag, New York (1995) 67, Elsevier Science, (2005)

33