Interpolation with quintic Powell–Sabin splines

Applied Numerical Mathematics 62 (2012) 620–635

Contents lists available at SciVerse ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Interpolation with quintic Powell–Sabin splines Hendrik Speleers Department of Computer Science, Katholieke Universiteit Leuven, Celestijnenlaan 200A, B-3001 Leuven, Belgium

a r t i c l e

i n f o

Article history: Received 31 January 2011 Received in revised form 1 December 2011 Accepted 16 January 2012 Available online 3 February 2012 Keywords: Hermite interpolation Quintic Powell–Sabin splines C 2 macro-elements Normalized B-spline representation

a b s t r a c t We discuss local Hermite interpolation by C 2 quintic Powell–Sabin splines represented in a normalized B-spline basis. We derive explicit formulae for the spline coefficients in this B-spline representation to interpolate given Hermite data. As part of the analysis, we show how tensor algebra can be used to describe polynomials in Bernstein–Bézier form and to simplify their manipulation. © 2012 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction The construction of smooth functions interpolating scattered (Hermite) data is a problem that arises in application areas like visualization or reconstruction of smooth objects for numerical simulation. This problem can be addressed by working with smooth polynomial spline functions defined on a triangulation whose vertices are located at the data points. For the construction of smooth splines with a low polynomial degree one often considers triangulations with a particular refinement. For example, the Powell–Sabin refinement splits every triangle into six subtriangles [13]. Interpolation and approximation schemes based on C 1 quadratic Powell–Sabin splines can be found in [4,11,16,17] and shape-preserving schemes in [12,18,21,22]. In this paper we address the problem of Hermite interpolation using suitable subspaces of C 2 quintic splines defined over triangulations with a Powell–Sabin refinement. Such spline spaces have been constructed and analyzed by several authors, see [1,7–10,15]. The standard Hermite spline basis is the natural choice of basis to face the addressed interpolation problem. An alternative basis for particular quintic Powell–Sabin splines has been developed recently in [20], possessing some interesting properties. The basis functions (so-called Powell–Sabin B-splines) have a local support, they form a partition of unity and are nonnegative. This B-spline representation allows a natural definition of control polynomials and control points, which are useful for computer aided geometric design. A spline in such a representation can also be evaluated in a stable way, using only a sequence of simple convex combinations. These B-splines can be seen as specific C 2 triangular finite (or macro-) elements. A similar basis has been constructed for quadratic Powell–Sabin splines in [3]. The purpose of this paper is to solve the Hermite interpolation problem with respect to the quintic Powell–Sabin Bspline basis. Basically, we need to find the conversion between the Hermite basis and the B-spline basis. To simplify the computation, we make use of (multilinear) tensor algebra. We provide a compact representation of bivariate polynomials by means of tensors. Tensors are a natural generalization of matrices. This tensor representation allows us to obtain explicit compact formulae for the coefficients of the spline interpolant represented in terms of Powell–Sabin B-splines. This interpolation problem requires data at the vertices of the triangulation and at the Powell–Sabin split points. We will also discuss

E-mail address: [email protected]. 0168-9274/$36.00 © 2012 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2012.01.008

H. Speleers / Applied Numerical Mathematics 62 (2012) 620–635

621

Fig. 1. Schematic representation of the Bézier ordinates of a cubic bivariate polynomial.

interpolation using some reduced spline spaces. The corresponding reduced interpolation problems only need data at the vertices of the triangulation. The paper is organized as follows. In Section 2 we recall some relevant concepts of Bernstein–Bézier polynomials on triangles, and we point out the relation to tensors in the next section. Section 4 summarizes the construction of the quintic Powell–Sabin B-spline basis. In Section 5 we solve the Hermite interpolation problem with respect to the quintic Powell– Sabin B-splines, and in Section 6 we discuss interpolation in some reduced spline spaces. Section 7 presents a numerical example. We end in Section 8 with some concluding remarks. 2. Bivariate polynomials on triangles 2.1. Bernstein–Bézier representation of polynomials Let T ( V 1 , V 2 , V 3 ) be a non-degenerate triangle. Any point P in the plane of the triangle can be uniquely expressed in terms of the barycentric coordinates τ = (τ1 , τ2 , τ3 ), such that

P=

3 

τi V i , and τ1 + τ2 + τ3 = 1.

(2.1)

i =1

Given two points P 1 and P 2 in the plane of the triangle, the barycentric direction δ = (δ1 , δ2 , δ3 ) of the vector P 2 − P 1 with respect to T is defined as the difference of the barycentric coordinates of both points. If the Euclidean distance  P 2 − P 1 2 = 1, then δ is called a unit barycentric direction. A barycentric direction δ always satisfies δ1 + δ2 + δ3 = 0. Let Πd denote the linear space of bivariate polynomials of total degree less than or equal to d. Any polynomial p d ∈ Πd on triangle T has a unique Bernstein–Bézier representation



pd (τ ) =

b i jk B dijk (τ ),

(2.2)

i + j +k=d

with

B dijk (τ ) =

d! i ! j !k!

(τ1 )i (τ2 ) j (τ3 )k

(2.3)

the Bernstein basis polynomials of degree d. The coefficients b i jk are called Bézier ordinates. The Bézier domain points ξi jk j

are defined as the points with barycentric coordinates ( di , d , dk ). By associating the Bézier ordinates b i jk with the corresponding Bézier domain points ξi jk , one can display the Bernstein–Bézier representation schematically as in Fig. 1 for the case d = 3. We refer to [5] for more details. 2.2. The blossom of polynomials Let Pd (τ 1 , τ 2 , . . . , τ d ) be the blossom (or polar form) of polynomial pd (τ ). The blossom is completely characterized by the following three properties (see, e.g., [14,19] for more details): 1. The blossom is symmetric with respect to any permutation π of its d arguments. 2. The blossom is multi-affine, i.e. affine in each of its d arguments. 3. The blossom evaluated on its diagonal equals the polynomial p d , i.e., Pd (τ , τ , . . . , τ ) = pd (τ ).

622

H. Speleers / Applied Numerical Mathematics 62 (2012) 620–635

We consider the case d = 3 in more detail, as we are particularly interested in it further on. Let 1, 2, 3, then the blossom of p 3 on domain triangle T ( V 1 , V 2 , V 3 ) can be expressed as

τ l = (τl,1 , τl,2 , τl,3 ), l =

P3 (τ 1 , τ 2 , τ 3 )

= b300 (τ1,1 τ2,1 τ3,1 ) + b030 (τ1,2 τ2,2 τ3,2 ) + b003 (τ1,3 τ2,3 τ3,3 ) + b210 (τ1,2 τ2,1 τ3,1 + τ1,1 τ2,2 τ3,1 + τ1,1 τ2,1 τ3,2 ) + b201 (τ1,3 τ2,1 τ3,1 + τ1,1 τ2,3 τ3,1 + τ1,1 τ2,1 τ3,3 ) + b021 (τ1,3 τ2,2 τ3,2 + τ1,2 τ2,3 τ3,2 + τ1,2 τ2,2 τ3,3 ) + b120 (τ1,1 τ2,2 τ3,2 + τ1,2 τ2,1 τ3,2 + τ1,2 τ2,2 τ3,1 ) + b102 (τ1,1 τ2,3 τ3,3 + τ1,3 τ2,1 τ3,3 + τ1,3 τ2,3 τ3,1 ) + b012 (τ1,2 τ2,3 τ3,3 + τ1,3 τ2,2 τ3,3 + τ1,3 τ2,3 τ3,2 ) + b111 (τ1,1 τ2,2 τ3,3 + τ1,1 τ2,3 τ3,2 + τ1,2 τ2,3 τ3,1 + τ1,2 τ2,1 τ3,3 + τ1,3 τ2,1 τ3,2 + τ1,3 τ2,2 τ3,1 ),

(2.4)

with b i jk , i + j + k = 3, the Bézier ordinates of p 3 on T . With the aid of the blossom one can easily find the Bernstein–Bézier representation of polynomial p d with respect to another (finer) triangle. Consider a second triangle, whose vertices W l , l = 1, 2, 3, have σ l = (σl,1 , σl,2 , σl,3 ) as barycentric coordinates with respect to the former triangle T ( V 1 , V 2 , V 3 ). The Bézier ordinates di jk of pd related to the new triangle can then be obtained as

di jk = Pd (σ 1 , . . . , σ 1 , σ 2 , . . . , σ 2 , σ 3 , . . . , σ 3 ),





 



i times

 

j times



(2.5)



k times

for all i + j + k = d. Thus, subdivision of a polynomial onto a partition of the original domain triangle can be compactly described through the blossom. Blossoming can also be used to describe directional derivatives of a polynomial (2.2). The r-th order directional derivative of polynomial pd with respect to the unit barycentric directions δ i , i = 1, . . . , r, on triangle T ( V 1 , V 2 , V 3 ) can be compactly expressed as

D rδ 1 ,...,δr pd (τ ) =

d! Pd (τ , . . . , τ , δ 1 , . . . , δ r ). (d − r )!

(2.6)

Here, the r unit barycentric directions δ i and (d − r ) times the barycentric coordinates blossom.

τ are taken as the arguments of the

3. Polynomials and tensors 3.1. Tensor algebra A tensor, also known as a multiway array or a multidimensional matrix, is a higher-order generalization of a vector (first order tensor) and a matrix (second order tensor). The order of tensor A ∈ R I 1 × I 2 ×···× I N is N. An element of A is denoted by ai 1 i 2 ...i N where 1  in  I n . The n-mode vectors of A are the I n -dimensional vectors obtained from A by varying the index in while keeping the other indices fixed. In this terminology, column vectors of a matrix are referred to as 1-mode vectors and row vectors as 2-mode vectors. A tensor can be represented in matrix form by unfolding the tensor along a single dimension. The matrix unfolding A (n) ∈ R I n ×( I n+1 ... I N I 1 ... I n−1 ) of a tensor A ∈ R I 1 ×···× I N is a matrix composed of all the n-mode vectors of the tensor. It contains the element ai 1 i 2 ...i N at the position with row number in and column number equal to

1+

 n −1 

( i j − 1)

j =1

n −1 k = j +1



Ik

+

N 



( i j − 1)

j =n+1

N  k = j +1

Ik

n −1



Ik .

k =1

A generalization of the product of two matrices is the product of a tensor and a matrix. The n-mode product of a tensor A ∈ R I 1 ×···× I N and a matrix U ∈ R J n × I n is defined as the tensor B ∈ R I 1 ×···× I n−1 × J n × I n+1 ×···× I N of which the elements are given by

b i 1 ...in−1 jn in+1 ...i N =

In 

ai 1 ...in−1 in in+1 ...i N u jn in .

(3.1)

i n =1

This product is denoted by

B = A ×n U .

(3.2)

Using matrix unfolding, (3.2) can be written as a classical matrix–matrix product, i.e.,

B (n) = U A (n) . More details on tensor algebra can be found in [2].

(3.3)

H. Speleers / Applied Numerical Mathematics 62 (2012) 620–635

623

3.2. A tensor representation for blossoms of polynomials We now explain how tensor algebra can be used to describe and manipulate polynomials in Bernstein–Bézier form. We just focus on the case d = 3, but the approach can be generalized to any degree d. We define the symmetric (3 × 3 × 3)-tensor B by its matrix unfolding B (1) as

⎡

b300 B (1) = ⎣ b210 b201

b210 b120 b111

b201 b210 b111 b120 b102 b111

b120 b030 b021

b111 b201 b021 b111 b012 b102

b111 b021 b012

⎤

b102 b012 ⎦ . b003

(3.4)

This tensor consists of the Bézier ordinates b i jk , i + j + k = 3, of the polynomial p 3 (τ ) in (2.2) defined on domain triangle T ( V 1 , V 2 , V 3 ). Note that the three matrix unfoldings of B are equal, i.e., B (1) = B (2) = B (3) . The barycentric coordinates τ l = (τl,1 , τl,2 , τl,3 ), l = 1, 2, 3, with respect to triangle T can be arranged into the row vectors T l = [τl,1 τl,2 τl,3 ], l = 1, 2, 3. It is then easy to check that the blossom (2.4) can be written as the following product:

P3 (τ 1 , τ 2 , τ 3 ) = B ×1 T 1 ×2 T 2 ×3 T 3 .

(3.5)

From the diagonality property of the blossom, it follows that a polynomial p d in Bernstein–Bézier form (2.2) can be written as

pd (τ 1 ) = B ×1 T 1 ×2 T 1 ×3 T 1 .

(3.6)

We now use this tensor representation to solve the following Hermite interpolation problem: find the Bézier ordinates b i jk , i + j + k = 3, of the polynomial p 3 (x, y ) that satisfies

∂ a+b p 3 ( P ) = h xa yb , ∂ xa ∂ yb

a  0, b  0, a + b  3 ,

(3.7)

for a given point P and a given set of h xa yb -values. We could determine the ten unknowns b i jk by formulating and solving a linear system of ten equations. Alternatively, we use the following tensor approach. Let

hˆ xa yb =

(3 − a − b)! 6

(3.8)

h xa yb ,

we define the (3 × 3 × 3)-tensor H by its matrix unfolding H (1) as

⎡ ⎢

hˆ

H (1) = ⎣ hˆ x

hˆ x

hˆ x2

hˆ y

hˆ x

hˆ x2

hˆ xy hˆ x2

hˆ xy

hˆ x3

hˆ y

hˆ x2 y hˆ xy

hˆ y hˆ xy hˆ y 2 hˆ xy hˆ x2 y hˆ xy 2 hˆ y 2

hˆ xy

hˆ y 2

hˆ x2 y hˆ

⎤ ⎥

hˆ xy 2 ⎦ . hˆ

xy 2

(3.9)

y3

Assume the point P has barycentric coordinates τ = (τ1 , τ2 , τ3 ) with respect to the domain triangle. We combine these barycentric coordinates together with the unit barycentric directions δ x = (δx,1 , δx,2 , δx,3 ) and δ y = (δ y ,1 , δ y ,2 , δ y ,3 ) along the x- and y-direction into the matrix M as

⎡

τ1

M = ⎣ δx,1

τ2

τ3

⎤

δx,2 δx,3 ⎦ . δ y ,1 δ y ,2 δ y ,3

(3.10)

Using (2.6) and (3.5), we then obtain that interpolation problem (3.7) can be reformulated as

H = B ×1 M ×2 M ×3 M .

(3.11)

Because

B = B ×n M ×n M −1 ,

n = 1, 2, 3,

the solution of our problem is found as

B = H ×1 M −1 ×2 M −1 ×3 M −1 .

(3.12)

Instead of solving a linear system of ten equations, by using the tensor approach we only need to compute the inverse of a (3 × 3)-matrix M and the tensor product (3.12). It turns out that the inverse of M can be easily obtained. Let (x P , y P ) be the Cartesian coordinates of the point P . From the definition of barycentric coordinates, see (2.1), we derive the relation

⎡

τ1

τ2

⎣ δx,1 δx,2 δ y ,1 δ y ,2

τ3

⎤⎡

x1 δx,3 ⎦ ⎣ x2 δ y ,3 x3

y1 y2 y3

⎤

⎡

1 xP 1⎦=⎣ 1 1 0

yP 0 1

⎤

1 0 ⎦. 0

624

H. Speleers / Applied Numerical Mathematics 62 (2012) 620–635

We then obtain that the inverse of M can be computed as

⎡

x1 M −1 = ⎣ x2 x3

⎤⎡

y1 y2 y3

1 xP 1 ⎦⎣ 1 1 0

yP 0 1

⎤−1

1 0⎦ 0

⎡

1 =⎣1 1

x1 − x P x2 − x P x3 − x P

⎤

y1 − y P y2 − y P ⎦ . y3 − y P

(3.13)

Combining (3.12) and (3.13), it follows that the solution of problem (3.7) can be compactly written as a tensor–matrix product. Calculating this product gives us a more condensed analytical formula than in case one would solve problem (3.7) in the classical way by using a linear system of ten equations. We provide an explicit expression of such a solution in Section 5.1 in the context of quintic Powell–Sabin spline interpolation. 4. Quintic Powell–Sabin splines 4.1. The PS5-spline space Given a simply connected subset Ω ⊂ R2 with polygonal boundary ∂Ω , let be a conforming triangulation of Ω . Let n v , nt and ne be the number of vertices, triangles and edges in , respectively. A Powell–Sabin (PS-) refinement ∗ of partitions each triangle T j ∈ into six smaller triangles in the following way: 1. Choose an interior point Z j in each triangle T j , so that if two triangles Ti and T j have a common edge, then the line joining Z i and Z j intersects the common edge at a point between its vertices. 2. Join each point Z j to the vertices of T j . 3. For each edge of the triangle T j (a) which is common to a triangle Ti : join Z j to the intersection point of that edge and the line Z i − Z j . (b) which belongs to the boundary ∂Ω : join Z j to an arbitrary point on that edge. The obtained subtriangles are denoted by T ∗ ∈ ∗ . The space of piecewise quintic polynomials on ∗ with global C 2 continuity is given by









S 52 ∗ = s ∈ C 2 (Ω): s|T ∗ ∈ Π5 , T ∗ ∈ ∗ .

(4.1)

We consider a particular subspace of S 52 ( ∗ ) with additional smoothness around some vertices and edges. Let nt ∗ ∗ be the set of vertices in , let Z ∗ = { Z i }i = 1 be the set of triangle split points in , and let E be the set of ∗

that connect a triangle split point Z i to an edge split point. The quintic Powell–Sabin (PS5-) spline space is

















Sˆ 52 ∗ = s ∈ S 52 ∗ : s ∈ C 3 ( W ), W ∈ V ∪ Z ∗ ; s ∈ C 3 (e ), e ∈ E ∗ .

V = { V i }ni =v 1

all edges in defined as

(4.2)

Here, C μ ( W ) means that the polynomials on triangles in ∗ sharing the vertex W have common derivatives up to order μ at that vertex. Analogously, C μ (e ) means that the polynomials on triangles in ∗ sharing the edge e have common derivatives up to order μ along that edge. In [9,10] it is shown that the dimension of the space Sˆ 52 ( ∗ ) is equal to 10n v + nt . The following Hermite interpolation

problem can then be considered [10]: there exists a unique spline s(x, y ) ∈ Sˆ 52 ( ∗ ) such that

∂ a+b s( V l ) = f xa yb ,l , ∂ xa ∂ yb

l = 1 , . . . , n v , a  0, b  0, a + b  3 ,

(4.3a)

and

s ( Z m ) = gm ,

m = 1, . . . , nt ,

(4.3b)

for any given set of f xa yb ,l -values and gm -values. 4.2. A normalized B-spline representation A procedure for the construction of a normalized basis for Sˆ 52 ( ∗ ) was developed in [20]. Ten basis functions { B iv, j (x, y ), j = 1, . . . , 10} are associated with each vertex V i of and one basis function B kt (x, y ) is associated with each triangle Tk of , so that every PS5-spline s(x, y ) ∈ Sˆ 2 ( ∗ ) can be represented as 5

s(x, y ) =

nv  10  i =1 j =1

c iv, j B iv, j (x, y ) +

nt 

ckt B kt (x, y ).

k =1

In addition, the basis functions B iv, j (x, y ) and B kt (x, y ) have a local support, and

(4.4)

H. Speleers / Applied Numerical Mathematics 62 (2012) 620–635

B iv, j (x, y )  0, B kt (x, y )  0,

nv  10 

B iv, j (x, y ) +

i =1 j =1

nt 

625

B kt (x, y ) = 1,

(4.5)

k =1

for all (x, y ) ∈ Ω . A system satisfying these properties is often called a blending system. The functions B iv, j (x, y ) and B kt (x, y ) will be referred to as Powell–Sabin B-splines with respect to vertex V i and triangle Tk , respectively. We now summarize from [20] the construction of such a normalized basis. We first define for each vertex V i ten linearly independent decuplets {αiab , j , 0  a + b  3}, j = 1, . . . , 10. Let M i be the molecule (also called 1-ring) of vertex V i , defined as the union of all triangles in that contain V i . The decuplets {αiab , j , 0  a + b  3} are constructed as follows: 1. For each vertex V i in , identify the corresponding PS5-points. They are defined as

S il =

2 5

Vi +

3 5

Vl,

(4.6)

for all vertices V l that are situated at the boundary of the molecule M i of V i . The vertex V i itself is also a PS5-point. 2. For each vertex V i , find a triangle t i ( Q iv,1 , Q iv,2 , Q iv,3 ) that contains all the PS5-points of V i . Denote its vertices by Q iv, j = ( X iv, j , Y iv, j ). The triangles t i , i = 1, . . . , n v , are called PS5-triangles. We remark that PS5-triangles are not uniquely defined. In [20] an optimization strategy was proposed to select triangles of minimal area. 3. For each vertex V i , consider the ten Bernstein basis polynomials (2.3) of degree three defined on PS5-triangle t i . Each decuplet {αiab , j , 0  a + b  3}, j = 1, . . . , 10, is related to the function and derivative values up to order three of one of these cubic Bernstein basis polynomials evaluated at the point V i :

αiab,1 = αiab,2 = αiab,3 = αiab,4 = ...,

 a+b

20

3

(5 − a − b) (4 − a − b) 5  a+b 20

3

(5 − a − b) (4 − a − b) 5  a+b 20

3

(5 − a − b) (4 − a − b) 5  a+b 20

3

(5 − a − b) (4 − a − b) 5

αiab,10 =

∂ a+b 3 B ( V i ), ∂ xa ∂ yb 300 ∂ a+b 3 B ( V i ), ∂ xa ∂ yb 030 ∂ a+b 3 B ( V i ), ∂ xa ∂ yb 003 ∂ a+b 3 B ( V i ), ∂ xa ∂ yb 210

 a+b

20

3

(5 − a − b) (4 − a − b) 5

∂ a+b 3 B ( V i ), ∂ xa ∂ yb 111

(4.7)

for all 0  a + b  3. t v Given the sets {αiab , j , 0  a + b  3}, the basis functions B i , j (x, y ) and B k (x, y ) are then constructed in the following way:

1. The B-spline B iv, j (x, y ) with respect to vertex V i is defined as the unique solution of interpolation problem (4.3) with all f xa yb ,l = 0, except for l = i, where f xa yb ,i = αiab , j , and with all gm = 0, except for any m such that Tm ∈ M i , where m gm = βim, j . The values αiab , j are determined by the PS5-triangle t i , see (4.7). The values βi , j are specified with the aid of the Bernstein–Bézier representation of the B-spline as follows. We consider the macro-triangle Tk ( V 1 , V 2 , V 3 ) in the molecule M 1 of V 1 with its PS-refinement as shown in Fig. 2, and the B-spline B 1v, j (x, y ) associated with vertex V 1 . On each of the six subtriangles, the spline is a quintic polynomial that can be represented in its Bernstein–Bézier formulation, i.e. with d = 5 in Eqs. (2.2) and (2.3). The corresponding Bézier ordinates are schematically represented in Fig. 3(a). Bézier ordinates that are known to be zero are indicated by open bullets !. Because of the C 3 -smoothness at vertex V 1 , the Bézier ordinates in the 3-disk around vertex V 1 are completely determined by the values {α1ab, j , 0  a + b  3}. The Bézier ordinates in the 3-disk around the split point Z k are found by subdivision of a single cubic polynomial p 3v defined on the triangle with vertices

P1 =

3 5

V1 +

2 5

Zk ,

P2 =

3 5

V2 +

2 5

Zk ,

P3 =

3 5

V3 +

2 5

Zk .

(4.8)

The Bernstein–Bézier representation of this polynomial p 3v has Bézier ordinates that are all zero, except for the three ordinates in the neighbourhood of P 1 (see Fig. 3(a) right). Their values are also determined by 3

α1ab, j . The values of the

remaining non-zero Bézier ordinates in Fig. 3(a) can be computed from the C -smoothness conditions across the interior edges Z k − R 12 and Z k − R 31 . The value of β1k, j is obtained by evaluating the cubic polynomial p 3v at split point Z k .

626

H. Speleers / Applied Numerical Mathematics 62 (2012) 620–635

Fig. 2. A Powell–Sabin refinement of a triangle Tk ( V 1 , V 2 , V 3 ) drawn with dashed lines.

Fig. 3. Schematic representation of the Bézier ordinates of (a) a B-spline with respect to a vertex and (b) a B-spline with respect to a triangle. Bézier ordinates that are known to be zero are indicated by open bullets !; other Bézier ordinates are indicated by filled bullets ". The Bernstein–Bézier representations of the cubic polynomials p 3v and pt3 used for the computation of the ordinates in the 3-disk around the triangle split point are illustrated on the right.

2. The B-spline B kt (x, y ) with respect to triangle Tk is defined as the unique solution of interpolation problem (4.3) with all f xa yb ,l = 0 and with all gm = 0, except for m = k, where gk = βk = 0. The value βk is determined with the aid of the Bernstein–Bézier representation of the B-spline. We consider again the triangle Tk ( V 1 , V 2 , V 3 ) shown in Fig. 2. The corresponding Bézier ordinates are schematically represented in Fig. 3(b). Most of these ordinates are zero. The Bézier ordinates in the 3-disk around the split point Z k are found by subdivision of a cubic polynomial pt3 defined on the triangle with vertices (4.8). The Bézier ordinates in the Bernstein–Bézier representation of this polynomial pt3 are all zero, except for the central ordinate b111 = 1. The value of βk is obtained by evaluating pt3 at split point Z k . For more details we refer to [20]. It is easy to see that a B-spline B iv, j is zero outside the molecule M i of vertex V i , and that a B-spline B kt vanishes outside triangle Tk . The specific choice of the values αiab , j in (4.7) and the fact that the PS5-triangles contain the PS5-points are the main ingredients to obtain the properties (4.5), see [20].

H. Speleers / Applied Numerical Mathematics 62 (2012) 620–635

627

Fig. 4. Schematic representation of the B-spline coefficients c iv, j , j = 1, . . . , 10, with respect to PS5-triangle t i ( Q iv,1 , Q iv,2 , Q iv,3 ).

Since each B-spline B iv, j (x, y ), j = 1, . . . , 10, with respect to vertex V i is related to a Bernstein basis polynomial of degree three, see (4.7), the coefficients c iv, j of a PS5-spline given by (4.4) can be represented schematically as in Fig. 4 with respect to PS5-triangle t i , in correspondence to the representation in Fig. 1. The cubic polynomial defined on t i with the coefficients c iv, j as Bézier ordinates is called the control polynomial with respect to vertex V i , and it is denoted by T i (x, y ). From the local support of the B-splines and from (4.7) it follows that a PS5-spline s(x, y ) is related to its control polynomial T i (x, y ) by

 a+b a+b ∂ a+b 20 3 ∂ s( V i ) = T i ( V i ), a b (5 − a − b) (4 − a − b) 5 ∂x ∂ y ∂ xa ∂ yb

(4.9)

implying that this control polynomial is tangent to the PS5-spline at vertex V i . In the next section we will derive explicit expressions for the PS5-spline coefficients in B-spline representation (4.4) to satisfy the Hermite interpolation problem (4.3). 5. PS5-spline interpolation 5.1. Hermite interpolation at the vertices We first show how to satisfy the interpolation conditions (4.3a) at the vertices. Since the spline coefficients c iv, j can be interpreted as Bézier ordinates defined over the PS5-triangles, we can use the tensor approach described in Section 3.2 to solve this interpolation problem. We arrange the spline coefficients c iv, j , j = 1, . . . , 10, associated with vertex V i = (xi , y i ) into a (3 × 3 × 3)-tensor C. Its matrix unfolding C (1) is given by

⎡

⎤

c iv,1

c iv,4

c iv,5

c iv,4

c iv,7

c iv,10

c iv,5

c iv,10

c iv,8

C (1) = ⎣ c iv,4

c iv,7

c iv,10

c iv,7

c iv,2

c iv,6

c iv,10

c iv,6

c iv,9 ⎦ .

c iv,5

c iv,10

c iv,8

c iv,10

c iv,6

c iv,9

c iv,8

c iv,9

c iv,3

⎢

⎥

(5.1)

The barycentric coordinates of V i with respect to its PS5-triangle t i are denoted by (γi ,1 , γi ,2 , γi ,3 ). The unit barycentric y y y directions along the x- and y-direction with respect to t i are denoted by (γix,1 , γix,2 , γix,3 ) and (γi ,1 , γi ,2 , γi ,3 ), respectively. They can be arranged into the matrix G as

⎤

⎡

γi,1 γi,2 γi,3 ⎢γx γx γx ⎥ G = ⎣ i ,1 i ,2 i ,3 ⎦ . γi,y1 γi,y2 γi,y3

(5.2)

By the definition of PS5-triangle t i and by (3.13), the inverse of G can be written as

⎡ ⎢

1

G −1 = ⎣ 1 1 Let

⎤

X iv,1 − xi

Y iv,1 − y i

X iv,2 − xi

Y iv,2 − y i ⎦ . Y iv,3 − y i

X iv,3

− xi

⎥

(5.3)

628

H. Speleers / Applied Numerical Mathematics 62 (2012) 620–635

ˆf

xa y b

=

 a+b (5 − a − b)! 5 120

3

f xa yb ,i ,

(5.4)

we define the (3 × 3 × 3)-tensor F by its matrix unfolding F (1) as

⎡

ˆf

⎢

ˆ F (1 ) = ⎢ ⎣ fx

ˆf y

ˆf x

ˆf y

ˆf x

ˆf

x2

ˆf xy

ˆf y

ˆf xy

ˆf

ˆf

ˆf xy

ˆf

ˆf

x3

ˆf

x2 y

ˆf xy

ˆf

x2 y

ˆf

ˆf

ˆf xy

x2 y

ˆf

xy 2

ˆf

ˆf

xy 2

ˆf

x2

ˆf xy

y2

x2

ˆf

y2

y2

xy 2

⎤ ⎥ ⎥. ⎦

(5.5)

y3

Note that we need to use a different factor in (5.4) than in (3.8), because we also have to take into account the constant in the relation (4.7) between the derivatives of the Powell–Sabin B-splines and the Bernstein basis polynomials defined over the PS5-triangle t i . Combining (3.7)–(3.11) with (4.9), it follows that interpolation problem (4.3a) at vertex V i can be reformulated as

F = C ×1 G ×2 G ×3 G .

(5.6)

Thus,

C = F ×1 G −1 ×2 G −1 ×3 G −1 .

(5.7)

c iv, j

The values of the coefficients can then be computed using (3.1), (5.1), (5.3), (5.5) and (5.7). For a compact notation of these expressions, we make use of a function R i ( P 1 , P 2 , P 3 ) with three points P l = (x P l , y P l ), l = 1, 2, 3, as arguments. This function is defined as

R i ( P 1, P 2, P 3) = f i +

+ + +

1

3 1

 (x P 1 − xi ) + (x P 2 − xi ) + (x P 3 − xi ) f x,i

 ( y P 1 − y i ) + ( y P 2 − y i ) + ( y P 3 − y i ) f y ,i

3 5  36 5 

 (x P 1 − xi )(x P 2 − xi ) + (x P 2 − xi )(x P 3 − xi ) + (x P 3 − xi )(x P 1 − xi ) f x2 ,i

(x P 1 − xi )( y P 2 − y i ) + (x P 2 − xi )( y P 3 − y i ) + (x P 3 − xi )( y P 1 − y i ) 36  + ( y P 1 − y i )(x P 2 − xi ) + ( y P 2 − y i )(x P 3 − xi ) + ( y P 3 − y i )(x P 1 − xi ) f xy ,i

+ +

5  36 25

 ( y P 1 − y i )( y P 2 − y i ) + ( y P 2 − y i )( y P 3 − y i ) + ( y P 3 − y i )( y P 1 − y i ) f y 2 ,i

(x P 1 − xi )(x P 2 − xi )(x P 3 − xi ) f x3 ,i 324 25  (x P 1 − xi )(x P 2 − xi )( y P 3 − y i ) + ( y P 1 − y i )(x P 2 − xi )(x P 3 − xi ) + 324  + (x P 1 − xi )( y P 2 − y i )(x P 3 − xi ) f x2 y ,i 25  ( y P 1 − y i )( y P 2 − y i )(x P 3 − xi ) + (x P 1 − xi )( y P 2 − y i )( y P 3 − y i ) 324  + ( y P 1 − y i )(x P 2 − xi )( y P 3 − y i ) f xy 2 ,i

+

+

25 324

( y P 1 − y i )( y P 2 − y i )( y P 3 − y i ) f y 3 ,i .

(5.8)

We then arrive at the following result. Theorem 5.1. A PS5-spline s(x, y ) given by (4.4), with the coefficients c iv, j , j = 1, . . . , 10, corresponding to vertex V i chosen as





































c iv,1 = R i Q iv,1 , Q iv,1 , Q iv,1 ,

c iv,2 = R i Q iv,2 , Q iv,2 , Q iv,2 ,

c iv,3 = R i Q iv,3 , Q iv,3 , Q iv,3 ,

c iv,4 = R i Q iv,1 , Q iv,1 , Q iv,2 ,

c iv,5 = R i Q iv,1 , Q iv,1 , Q iv,3 ,

c iv,6 = R i Q iv,2 , Q iv,2 , Q iv,3 ,

c iv,7 = R i Q iv,2 , Q iv,2 , Q iv,1 ,

c iv,8 = R i Q iv,3 , Q iv,3 , Q iv,1 ,

c iv,9 = R i Q iv,3 , Q iv,3 , Q iv,2 ,





c iv,10 = R i Q iv,1 , Q iv,2 , Q iv,3 ,

(5.9)

satisfies interpolation problem (4.3a) at vertex V i . Note that the expressions in (5.9) only require to know the Cartesian coordinates of vertex V i and the three points Q iv, j , j = 1, 2, 3, of the corresponding PS5-triangle t i .

H. Speleers / Applied Numerical Mathematics 62 (2012) 620–635

629

Fig. 5. Schematic representation of the Bézier ordinates of a PS5-spline.

5.2. Interpolation at the split points We now show how to satisfy (4.3b). To compute the value of a PS5-spline at a split point Z k we make use of the Bernstein–Bézier representation. We consider the macro-triangle Tk ( V 1 , V 2 , V 3 ), as shown in Fig. 2, and we assume that the points indicated in the figure have the following barycentric coordinates:

V 1 = (1, 0, 0),

V 2 = (0, 1, 0),

R 12 = (λ12 , λ21 , 0),

V 3 = (0, 0, 1),

R 23 = (0, λ23 , λ32 ),

Z k = ( z1 , z2 , z3 ),

R 31 = (λ13 , 0, λ31 ).

(5.10)

The Bézier ordinates are schematically represented in Fig. 5. In [20] it is shown how these ordinates can be computed from the spline coefficients c iv, j and ckt in a stable way, using a sequence of convex combinations. We will briefly review the computation steps without going into the details. We will use the same notation as in [20] for the representation of all (intermediate) ordinates. The Bézier ordinates d1 , . . . , d16 in the neighbourhood of V 1 are completely determined by control polynomial T 1 (x, y ). These ordinates can be found through subdivision of T 1 (x, y ). A two-stage subdivision approach was proposed in [20]. In the first stage the Bernstein–Bézier form of the control polynomial is computed onto the triangle with the PS5-points V 1 , S 12 and S 13 as its three vertices. The Bézier ordinates of the control polynomial on this finer triangle are denoted by e 1, j , j = 1, . . . , 10, as illustrated in Fig. 6. In the second stage the values of d1 , . . . , d16 are computed from the ordinates e 1, j through subdivision. In a similar way, the values of Bézier ordinates d17 , . . . , d48 in the neighbourhood of vertices V 2 and V 3 can be computed from the intermediate ordinates e 2, j and e 3, j , j = 1, . . . , 10. Here, e 2, j are the Bézier ordinates of control polynomial T 2 (x, y ) restricted on the triangle with vertices ( V 2 , S 23 , S 21 ), and e 3, j are the Bézier ordinates of control polynomial T 3 (x, y ) restricted on the triangle with vertices ( V 3 , S 31 , S 32 ). The values of the Bézier ordinates d49 , . . . , d75 are obtained from the C 3 -smoothness conditions of the PS5-spline across the interior edges Z k − R 12 , Z k − R 23 and Z k − R 31 . The remaining Bézier ordinates d76 , . . . , d91 are found by subdivision of a single cubic polynomial p 3 defined on the triangle with vertices (4.8). This polynomial p 3 is defined by the following ten Bézier ordinates:

b300 = d7 , b003 = d39 ,

b210 = d˜ 12 , b102 = d˜ 44 ,

b201 = d˜ 14 , b012 = d˜ 46 ,

b030 = d23 ,

b021 = d˜ 28 ,

b120 = d˜ 30 ,

b111 = ckt .

(5.11)

The values of d7 , d˜ 12 and d˜ 14 are given by

d7 = z1 ( z1 e 1,1 + z2 e 1,4 + z3 e 1,5 ) + z2 ( z1 e 1,4 + z2 e 1,7 + z3 e 1,10 ) + z3 ( z1 e 1,5 + z2 e 1,10 + z3 e 1,8 ),

(5.12a)

d˜ 12 = z1 ( z1 e 1,4 + z2 e 1,7 + z3 e 1,10 ) + z2 ( z1 e 1,7 + z2 e 1,2 + z3 e 1,6 ) + z3 ( z1 e 1,10 + z2 e 1,6 + z3 e 1,9 ),

(5.12b)

d˜ 14 = z1 ( z1 e 1,5 + z2 e 1,10 + z3 e 1,8 ) + z2 ( z1 e 1,10 + z2 e 1,6 + z3 e 1,9 ) + z3 ( z1 e 1,8 + z2 e 1,9 + z3 e 1,3 ).

(5.12c)

630

H. Speleers / Applied Numerical Mathematics 62 (2012) 620–635

Fig. 6. A PS5-triangle t 1 ( Q 1v,1 , Q 1v,2 , Q 1v,3 ) of vertex V 1 containing the PS5-points V 1 , S 12 and S 13 on triangle T ( V 1 , V 2 , V 3 ), together with the schematic representation of the Bézier ordinates e 1, j , j = 1, . . . , 10, of the subdivided control polynomial T 1 (x, y ) onto the triangle with as vertices the points V 1 , S 12 and S 13 .

The values of d23 , d˜ 28 and d˜ 30 are computed from e 2, j , j = 1, . . . , 10, and their expressions are similar to (5.12). Likewise,

the expressions for d39 , d˜ 44 and d˜ 46 are based on e 3, j , j = 1, . . . , 10. For more details we refer to [20]. The value of PS5-spline s(x, y ) at split point Z k can be computed by evaluating the polynomial p 3 with Bézier ordinates (5.11) at Z k , i.e.,

s( Z k ) = d91 = ( z1 )3 d7 + 3( z1 )2 z2 d˜ 12 + 3( z1 )2 z3 d˜ 14 + ( z2 )3 d23 + 3( z2 )2 z3 d˜ 28

+ 3( z2 )2 z1 d˜ 30 + ( z3 )3 d39 + 3( z3 )2 z1 d˜ 44 + 3( z3 )2 z2 d˜ 46 + 6z1 z2 z3 ckt .

(5.13)

The following theorem follows directly from (5.13). Theorem 5.2. A PS5-spline s(x, y ) given by (4.4), with the coefficient ckt corresponding to triangle Tk ( V 1 , V 2 , V 3 ) ∈ chosen as

ckt =

1 6z1 z2 z3



gk − ( z1 )3 d7 − 3( z1 )2 z2 d˜ 12 − 3( z1 )2 z3 d˜ 14 − ( z2 )3 d23 − 3( z2 )2 z3 d˜ 28

 − 3( z2 )2 z1 d˜ 30 − ( z3 )3 d39 − 3( z3 )2 z1 d˜ 44 − 3( z3 )2 z2 d˜ 46 ,

(5.14)

satisfies interpolation problem (4.3b) at the split point Z k in triangle Tk . Combining Theorems 5.1 and 5.2 gives us the full Hermite interpolation scheme for PS5-splines in representation (4.4). In the next section we consider some alternative interpolation schemes requiring only data at the vertices in and still preserving the same approximation order. 6. Reduced PS5-spline interpolants 6.1. C 3 -smoothness across an interior edge Imposing additional local supersmoothness is a possible way to reduce the number of degrees of freedom [10]. We now derive the restriction on the spline coefficient ckt such that the PS5-spline will be C 3 -continuous across the interior edge Z k − V 1 inside macro-triangle Tk ( V 1 , V 2 , V 3 ), see Fig. 2. This condition is satisfied when the values of d3 , d6 , d7 , d8 , d11 , d12 , d13 , d14 , d15 , d52 , d55 , d72 , d75 , d76 , d84 and d85 , see Fig. 5, can be regarded as Bézier ordinates of a single cubic polynomial after subdivision. We consider the cubic polynomial with Bézier ordinates f abc , a + b + c = 3, defined on the triangle with vertices

P1 =

4 5

V1 +

1 5

Zk ,

P2 =

1 5

V1 +

3 5

V2 +

1 5

Zk ,

P3 =

1 5

V1 +

3 5

V3 +

1 5

Zk .

(6.1)

In case of the required additional C 3 -smoothness, it follows that the values of d3 , d6 , d7 , d8 , d11 , d12 , d13 , d14 , d15 , d52 , d55 , d72 , d75 , d76 , d84 and d85 can be computed from f abc using the formulae (2.4)–(2.5). This gives us new expressions for these

H. Speleers / Applied Numerical Mathematics 62 (2012) 620–635

631

di in terms of f abc . Comparing the new expressions with the expressions given in [20], we get the restriction on the spline coefficient ckt . We recall from [20] the following relations:

d3 = z1 e 1,1 + z2 e 1,4 + z3 e 1,5 , d˜ 6 =

(6.2a)

d6 − λ12 d3

= z1 e 1,4 + z2 e 1,7 + z3 e 1,10 , λ21 d8 − λ13 d3 d˜ 8 = = z1 e 1,5 + z2 e 1,10 + z3 e 1,8 , λ31 d11 − λ12 d6 − λ12 λ21 d˜ 6 ˜ d˜ 11 = = z1 e 1,7 + z2 e 1,2 + z3 e 1,6 , λ21 2 d15 − λ13 d8 − λ13 λ31 d˜ 8 ˜ d˜ 15 = = z1 e 1,8 + z2 e 1,9 + z3 e 1,3 , λ31 2 ˜ d12 − λ12 d7 − λ21 z1 d˜ 6 − λ21 z2 d˜ 11 d˜ 13 = = z1 e 1,10 + z2 e 1,6 + z3 e 1,9 . λ21 z3

(6.2b) (6.2c) (6.2d) (6.2e)

(6.2f)

In case of C 3 -smoothness across the interior edge Z k − V 1 , it is easy to verify that we must have

d˜ 6 = f 210 ,

d3 = f 300 ,

d˜ 8 = f 201 ,

˜

d˜ 11 = f 120 ,

˜

d˜ 15 = f 102 ,

d˜ 13 = f 111 .

(6.3)

We also obtain from [20] that

d52 − λ12 d11

˜ ˜ = λ12 (λ12 d˜ 6 + λ21 d˜ 11 ) + λ21 (λ12 d˜ 31 + λ21 d˜ 24 ), λ21 d55 − λ12 d12 ˜ ˜ d˜ 56 = = λ12 ( z1 d˜ 6 + z2 d˜ 11 + z3 d˜ 13 ) + λ21 ( z1 d˜ 31 + z2 d˜ 24 + z3 d˜ 29 ), λ21 d72 − λ13 d15 ˜ ˜ d˜ 71 = = λ13 (λ13 d˜ 8 + λ31 d˜ 15 ) + λ31 (λ13 d˜ 43 + λ31 d˜ 38 ), λ31 d75 − λ13 d14 ˜ ˜ d˜ 74 = = λ13 ( z1 d˜ 8 + z2 d˜ 13 + z3 d˜ 15 ) + λ31 ( z1 d˜ 43 + z2 d˜ 45 + z3 d˜ 38 ), λ31 d˜ 53 =

(6.4a) (6.4b) (6.4c) (6.4d)

where

d˜ 24 =

˜

d˜ 31 = d˜ 29 =

d24 − λ21 d19

λ12

= z1 e 2,8 + z2 e 2,5 + z3 e 2,10 ,

d31 − λ21 d24 − λ21 λ12 d˜ 24

λ12 2

= z1 e 2,3 + z2 e 2,8 + z3 e 2,9 ,

˜

d30 − λ21 d23 − λ12 z1 d˜ 31 − λ12 z2 d˜ 24

λ12 z3

= z1 e 2,9 + z2 e 2,10 + z3 e 2,6 ,

(6.5a) (6.5b)

(6.5c)

and

d˜ 38 =

d38 − λ31 d35

˜

d43 − λ31 d38 − λ31 λ13 d˜ 38

d˜ 43 = d˜ 45 =

λ13

= z1 e 3,7 + z2 e 3,10 + z3 e 3,4 ,

λ13 2 ˜

= z1 e 3,2 + z2 e 3,6 + z3 e 3,7 ,

d44 − λ31 d39 − λ13 z1 d˜ 43 − λ13 z3 d˜ 38

λ13 z2

= z1 e 3,6 + z2 e 3,9 + z3 e 3,10 .

(6.6a) (6.6b)

(6.6c)

In case of C 3 -smoothness across Z k − V 1 , we find that

d˜ 53 = λ12 (λ12 f 210 + λ21 f 120 ) + λ21 (λ12 f 120 + λ21 f 030 ),

(6.7a)

d˜ 56 = λ12 ( z1 f 210 + z2 f 120 + z3 f 111 ) + λ21 ( z1 f 120 + z2 f 030 + z3 f 021 ),

(6.7b)

d˜ 71 = λ13 (λ13 f 201 + λ31 f 102 ) + λ31 (λ13 f 102 + λ31 f 003 ),

(6.7c)

d˜ 74 = λ13 ( z1 f 201 + z2 f 111 + z3 f 102 ) + λ31 ( z1 f 102 + z2 f 012 + z3 f 003 ).

(6.7d)

632

H. Speleers / Applied Numerical Mathematics 62 (2012) 620–635

Combining (6.3), (6.4) and (6.7), we find the following necessary conditions for the additional C 3 -smoothness

˜ λ12 d˜ 31 + λ21 d˜ 24 = λ12 f 120 + λ21 f 030 , ˜ z1 d˜ 31 + z2 d˜ 24 + z3 d˜ 29 = z1 f 120 + z2 f 030 + z3 f 021 , ˜ λ13 d˜ 43 + λ31 d˜ 38 = λ13 f 102 + λ31 f 003 , ˜ z1 d˜ 43 + z2 d˜ 45 + z3 d˜ 38 = z1 f 102 + z2 f 012 + z3 f 003 .

(6.8a) (6.8b) (6.8c) (6.8d)

We can reformulate (6.8) as

λ12 ˜˜ ˜ (d31 − d˜ 11 ), λ21   z1 z2 λ12 ˜ ˜ f 021 = d˜ 29 + − (d˜ 31 − d˜ 11 ), z3 z3 λ21 λ13 ˜˜ ˜ (d43 − d˜ 15 ), f 003 = d˜ 38 + λ31   z1 z3 λ13 ˜˜ ˜ ˜ (d43 − d˜ 15 ). − f 012 = d45 + z2 z2 λ31 f 030 = d˜ 24 +

(6.9a) (6.9b) (6.9c) (6.9d)

Finally, we deduce from [20] that

d˜ 77 =

d76 − λ12 d13

λ21

˜ ˜ = z1 ( z1 d˜ 6 + z2 d˜ 11 + z3 d˜ 13 ) + z2 ( z1 d˜ 31 + z2 d˜ 24 + z3 d˜ 29 ) + z3 ckt .

(6.10)

The additional C 3 -smoothness across Z k − V 1 implies that

d˜ 77 = z1 ( z1 f 210 + z2 f 120 + z3 f 111 ) + z2 ( z1 f 120 + z2 f 030 + z3 f 021 ) + z3 ( z1 f 111 + z2 f 021 + z3 f 012 ),

(6.11)

resulting in the following condition on the spline coefficient ckt :

ckt = z1 f 111 + z2 f 021 + z3 f 012 .

(6.12)

Filling (6.2)–(6.3), (6.5)–(6.6) and (6.9) into (6.12), we get an expression for in terms of e i , j , i = 1, 2, 3 and j = 1, . . . , 10. Since the ordinates e i , j can be computed through subdivision from the three control polynomials T i (x, y ) corresponding to the three vertices V i , i = 1, 2, 3, we have an expression for ckt in terms of the vertex coefficients c iv, j , i = 1, 2, 3 and j = 1, . . . , 10. A similar procedure can be used to derive the restrictions on ckt for additional C 3 -smoothness across the other interior edges Z k − V 2 or Z k − V 3 . ckt

6.2. Reduced PS5-splines In this section we consider some reduced PS5-spline spaces that preserve the approximation order. The B-spline coefficients ckt , k = 1, . . . , nt , in representation (4.4) are determined in a specific way such that the obtained PS5-splines are still able to reproduce quintic polynomials. A first approach is requiring additional C 3 -smoothness across an interior edge for all macro-triangles in the triangulation. Such a reduced PS5-spline space was described in [1,10]. For each macro-triangle Tk ( V 1 , V 2 , V 3 ) in , there are three possible choices to impose the additional smoothness, i.e., across the edge Z k − V 1 , Z k − V 2 or Z k − V 3 . In case of Z k − V 1 , we get the following result. Theorem 6.1. PS5-splines s(x, y ) given by (4.4), with the coefficient ckt corresponding to triangle Tk ( V 1 , V 2 , V 3 ) ∈ chosen as

ckt = z1 d˜ 13 + z2 d˜ 29 + z3 d˜ 45 +

z2 z3



z1 − z2

   λ12 ˜˜ z3 λ13 ˜˜ ˜ ˜ (d31 − d˜ 11 ) + (d43 − d˜ 15 ), z1 − z3 λ21 z2 λ31

(6.13)

can reproduce quintic polynomials. Proof. From (6.3), (6.9) and (6.12), it follows that the spline s(x, y ) is C 3 -continuous across interior edge Z k − V 1 . From [1] we know that such a spline is able to reproduce quintic polynomials. 2 The next theorem provides a symmetric condition on the coefficient ckt . This reduced spline space was developed in [20].

H. Speleers / Applied Numerical Mathematics 62 (2012) 620–635

633

Theorem 6.2. PS5-splines s(x, y ) given by (4.4), with the coefficient ckt corresponding to triangle Tk ( V 1 , V 2 , V 3 ) ∈ chosen as

ckt = z1 ( z1 e 1,10 + z2 e 1,6 + z3 e 1,9 ) + z2 ( z1 e 2,9 + z2 e 2,10 + z3 e 2,6 ) + z3 ( z1 e 3,6 + z2 e 3,9 + z3 e 3,10 ),

(6.14)

can reproduce quintic polynomials. From [20] we know that the values e i , j , i = 1, 2, 3 and j = 1, . . . , 10, are computed as convex combinations of the B-spline coefficients c iv, j corresponding to the three vertices V i , i = 1, 2, 3, of the macro-triangle Tk . It follows that (6.14) is also a convex combination of these B-spline coefficients c iv, j . Because of the local support and the convex partition of unity property (4.5) of the Powell–Sabin B-splines, each value of a reduced PS5-spline s(x, y ) constructed by Theorem 6.2 on the triangle Tk is a local convex combination of its B-spline coefficients c iv, j , i = 1, 2, 3 and j = 1, . . . , 10. We now consider reduced PS5-splines defined on a uniform PS-triangulation ∗ , where the barycentric coordinates in (5.10) satisfy

z1 = z2 = z3 = 1/3 and

λ12 = λ21 = λ23 = λ32 = λ31 = λ13 = 1/2

(6.15)

on all triangles. Theorem 6.3. Let ∗ be a uniform PS-refinement of a uniform triangulation , satisfying (6.15), then the reduced PS5-splines described in Theorems 6.1 and 6.2 are equivalent on ∗ . Moreover, these splines are locally C 3 -continuous on each macro-triangle in . Proof. The first three terms in (6.13) are equivalent to (6.14). The remaining terms in (6.13) are zero because of the uniformity property (6.15) of the PS-triangulation. The local C 3 -continuity on each macro-triangle follows from the symmetry of formula (6.14). 2 The reduced spline space described in Theorem 6.3 was already considered in [7] in the case of uniform three-directional meshes. Theorem 5.1 can be used to get an interpolation scheme based on the reduced PS5-splines described in Theorems 6.1 and 6.2. These reduced spline interpolants have the advantage that only data are needed located at the vertices of the triangulation . 7. A numerical example We consider Franke’s test function [6],

f (x, y ) =

3 4

 exp −

+

1 2

(9x − 2)2 + (9 y − 2)2

 exp −

4

 + 2

(9x − 7) + (9 y − 3) 2

4

3 4

  (9x + 1)2 9 y + 1 exp − −

−

49

1 5

10





exp −(9x − 4)2 − (9 y − 7)2 ,

(7.1)

n

on the domain Ω = [0, 1] × [0, 1]. Let { V i }i =v 1 be the set of 25 data points given in Fig. 7. The figure shows the corresponding Delaunay triangulation with a particular PS-refinement. The dimension of the full PS5-spline space on this mesh is 282, and the dimension of the reduced spline space is 250. The data f xa yb ,l and gm in (4.3) are sampled from the function f (x, y ). To measure the accuracy of spline fit s(x, y ), we consider the maximum error and the mean error on a 400 × 400 uniform grid G on Ω , i.e.,





E max = max  f (x, y ) − s(x, y ), (x, y )∈G



E mean =

(x, y )∈G | f (x, y ) − s(x, y )|

400 × 400

.

The PS5-spline interpolant, defined by Theorems 5.1 and 5.2 on the mesh shown in Fig. 7, has a maximum error equal to 0.01018, and its mean error is 0.00040. The reduced PS5-spline interpolant, defined by Theorems 5.1 and 6.2, has a maximum error of 0.01503 and a mean error of 0.00060. This reduced PS5-spline interpolant is depicted in Fig. 8. In order to test the approximation power of the methods numerically, we computed PS5-splines interpolating (7.1) on several uniform three-directional meshes. The maximum and mean errors evaluated on a 400 × 400 uniform grid on Ω are shown in Table 1 for the full PS5-spline interpolants defined by Theorems 5.1 and 5.2, and for the reduced PS5-spline interpolants defined by Theorems 5.1 and 6.2. Note that these reduced PS5-splines also satisfy (6.13) as stated in Theorem 6.3. Because PS5-splines are quintic splines, the expected optimal approximation order is six, meaning that the error converges asymptotically like O (h6 ), with h the mesh size of the triangulation. By halving the mesh size in each step, we may then expect that the error decreases asymptotically by a factor of about 64. In Table 1 we see that the mean error of the full PS5-spline interpolants decreases successively by factors 20.8, 59.1, 49.1, 57.6 and 62.1 for this numerical example. The mean error of the reduced PS5-spline interpolants decreases successively by factors 17.4, 57.5, 51.8, 58.3 and 62.2.

634

H. Speleers / Applied Numerical Mathematics 62 (2012) 620–635

Fig. 7. A triangulation with PS-refinement.

Fig. 8. A reduced PS5-spline (with the triangular mesh lines) and its contour plot. The spline is defined on the mesh in Fig. 7, and it interpolates the function (7.1) at the vertices.

Table 1 The degrees of freedom, the maximum errors and mean errors of the full and reduced PS5-spline interpolants to function (7.1) defined on uniform threedirectional meshes consisting of n v vertices and nt triangles. Full PS5

Reduced PS5

nv

nt

dof

E max

E mean

dof

E max

E mean

9 25 81 289 1089 4225

8 32 128 512 2048 8192

98 282 938 3402 12 938 50 442

1.328e−01 1.757e−02 3.675e−04 1.569e−05 3.763e−07 5.945e−09

1.254e−02 6.025e−04 1.020e−05 2.076e−07 3.602e−09 5.801e−11

90 250 810 2890 10 890 42 250

1.328e−01 2.460e−02 6.080e−04 1.900e−05 4.291e−07 7.010e−09

1.726e−02 9.893e−04 1.722e−05 3.325e−07 5.708e−09 9.173e−11

8. Concluding remarks In this paper we have provided local Hermite interpolation rules for quintic Powell–Sabin splines represented in the normalized B-spline basis. The general interpolation scheme uses data at the vertices of the triangulation and at the

H. Speleers / Applied Numerical Mathematics 62 (2012) 620–635

635

Powell–Sabin split points. We have also discussed interpolation based on several reduced spline spaces. These reduced spline interpolants only need data at the vertices of the triangulation. Instead of using the full PS5-spline basis, the reduced PS5-splines can, of course, also be represented in a normalized basis with a reduced number of basis functions. Such a basis function can be defined by setting in the spline representation (4.4), with ckt chosen as in Theorem 6.1 or 6.2, a single coefficient c iv, j to be one and all other coefficients c iv, j to be zero. These basis functions have a local support, they form a partition of unity, and in the case of Theorem 6.2, they are also all nonnegative. The presented Hermite interpolation schemes make use of derivative information (up to order three) at the vertices. When no derivative information is available, these interpolation schemes can be incorporated into a two-stage method. In the first stage the needed derivative information is then locally estimated from the data, for example by doing local least squares with polynomials or radial basis functions. A good interpolation scheme must preserve locality and approximation order, and it is expected to have some shape fidelity. Acknowledgements Hendrik Speleers is a Postdoctoral Fellow of the Research Foundation Flanders (Belgium). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

P. Alfeld, L.L. Schumaker, Smooth macro-elements based on Powell–Sabin triangle splits, Adv. Comput. Math. 16 (2002) 29–46. L. De Lathauwer, B. De Moor, J. Vandewalle, A multilinear singular value decomposition, SIAM J. Matrix Anal. Appl. 21 (2000) 1253–1278. P. Dierckx, On calculating normalized Powell–Sabin B-splines, Comput. Aided Geom. Design 15 (1997) 61–78. P. Dierckx, S. Van Leemput, T. Vermeire, Algorithms for surface fitting using Powell–Sabin splines, IMA J. Numer. Anal. 12 (1992) 271–299. G. Farin, Triangular Bernstein–Bézier patches, Comput. Aided Geom. Design 3 (1986) 83–127. R. Franke, Scattered data interpolation: tests of some methods, Math. Comp. 38 (1982) 181–200. M. Laghchim-Lahlou, P. Sablonnière, C r -finite elements of Powell–Sabin type on the three direction mesh, Adv. Comput. Math. 6 (1996) 191–206. M.J. Lai, On C 2 quintic spline functions over triangulations of Powell–Sabin’s type, J. Comput. Appl. Math. 73 (1996) 135–155. M.J. Lai, L.L. Schumaker, Macro-elements and stable local bases for splines on Powell–Sabin triangulations, Math. Comp. 72 (2003) 335–354. M.J. Lai, L.L. Schumaker, Spline Functions on Triangulations, Encyclopedia of Mathematics and its Applications, vol. 110, Cambridge University Press, 2007. C. Manni, P. Sablonnière, Quadratic spline quasi-interpolants on Powell–Sabin partitions, Adv. Comput. Math. 26 (2007) 283–304. B. Mulansky, J.W. Schmidt, Powell–Sabin splines in range restricted interpolation of scattered data, Computing 53 (1994) 137–154. M.J.D. Powell, M.A. Sabin, Piecewise quadratic approximations on triangles, ACM Trans. Math. Softw. 3 (1977) 316–325. L. Ramshaw, Blossoming: a connect-the-dots approach to splines, Tech. Rep. 19, Digital Systems Research Center, 1987. P. Sablonnière, Composite finite elements of class C 2 , in: C.K. Chui, L.L. Schumaker, F.I. Utreras (Eds.), Topics in Multivariate Approximation, Academic Press, 1987, pp. 207–217. P. Sablonnière, Error bounds for Hermite interpolation by quadratic splines on an α -triangulation, IMA J. Numer. Anal. 7 (1987) 495–508. D. Sbibih, A. Serghini, A. Tijini, Polar forms and quadratic spline quasi-interpolants on Powell–Sabin partitions, Appl. Numer. Math. 59 (2009) 938–958. L.L. Schumaker, H. Speleers, Nonnegativity preserving macro-element interpolation of scattered data, Comput. Aided Geom. Design 27 (2010) 245–261. H.P. Seidel, An introduction to polar forms, IEEE Comp. Graph. Appl. 13 (1993) 38–46. H. Speleers, A normalized basis for quintic Powell–Sabin splines, Comput. Aided Geom. Design 27 (2010) 438–457. K. Willemans, P. Dierckx, Surface fitting using convex Powell–Sabin splines, J. Comput. Appl. Math. 56 (1994) 263–282. K. Willemans, P. Dierckx, Smoothing scattered data with a monotone Powell–Sabin spline surface, Numer. Algorithms 12 (1996) 215–232.