Construction and analysis of cubic Powell–Sabin B-splines

Computer Aided Geometric Design 57 (2017) 1–22

Contents lists available at ScienceDirect

Computer Aided Geometric Design www.elsevier.com/locate/cagd

Construction and analysis of cubic Powell–Sabin B-splines Jan Grošelj a,∗ , Hendrik Speleers b a b

Faculty of Mathematics and Physics, University of Ljubljana, Jadranska 19, 1000 Ljubljana, Slovenia Department of Mathematics, University of Rome ‘Tor Vergata’, Via della Ricerca Scientifica, 00133 Rome, Italy

a r t i c l e

i n f o

Article history: Available online 24 May 2017 Keywords: C 1 cubic splines Powell–Sabin splines Normalized B-splines Stable basis Control structure

a b s t r a c t We consider a new B-spline representation for the space of C 1 cubic splines defined on a triangulation with a Powell–Sabin refinement. The construction is based on lifting particular triangles and line segments from the domain. We prove that the B-splines form a locally supported stable basis and a convex partition of unity. Furthermore, we provide explicit expressions for the B-spline coefficients of any element of the cubic spline space and show how to compute the Bernstein–Bézier form of such a spline in a stable way. The B-spline representation induces a natural control structure that is useful for geometric modelling. Finally, we explore how classical C 1 quadratic Powell–Sabin splines and C 1 cubic Clough–Tocher splines can be expressed in the new B-spline representation. © 2017 Elsevier B.V. All rights reserved.

1. Introduction Bivariate spline functions are a mature and well-researched field in numerical analysis. A bivariate spline on a triangulation is a piecewise function consisting of bivariate polynomials of a particular degree on every triangle of the triangulation that are joined together with a certain order of smoothness across the edges of neighbouring triangles. A general construction of splines of a low degree with respect to smoothness is difficult because it depends on the geometry of the underlying triangulation. Although many strategies have been developed to overcome this problem, one of the most popular and successful approaches is to refine the given triangulation into a finer triangulation by applying a certain split on every triangle. The most common splitting techniques are the Clough–Tocher and Powell–Sabin splits that refine every triangle into three and six smaller triangles, respectively (see, e.g., Lai and Schumaker, 2007). The construction of C 1 cubic splines on Clough–Tocher triangulations was pioneered by Clough and Tocher (1965) in the context of finite elements. Due to the local nature of their definition, these splines were also successfully applied to scattered data interpolation (see Farin, 1985; Mann, 1999). Several polynomial splines of higher degree and smoothness were developed on Clough–Tocher triangulations (see Lai and Schumaker, 2007, for an overview), but their constructions were made possible with the help of additional smoothness constraints resulting in local super-smoothness. The story of splines on Powell–Sabin triangulations is similar. Powell and Sabin (1977) proposed C 1 quadratic splines for contour plotting, and later families of polynomial splines of higher degree and smoothness were studied (see, again, Lai and Schumaker, 2007). However, Powell–Sabin splines have become increasingly interesting since the development of the (normalized) B-spline representation of quadratic splines by Dierckx (1997). This representation has shown to be fruitful in many applications, such as surface modelling and compression (Dierckx, 1997; Maes and Bultheel, 2006; Vanraes et al., 2004), interpolation and approximation (Manni and Sablonnière, 2007; Sbibih et al., 2009), and solving partial dif-

*

Corresponding author. E-mail addresses: [email protected] (J. Grošelj), [email protected] (H. Speleers).

http://dx.doi.org/10.1016/j.cagd.2017.05.003 0167-8396/© 2017 Elsevier B.V. All rights reserved.

2

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

ferential equations numerically (Speleers et al., 2006, 2012). The B-spline construction has been adapted in several ways to Powell–Sabin splines of higher degree; in particular, C 1 cubic splines (Grošelj and Krajnc, 2016; Lamnii et al., 2014; Speleers, 2015) and splines of arbitrary smoothness (Grošelj, 2016; Speleers, 2010a, 2013a), which have certain additional local super-smoothness. The C 1 quadratic construction has also been extended to the trivariate setting (Sbibih et al., 2012) and the general multivariate setting (Speleers, 2013b). Other quadratic spline representations on triangulations endowed with a split into 12 subtriangles can be found in Cohen et al. (2013); Lamnii et al. (2015). The presence of local super-smoothness in spline spaces helps to reduce the number of degrees of freedom and to simplify the construction of basis functions. However, there is also a drawback: it complicates (and often prohibits) the definition of sequences of nested spaces. Nestedness is an important feature in computer aided geometric design and isogeometric analysis. For example, it allows the construction of subdivision schemes, as explored by Vanraes et al. (2004); Windmolders and Dierckx (1999) for C 1 quadratic Powell–Sabin splines. In this work, we consider a B-spline representation of pure C 1 cubic splines on a Powell–Sabin triangulation with no additional local super-smoothness. A different B-spline representation for these splines has already been provided by Grošelj and Krajnc (2016) but is based on B-splines that are non-negative only under certain restrictions on the geometry of the underlying Powell–Sabin triangulation. The B-splines presented herein are non-negative regardless of the triangulation and are a generalization of the B-splines that were developed by Speleers (2015) for a particular subspace of C 1 cubic splines with some additional smoothness properties in the interior of the triangles. We show that the cubic B-splines we introduce have many desirable properties from the perspective of computer aided geometric design and approximation theory, i.e., they are local, stable, and form a convex partition of unity. They are also attractive because they are able to represent both C 1 quadratic Powell–Sabin splines and C 1 cubic Clough–Tocher splines in a unified context. The paper is organized as follows. In Section 2, we first review some preliminaries on bivariate polynomials on triangles and splines on triangulations. Then, we provide a minimal determining set for the space of C 1 cubic Powell–Sabin splines and discuss different interpolation problems that characterize the elements of the space. Section 3 covers a B-spline construction that is based on lifting particular triangles and line segments from the domain. In Section 4, we examine the corresponding B-spline representation. We prove that the B-splines form a locally supported stable basis for C 1 cubic Powell–Sabin splines and a convex partition of unity. Furthermore, we give explicit expressions for the Greville points and the coefficients of the B-spline representation, and we present a stable computation of the Bernstein–Bézier form. In Section 5, we propose a control structure for splines in the B-spline representation. Section 6 explores the relationship of the B-spline representation with other macro-element spaces; in particular, C 1 quadratic Powell–Sabin splines and C 1 cubic Clough–Tocher splines. In Section 7, we end with a few concluding remarks. 2. C 1 cubic splines on Powell–Sabin triangulations In this section, we first review a few basic definitions and properties of bivariate polynomials on triangles and splines on triangulations. We then introduce the space of C 1 cubic splines on Powell–Sabin triangulations and discuss the characterization of its elements. 2.1. Bernstein–Bézier form of bivariate polynomials and splines Let T =  V 1 , V 2 , V 3  be a non-degenerate triangle in R2 with vertices V 1 , V 2 , V 3 . Every point P = (x, y ) ∈ R2 can be uniquely described by a triplet of barycentric coordinates (τ1 , τ2 , τ3 ) with respect to T . This triplet is uniquely determined by the conditions

P = τ1 V 1 + τ2 V 2 + τ3 V 3 ,

τ1 + τ2 + τ3 = 1.

The values τ1 , τ2 , τ3 are all non-negative whenever P belongs to T . The set of Bernstein basis polynomials { B di, j ,k }i + j +k=d of degree d with respect to T is defined by

B di, j ,k (x, y ) = B di, j ,k ( P ) =

d! i ! j !k!

τ1i τ2j τ3k , i + j + k = d.

These basis polynomials form a convex partition of unity on T . We say that the multi-indices (i , j , k) and (i  , j  , k ) are at distance  if

|i − i  | + | j − j  | + |k − k | = 2.

(1)

Let Pd be the space of bivariate polynomials of total degree less than or equal to d. Every polynomial p ∈ Pd can be uniquely expressed in terms of the Bernstein basis polynomials of degree d, i.e.,

p (x, y ) =



b i , j ,k B di, j ,k (x, y ),

(2)

i + j +k=d

for some values b i , j ,k ∈ R called the Bézier ordinates of p. The expression (2) is referred to as the Bernstein–Bézier representation of the polynomial p. There is a natural association between the ordinate b i , j ,k and the Bézier domain point j

B i , j ,k ∈ R2 specified by the barycentric coordinates ( di , d , dk ) with respect to T . Since

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22



B i , j ,k B di, j ,k ( P ) = P ,

3

(3)

i + j +k=d

every point on the graph of p above T can be described as a convex combination of the control points ( B i , j ,k , b i , j ,k ), and every such point can be evaluated in a stable way as a series of simple convex combinations by means of the de Casteljau algorithm. Moreover, the control points can be organized into a triangular net, known as the Bézier control net, by connecting the neighbouring control points at distance 1 according to relation (1). This net provides a rough approximation of the graph and is tangent to the graph surface at the domain triangle’s vertices. Furthermore, there exists a constant K depending on d such that

 p ∞,T ≤ max |bi , j ,k | ≤ K  p ∞,T ,

(4)

i + j +k=d

which ensures that the Bernstein–Bézier representation (2) is stable. A comprehensive overview of the theory on Bernstein– Bézier representations can be found in, e.g., Farin (1986); Lai and Schumaker (2007). The Bézier ordinates b i , j ,k of p can be neatly expressed as values of the unique, symmetric, multi-affine polynomial B p : (R2 )d → R with the property

B p ( P [d]) = p ( P ). Here, P [k] denotes the repetition of P k times. The function B p is called the blossom (or polar form) of p, and

b i , j ,k = B p ( V 1 [i ], V 2 [ j ], V 3 [k]).

(5)

The blossom value B p ( P 1 , . . . , P d ) can be efficiently calculated with the help of the multi-affine de Casteljau algorithm, i.e., the de Casteljau algorithm that uses the barycentric coordinates of the point P k at the k-th step. This also justifies the derivation formula

D vk · · · D v 1 p( P ) =

d! B p ( P [d − k], v 1 , . . . , v k ). (d − k)!

(6)

We refer to Seidel (1993); Speleers (2011) for more insight into blossoming of triangular polynomials. A triangulation  of a polygonal domain  in R2 is a finite set of triangles {Tn }n such that any non-empty intersection of two triangles Tn and Tn is either a common vertex or a common edge. The space of C r splines of degree d on  is given by

  Sdr () = s ∈ C r () : s|T ∈ Pd , T ∈  . The Bernstein–Bézier form of s ∈ Sdr () is a notion for the Bernstein–Bézier representation of the restrictions of s to every triangle of . Let us denote the union of all the Bézier domain points by D d (). A subset  ⊆ D d () is a determining set for Sdr () if setting the Bézier ordinates of s ∈ Sdr () corresponding to the points in  to zero implies that s ≡ 0. A determining set M with the smallest possible cardinality is a minimal determining set. It characterizes the spline space since every s ∈ Sdr () can be uniquely described by the Bézier ordinates corresponding to the points in M. All the other Bézier ordinates can be expressed using the continuity conditions. In particular, a well-known requirement for a C 1 joint between two polynomial patches of degree d with Bézier ordinates b i , j ,k and  b i , j ,k on triangles T =  V 1 , V 2 , V 3  and T =  V 2 , V 3 , V 4  along the common edge  V 2 , V 3  is

 b j ,k,0 = b0, j ,k ,

j + k = d,  b j ,k,1 = σ1 b1, j ,k + σ2 b0, j +1,k + σ3 b0, j ,k+1 ,

(7a) j + k = d − 1,

(7b)

where (σ1 , σ2 , σ3 ) are the barycentric coordinates of V 4 with respect to T . More details can be found in, e.g., Farin (1986); Lai and Schumaker (2007). 2.2. Space of C 1 cubic splines on a Powell–Sabin refined triangulation A Powell–Sabin refinement PS of  is a refinement of  obtained by splitting every triangle of  into six smaller triangles in the following way (Powell and Sabin, 1977). 1. For each triangle Tn =  V i , V j , V k  ∈ , select a triangle split point Z n inside Tn and connect it to the vertices V i , V j, V k. 2. For each pair of triangles Tn , Tn ∈  with a common edge  V i , V j , connect the triangle split points Z n and Z n . The points Z n and Z n must be chosen so that the line segment  Z n , Z n  intersects the interior of the edge  V i , V j . The point of intersection is an edge split point and is denoted by R i , j . 3. For each boundary edge  V i , V j  of , the edge split point R i , j can be an arbitrary point in the interior of the edge. Connect it to the triangle split point Z n of the triangle Tn ∈  that contains  V i , V j .

4

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

Fig. 1. A Powell–Sabin refinement of a triangulation  with domain points D 3 (PS ). The coloured points represent a minimal determining set for S13 (PS ).

An example of a Powell–Sabin refinement is shown in Fig. 1. The conditions in the construction of PS can always be satisfied, e.g., by selecting the triangle split points to be the incentres of the triangles in . Hereafter, the barycentric coordinates of the triangle split point Z n with respect to Tn =  V i , V j , V k  will be denoted by (ξn,i , ξn, j , ξn,k ). The barycentric coordinates of R i , j can be described by (λi , j , λ j ,i ) with respect to the edge  V i , V j  or by (μn,n , μn ,n ) with respect to the line segment  Z n , Z n  if the edge  V i , V j  is shared by the triangles Tn , Tn ∈ . Clearly, all these coordinates are positive. In this paper, we focus on the spline space S13 (PS ). The dimension of S13 (PS ) is given by (see Speleers, 2015; Grošelj and Krajnc, 2016)

dim(S13 (PS )) = 3|V | + 4|E | = 3|V | + 6|| + 2|∂ E |, where V is the set of vertices, E the set of edges, and ∂ E the set of boundary edges of . This dimension formula is confirmed by the following proposition. Proposition 1. Let M =

| V |

i =1

M i be a subset of D 3 (PS ) constructed as follows.

+ 13 R i , j , 23 V i + 13 Z n to M i . Z n to M i . Additionally, add the point 31 V i + 23 R i , j if

1. For each vertex V i of , select a triangle  V i , R i , j , Z n  ∈ PS , and add the points V i , 2. For each triangle  V i , R i , j , Z n  of PS , add the point  V i , R i , j  is a boundary edge of PS .

1 V 3 i

1 3

+ R i, j +

1 3

2 V 3 i

Then, M is a minimal determining set for S13 (PS ). Proof. Given the Bézier ordinates of s ∈ S13 (PS ) corresponding to the domain points in M, we need to show that the Bézier ordinates of s corresponding to the points in D 3 (PS ) \ M are uniquely determined. Consider a vertex V i of . Since M i contains V i and two points at distance 1 from V i (according to relation (1)), all the Bézier ordinates of s corresponding to the domain points at distance 1 from V i are determined by the C 1 continuity conditions at V i . Furthermore, since M i contains 13 V i + 13 R i , j + 13 Z n for every triangle  V i , R i , j , Z n  of PS with a vertex at V i , all the Bézier ordinates of s corresponding to the domain points at distance 2 from V i are determined by the C 1 continuity conditions across the edges of PS connecting V i with triangle and edge split points. The only exception is the ordinate corresponding to a boundary edge  V i , R i , j , but 13 V i + 23 R i , j is explicitly added to M i . Finally, the Bézier ordinates corresponding to the domain points at distance 3 from V i are determined by the C 1 continuity conditions at the triangle split points and across the edges of PS connecting edge and triangle split points since all the other Bézier ordinates are already fixed. 2

The proof of Proposition 1 shows that the Bézier ordinate of s ∈ S13 (PS ) corresponding to a domain point in D 3 (PS ) \ M depends only on the Bézier ordinates associated with the points in M that are in a close neighbourhood of the domain point in D 3 (PS ) \ M. Also, verifying that the value of such a Bézier ordinate can be calculated in a stable way is straightforward. We say that M is local and stable (Lai and Schumaker, 2007, Definition 5.16). Since S13 (PS ) clearly contains P3 , it follows that S13 (PS ) has an optimal approximation order (Lai and Schumaker, 2007, Section 5.7). Fig. 1 shows an example of a minimal determining set for S13 (PS ) as described in Proposition 1.

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

5

Fig. 2. Schematic representations of possible interpolation problems for S13 (PS ) on a single macro-triangle.

It is common to characterize a macro-element space with a Hermite interpolation problem that uniquely determines its elements (see, e.g., Lai and Schumaker, 2007). There are many ways to do that for S13 (PS ). Suppose that the values and the first order derivatives are prescribed at the vertices of . Then, the Bézier ordinates of a spline corresponding to the domain points at most at distance 1 from the vertices of  can be calculated. To determine all the other ordinates, one of the following options can be used. 1. Prescribe a value and a directional derivative (in a direction not parallel to the edge) at every edge of PS connecting a vertex of  with an edge split point (see Fig. 2, left). 2. Prescribe a value and a directional derivative (in a direction not parallel to the edge) at every edge of PS connecting a vertex of  with a triangle split point and the value and the derivative (in the direction parallel to the edge) at the edge split point of every boundary edge of  (see Fig. 2, middle). 3. Prescribe the value and the first order derivatives at every triangle split point of PS , a directional derivative (in a direction not parallel to the edge) at every edge of PS connecting a vertex of  with a triangle split point, and the value and the derivative (in the direction parallel to the edge) at the edge split point of every boundary edge of  (see Fig. 2, right). We refer to Grošelj and Krajnc (2016, Section 3) for more details on this specific interpolation problem. 3. Construction of a B-spline representation In this section, we introduce a B-spline representation of s ∈ S13 (PS ),

s(x, y ) =

|V |  3 

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

i =1 r =1

| PS | 

t t cm Bm (x, y ) +

m =1

|∂ EPS |

c e B e (x, y ).

(8)

=1

There are three kinds of basis functions in this representation. The B-splines B iv,r , r = 1, 2, 3, are associated with the vertices t V i ∈ V , the B-splines B m are associated with the triangles Tm ∈ PS , and the B-splines B e are associated with the boundary edges ε ∈ ∂ EPS . We show how to construct the B-splines so that they are locally supported and form a convex partition of unity. In contrast to the approach considered by Grošelj and Krajnc (2016, Section 4), we do not require any restrictions on the geometry of the underlying Powell–Sabin triangulation. 3.1. A geometric approach to B-splines The entire construction of the B-spline functions in (8) is based on a particular choice of triangles and points in the domain (see Fig. 3). 1. For each vertex V i ∈ V , select a triangle Qiv =  Q iv,1 , Q iv,2 , Q iv,3 . 2. For each triangle Tm =  V i , R i , j , Z n  ∈ PS , select a point

2 t Qm = S ti,n + σmt ( V j − V i ), 3

σmt > 0,

where S ti,n = 23 V i + 13 Z n . 3. For each boundary edge ε =  V i , R i , j  ∈ ∂ EPS , select a point

2 Q e = S ei, j + σe ( V j − V i ), 3 where S ei, j =

2 V 3 i

σe > 0,

+ 13 R i , j .

The selection is constrained by the following properties that ensure the non-negativity of the B-splines as we verify later in this section.

6

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

Fig. 3. A choice of triangles and points satisfying Properties 1–3.

Property 1. For each vertex V i ∈ V , the triangle Qiv contains V i and all the domain points of D 3 (PS ) at distance 1 from V i . Property 2. For each vertex V i ∈ V and each macro-triangle Tn =  V i , V j , V k  ∈  with a vertex at V i , the triangle Qti,n = t t ,Qm  S ti,n , Q m   contains the domain points

1 3

Vi +

1 3

R i, j +

1 3

Z n,

1 3

Vi +

2 3

Z n,

1 3

Vi +

1 3

R i ,k +

1 3

Z n,

t t where Q m and Q m  are the points associated with the triangles Tm =  V i , R i , j , Z n  and Tm =  V i , R i ,k , Z n  of PS inside Tn and sharing the vertex V i .

Property 3. For each boundary edge ε =  V i , R i , j  ∈ ∂ EPS , the line segment Qe =  S ei, j , Q e  contains the domain point 2 R . 3 i, j

1 V 3 i

+

In the next two subsections, we detail the construction of the B-splines. The main idea is the determination of their Bernstein–Bézier forms based on lifting the constrained triangles and line segments from the domain. 3.2. B-splines associated with vertices The definition of the B-spline functions B iv,r , r = 1, 2, 3, associated with the vertex V i of  is based on a choice of suitable parameters for the interpolation of value and first order derivatives at the vertex V i . In terms of the Bernstein–Bézier form and the minimal determining set provided in Proposition 1, this is equivalent to specifying the Bézier ordinates corresponding to V i and two domain points at distance 1 from V i that are contained in M i . The Bézier ordinates corresponding to the points in all M j , j = i, are set to 0. As a consequence, the supports of B iv,r are contained in the triangles of  sharing the vertex V i . Suppose that

B iv,r ( V i ) = αi ,r ,

y

D x B iv,r ( V i ) = αix,r ,

D y B iv,r ( V i ) = αi ,r ,

(9)

such that

αi,1 + αi,2 + αi,3 = 1, αix,1 + αix,2 + αix,3 = 0, αiy,1 + αiy,2 + αiy,3 = 0.

(10)

As will be clear later, the B-splines r = 1, 2, 3, are the only basis functions in (8) with non-zero values at V i , so the conditions (10) are necessary for the B-splines to form a partition of unity. By choosing (αi ,1 , αi ,2 , αi ,3 ) to be the barycentric y y y coordinates of the vertex V i , and (αix,1 , αix,2 , αix,3 ) and (αi ,1 , αi ,2 , αi ,3 ) to be the barycentric coordinates of the vectors (1, 0) v v v v and (0, 1) with respect to the triangle Qi =  Q i ,1 , Q i ,2 , Q i ,3 , these conditions are naturally satisfied. More specifically, the parameters are uniquely determined by the matrix identity B iv,r ,

⎡

⎤⎡

⎤

⎤ ⎡ X iv,1 Y iv,1 1 αi,1 αi,2 αi,3 x yi 1 ⎢ αx αx αx ⎥ ⎢ X v Y v 1 ⎥ ⎣ i ⎣ i ,1 ⎦ = 1 0 0 ⎦, i ,2 i ,3 ⎦ ⎣ i ,2 i ,2 y y y v v 0 1 0 αi,1 αi,2 αi,3 X i ,3 Y i ,3 1 where (xi , y i ) are the Cartesian coordinates of V i , and ( X iv,r , Y iv,r ) are the Cartesian coordinates of Q iv,r .

(11)

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

7

Fig. 4. Schematic representation of the Bézier ordinates of a B-spline associated with a vertex.

Consider a B-spline B iv,r , r ∈ {1, 2, 3}, on a macro-triangle Tn =  V i , V j , V k  ∈ . Its Bernstein–Bézier form is depicted in Fig. 4. The choice of the interpolation parameters in (9) implies that

b1v,r = αi ,r , b2v,r = αi ,r + b3v,r = αi ,r + b4v,r = αi ,r +

(12a) 1 3 1 3 1 3

y

(12b)

y

(12c)

y

(12d)

(αix,r , αi ,r ) · ( R i , j − V i ), (αix,r , αi ,r ) · ( Z n − V i ), (αix,r , αi ,r ) · ( R i ,k − V i ),

where · denotes the dot product between vectors. Geometrically, these values can be viewed as the heights of the triangle that is obtained by lifting the triangle Qiv so that the height of the vertex above Q iv,r is 1, and the heights of the other two vertices are 0. More precisely, the ordinates b1v,r , b2v,r , b3v,r , b4v,r are the heights above the domain points V i , 2 V 3 i

1 3

2 V 3 i

2 V 3 i

+ 13 R i , j ,

Qiv

1 3

+ Zn, + R i ,k , respectively. Since by Property 1 the triangle contains all these domain points, the ordinates are non-negative. The construction now proceeds with the determination of the ordinates b6v,r , b7v,r , and b8v,r . Denote by Tm =  V i , R i , j , Z n  and Tm =  V i , R i ,k , Z n  the subtriangles of Tn in PS sharing the vertex V i . Suppose that the triangle Qti,n is lifted so that the height of the vertex above S ti,n is b3v,r , and the heights of the other two vertices are 0. Let b6v,r , b7v,r , and b8v,r be the heights of the lifted triangle above the points

1 3

Vi +

1 3

R i, j +

1 3 1 3

Vi +

1 3

Vi +

R i ,k +

1 3 2 3 1 3

Zn = 1 −

 Zn = 1 −

 Zn = 1 −

λ j ,i 2σmt

ξn, j t m

2σ

λk,i 2σmt 

it follows that

b6v,r

= 1−

λ j ,i 2σmt



− 

ξn,k

b3v,r ,

b7v,r

2σmt

2σ

S ti,n +

ξn, j 2σmt

S ti,n +

λk,i 2σmt 

−

1 V 3 i

+ 23 Z n , and

1 V 3 i

+ 13 R i ,k + 13 Z n , respectively. Since

t Qm ,



t m

= 1−

+ 13 R i , j + 13 Z n ,

λ j ,i

S ti,n +





1 V 3 i

(13a)

ξn, j t m

2σ

t Qm +

ξn,k 2σmt 

t Qm ,

(13b)

t Qm ,

ξn,k 2σmt 

(13c)



 b3v,r ,

b8v,r

= 1−

λk,i 2σmt 

 b3v,r .

(14a)

In this way, the C 1 condition (7) across the edge  V i , Z n  on the ordinates b3v,r , b6v,r , b7v,r , b8v,r is satisfied. Also, the ordinates are non-negative since by Property 2 the triangle Qti,n contains the four corresponding domain points.

The Bézier ordinate b5v,r is determined by the C 1 continuity conditions across the edge ε =  V i , R i , j . Let Tn ∈  be the triangle sharing the edge  V i , V j  with Tn . Then,

b5v,r = μn,n b6v,r + μn ,n b8 v,r ≥ 0,

(14b)

8

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

Fig. 5. The three B-splines associated with the interior vertex of .

where b8 v,r is the Bézier ordinate of B iv,r corresponding to the domain point of PS , then

b5v,r

= 1−

1 3

Vi +

2 3

+ 13 R i , j + 13 Z n . If ε is a boundary edge



λ j ,i 2σe

Since

1 V 3 i

b2v,r .

R i, j = 1 −

(14c)

λ j ,i



2σe

S ei, j +

λ j ,i 2σe

Q e ,

(15)

this corresponds to the height of the line segment Qe =  S ei, j , Q e  lifted so that the height of the vertex above S ei, j is b2v,r , and the height of the vertex above Q e is 0. Thus, b5v,r ≥ 0 by Property 3. The value of b9v,r is specified similarly. The Bézier ordinates corresponding to the domain points at distance 3 from V i are determined by the C 1 continuity conditions across the edges  R i , j , Z n  and  R i ,k , Z n , v v b10 ,r = λi , j b 5,r ,

v v b11 ,r = λi , j b 6,r ,

v v b12 ,r = λi , j b 7,r ,

(16a)

v b16 ,r

v b15 ,r

v b14 ,r

(16b)

=

λi ,k b9v,r ,

=

λi ,k b8v,r ,

=

λi ,k b7v,r ,

and at the triangle split point Z n , v v b13 ,r = ξn,i b 7,r .

(16c)

They are clearly all non-negative. The remaining ordinates are zero (see Fig. 4). In summary, all the Bézier ordinates of B iv,r are non-negative, which ensures that the B-spline is non-negative. Fig. 5 depicts an example of the three B-splines associated with a vertex. 3.3. B-splines associated with triangles and boundary edges t The B-spline B m associated with the triangle Tm =  V i , R i , j , Z n  of PS is defined by specifying the Bézier ordinate corresponding to the domain point 13 V i + 13 R i , j + 13 Z n . The rest of the Bézier ordinates corresponding to the domain points t of the minimal determining set provided in Proposition 1 are set to 0, ensuring that B m has a local support. Let Tn =  V i , V j , V k  ∈  be the macro-triangle that contains Tm , and let Tm =  V i , R i ,k , Z n  be the second triangle t of PS inside Tn with a vertex at V i . Consider the Bernstein–Bézier form of B m on Tn depicted in Fig. 6 (left). The Bézier t t t ordinate bt3 is determined by lifting the triangle Qti,n =  S ti,n , Q m ,Qm  from the domain so that the vertex above Q m is  t t at height 1, but the other two vertices corresponding to S ti,n and Q m  stay in the domain at height 0. The value of b 3 is set 1 1 1 to the height of the lifted triangle above the point 3 V i + 3 R i , j + 3 Z n , which by (13a) implies that

bt3 =

λ j ,i 2σmt

(17a)

.

In this way, all the Bézier ordinates corresponding to the points of the minimal determining set are specified, and the remaining ordinates can be calculated in view of the continuity conditions. By (7), the value of bt5 can be interpreted as the height of the lifted triangle above the domain point

bt5 =

ξn, j 2σmt

.

1 V 3 i

+ 23 Z n , which by (13b) implies that

(17b)

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

9

Fig. 6. Schematic representation of the Bézier ordinates of B-splines associated with a triangle (left) and a boundary edge (right) of PS .

Since by Property 2 the triangle Qti,n contains the domain points 13 V i + 13 R i , j + 13 Z n and 13 V i + 23 Z n , the ordinates bt3 and bt5 are non-negative. Furthermore, the C 1 continuity conditions across the edges  R i , j , Z n  and  R i ,k , Z n  and at the split point Z n imply

bt4 = λi , j bt3 ,

bt6 = λi , j bt5 ,

bt7 = λi ,k bt5 ,

bt8 = ξn,i bt5 ,

(17c)

and all these ordinates are clearly non-negative. The values of bt1 and bt2 are determined according to the type of the edge  V i , V j . If  V i , V j  is a boundary edge, the ordinates are both equal to 0; otherwise

bt1 = μn,n bt3 ,

bt2 = μn,n bt4 ,

(17d)

where n is the index of the second triangle

Tn ∈  with the edge  V i , V j . All the remaining ordinates are zero. We t conclude that B m is non-negative since all the Bézier ordinates are non-negative. The B-spline B e associated with the boundary edge ε =  V i , R i , j  of PS is defined similarly to a B-spline associated with a triangle and may in fact be seen as its degenerate case. Since there is only one triangle of PS with the edge ε , B e takes the role of the B-spline that would be associated with the non-existing second triangle of PS sharing the edge ε . The Bernstein–Bézier form of B e on a macro-triangle that contains ε is depicted in Fig. 6 (right). The Bézier ordinate be1 that is part of the minimal determining set provided in Proposition 1 is specified by lifting the line segment Qe =  S ei, j , Q e . Suppose that the endpoint corresponding to Q e is at height 1, and the endpoint corresponding to S ei, j stays in the domain at height 0. The value of be1 is determined as the height of the line segment above the point 13 V i + 23 R i , j . By (15), the Bézier ordinate is given by

be1 =

λ j ,i 2σe

(18a)

,

and is non-negative by Property 3. All the other Bézier ordinates corresponding to the points in the minimal determining set are set to 0, which results in all Bézier ordinates being equal to 0, except

be2 = λi , j be1 .

(18b) e

Thus, the B-spline B  is non-negative. Fig. 7 shows an example of B-splines associated with triangles and a boundary edge of PS . 4. Properties of the B-spline representation In what follows, we list the key properties of the B-spline functions and the B-spline representation (8) provided in Section 3. We show that the B-splines form a locally supported stable basis and a convex partition of unity. Furthermore, we derive the Greville points of the B-spline representation and show, in general, how the B-spline coefficients of any s ∈ S13 (PS ) can be determined. Finally, we present a calculation of the Bernstein–Bézier form of s expressed in the B-spline representation which allows a stable evaluation of its values. 4.1. Linear independence of the B-splines To show that the B-splines form a basis of S13 (PS ), it suffices to prove that they are linearly independent. Let s ∈ S13 (PS ) be equal to 0, and consider its B-spline representation (8). Evaluating s, D x s, and D y s at the vertex V i ∈ V implies

10

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

Fig. 7. B-splines associated with a triangle (left, middle) and boundary edge (right) of PS .

c iv,1 αi ,1 + c iv,2 αi ,2 + c iv,3 αi ,3 = 0, c iv,1 αix,1 + c iv,2 αix,2 + c iv,3 αix,3 = 0, y

y

y

c iv,1 αi ,1 + c iv,2 αi ,2 + c iv,3 αi ,3 = 0. y

By (11), the vectors (αi ,r , αix,r , αi ,r ), r = 1, 2, 3, are linearly independent, which implies that c iv,r = 0. Furthermore, consider t the coefficient cm associated with the triangle Tm =  V i , R i , j , Z n  of PS . The Bézier ordinate of s corresponding to the 1 domain point 3 V i + 13 R i , j + 13 Z n is equal to 3 

c iv,r

1−

r =1

λ j ,i

 t b3v,r + cm

t m

2σ

λ j ,i 2σmt

= 0.

t Since c iv,r = 0, r = 1, 2, 3, the coefficient cm is 0 as well. By the same line of arguments, c e = 0 for every boundary edge ε of PS . Therefore, the coefficients of the B-spline representation (8) of s ≡ 0 are all equal to 0, and the B-splines must be linearly independent.

4.2. Convex partition of unity of the B-splines The proof of the property that the B-splines form a partition of unity, i.e., |V |  3 

B iv,r (x, y ) +

i =1 r =1

| PS | 

t Bm (x, y ) +

m =1

|∂ EPS |

B e (x, y ) = 1,

(19)

=1

is similar to the proof of the linear independence property. Indeed, it amounts to showing that the coefficients of the B-spline representation (8) of s ≡ 1 are all equal to 1. The evaluation of s and its first order partial derivatives at V i ∈ V yields the system of linear equations

c iv,1 αi ,1 + c iv,2 αi ,2 + c iv,3 αi ,3 = 1, c iv,1 αix,1 + c iv,2 αix,2 + c iv,3 αix,3 = 0, y

y

y

c iv,1 αi ,1 + c iv,2 αi ,2 + c iv,3 αi ,3 = 0, t with the unique solution (c iv,1 , c iv,2 , c iv,3 ) = (1, 1, 1). The coefficient cm associated with the triangle Tm =  V i , R i , j , Z n  of PS can now be evaluated from the Bézier ordinate of s corresponding to the domain point 13 V i + 13 R i , j + 13 Z n , which is equal to

3  r =1

3

c iv,r

1−

λ j ,i 2σmt

 t b3v,r + cm

λ j ,i 2σmt

= 1.

e v v t Since r =1 b 3,r = 1 and c i ,r = 1, r = 1, 2, 3, the coefficient cm is equal to 1 as well. Similarly, c  = 1 for every boundary edge ε of PS . The partition of unity together with the non-negativity of the B-splines (shown in Section 3) implies that the value of s ∈ S13 (PS ) at any point of the domain  can be described as a convex combination of the coefficients in the B-spline representation of s.

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

11

4.3. Greville points of the B-spline representation Next, let us consider the representation of the polynomials x and y of P1 in terms of (8). Since by (3) the Bézier ordinates of the polynomials x and y are equal to the Cartesian coordinates of the domain points, this amounts to finding suitable expressions for the domain points by means of the Bézier ordinates of the B-splines. The construction of the B-splines B iv,r , r = 1, 2, 3, on the triangle Tn =  V i , V j , V k  ∈  implies that

b1v,1 Q iv,1 + b1v,2 Q iv,2 + b1v,3 Q iv,3 = V i , 2 1 b2v,1 Q iv,1 + b2v,2 Q iv,2 + b2v,3 Q iv,3 = V i + R i , j , 3 3 2 1 v v v v v v b3,1 Q i ,1 + b3,2 Q i ,2 + b3,3 Q i ,3 = V i + Z n , 3 3 2 1 b4v,1 Q iv,1 + b4v,2 Q iv,2 + b4v,3 Q iv,3 = V i + R i ,k . 3 3

(20a) (20b) (20c) (20d)

t Moreover, taking into account the B-spline B m associated with the triangle Tm =  V i , R i , j , Z n  ∈ PS , (13a) implies that

3 

t b6v,r S ti,n + bt3 Q m =

r =1

1 3

Vi +

1 3

R i, j +

1 3

(20e)

Z n.

Similarly, if ε =  V i , R i , j  is a boundary edge of PS , 3  r =1

1 2 b5v,r S ei, j + be1 Q e = V i + R i , j , 3 3

(20f)

as can be observed from (15). The equations in (20) show how to express all domain points of the minimal determining set provided in Proposition 1. They imply that |V |  3 

Q iv,r B iv,r (x, y ) +

i =1 r =1

| PS | 

t t Qm Bm (x, y ) +

m =1

|∂ EPS |

Q e B e (x, y ) = (x, y ).

(21)

=1

t Therefore, the points Q iv,r , Q m , Q e are the Greville points of the B-spline representation (8).

4.4. Coefficients of the B-spline representation So far, from (19) and (21) we know how any p ∈ P1 can be expressed in the B-spline representation (8). We now show how to calculate the coefficients of (8) for any s ∈ S13 (PS ) by exploiting the blossoming principle. The vertex B-spline interpolation property (9) implies that

c iv,1 αi ,1 + c iv,2 αi ,2 + c iv,3 αi ,3 = s( V i ), c iv,1 αix,1 + c iv,2 αix,2 + c iv,3 αix,3 = D x s( V i ), y

y

y

c iv,1 αi ,1 + c iv,2 αi ,2 + c iv,3 αi ,3 = D y s( V i ), for every V i ∈ V . Solving this system of equations with the help of (11) reveals that

c iv,r = s( V i ) + ∇ s( V i ) · ( Q iv,r − V i ),

r = 1, 2, 3,

(22)

where ∇ s denotes the gradient of s. Let sm ∈ P3 be the restriction of s to any triangle Tm ∈ PS with a vertex at V i . According to (6) and (22), the coefficient c iv,r can be expressed in terms of the blossom as

c iv,r = B sm ( V i , V i , V i ) + 3B sm ( V i , V i , Q iv,r − V i ) = B sm ( V i , V i , Q iv,r + 2( Q iv,r − V i )).

(23)

t Let us further consider the coefficient cm associated with a triangle Tm =  V i , R i , j , Z n  ∈ PS contained in the macrotriangle Tn =  V i , V j , V k  ∈ . Combining (5) and the Bernstein–Bézier form of the B-splines elaborated in Section 3 provides us with the relation

B sm ( V i , R i , j , Z n ) =

3 

t t c iv,r b6v,r + cm b3 .

r =1

A short calculation using (14a), (20c), (23), and properties of the blossom shows that

(24)

12

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

3  r =1

c iv,r b6v,r

= 1−

= 1−

λ j ,i

 3

2σmt

λ j ,i



2σmt

c iv,r b3v,r

r =1



B sm V i , V i ,

3 

 b3v,r ( Q iv,r

+ 2( Q

v i ,r

− V i ))

r =1



λ j ,i B sm ( V i , V i , Z n ). = 1− t 2σm

Inserting this into (24) and using (17a) gives t t t t cm = B sm ( V i , V i + 2σmt ( V j − V i ), Z n ) = B sm ( V i , Q m + (Q m − V i) + ( Q m − Z n ), Z n ).

(25)

If ε =  V i , R i , j  is a boundary edge of PS , the coefficient c  can be obtained similarly from e

B sm ( V i , R i , j , R i , j ) =

3 

c iv,r b5v,r + c e be1 .

r =1

With the use of (14c), (20b), (23), and (18a) the coefficient can be expressed as

c e = B sm ( V i , V i + 2σe ( V j − V i ), R i , j ) = B sm ( V i , Q e + ( Q e − V i ) + ( Q e − R i , j ), R i , j ).

(26)

Summarizing, thanks to (23), (25), and (26), the coefficients of the B-spline representation (8) of any s ∈ S13 (PS ) can be compactly written in terms of blossoming as

c iv,r = Sm ( V i , V i , Q iv,r ),

t t cm = Sm ( V i , Z n , Q m ),

c e = Sm ( V i , R i , j , Q e ),

(27)

where

Sm ( P 1 , P 2 , P 3 ) = B sm ( P 1 , P 2 , P 3 + ( P 3 − P 1 ) + ( P 3 − P 2 )), and sm is the restriction of s to the triangle Tm =  V i , R i , j , Z n  ∈ PS , and edge ε =  V i , R i , j  ∈ ∂ EPS . An explicit Hermite interpolation expression of bivariate cubic polynomials in Bernstein–Bézier form via blossoming can be found in Speleers (2011, Example 7.5). When simplifying this expression in our context, we get

Sm ( P 1 , P 2 , P 3 ) = sm ( P 1 ) +

1 3 1

∇ sm ( P 1 ) · (( P 2 − P 1 ) + 2( P 3 − P 1 ) + ( P 3 − P 2 ))

+ ( P 2 − P 1 ) · H sm ( P 1 ) · (2( P 3 − P 1 ) + ( P 3 − P 2 )) 6

1

= sm ( P 1 ) + ∇ sm ( P 1 ) · ( P 3 − P 1 ) + ( P 2 − P 1 ) · H sm ( P 1 ) · (3 P 3 − 2 P 1 − P 2 ), 6

where H sm denotes the Hessian matrix of sm , and



(u 1 , v 1 ) · H sm · (u 2 , v 2 ) = u 1



v 1 H sm



u2 v2

 .

4.5. Stable computation of the Bernstein–Bézier form To evaluate a spline s ∈ S13 (PS ) at a particular point of the domain, one can calculate the Bernstein–Bézier form of s on a triangle of PS that contains the evaluation point and then use the de Casteljau algorithm to obtain the value. We now show how to determine the Bézier ordinates of s expressed in the form (8) as convex combinations of the B-spline coefficients. Consider a triangle Tn =  V i , V j , V k  ∈ , and let the Bézier ordinates of s be denoted as in Fig. 8. They can be derived from the Bézier ordinates of the B-splines specified in (12), (14), (16), (17), and (18). It follows from (12) that

b1 = αi ,1 c iv,1 + αi ,2 c iv,2 + αi ,3 c iv,3 , b2 = τ2,1 c iv,1 + τ2,2 c iv,2 + τ2,3 c iv,3 , b3 = τ3,1 c iv,1 + τ3,2 c iv,2 + τ3,3 c iv,3 , b4 = τ4,1 c iv,1 + τ4,2 c iv,2 + τ4,3 c iv,3 , where (αi ,1 , αi ,2 , αi ,3 ), (τ2,1 , τ2,2 , τ2,3 ), (τ3,1 , τ3,2 , τ3,3 ), and (τ4,1 , τ4,2 , τ4,3 ) are the barycentric coordinates of the points t V i , 23 V i + 13 R i , j , 23 V i + 13 Z n , and 23 V i + 13 R i ,k with respect to the triangle Qiv =  Q iv,1 , Q iv,2 , Q iv,3 . Furthermore, let cm

t and cm  be the coefficients associated with the triangles Tm =  V i , R i , j , Z n  and Tm =  V i , R i ,k , Z n . Combining (17a)–(17b)

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

13

Fig. 8. Schematic representation of the Bézier ordinates of a cubic Powell–Sabin spline.

and (14a) reveals that t t b6 = τ6,1 b3 + τ6,2 cm + τ6,3 cm , t t b7 = τ7,1 b3 + τ7,2 cm + τ7,3 cm , t t b8 = τ8,1 b3 + τ8,2 cm + τ8,3 cm ,

where (τ6,1 , τ6,2 , τ6,3 ), (τ7,1 , τ7,2 , τ7,3 ), and (τ8,1 , τ8,2 , τ8,3 ) are the barycentric coordinates of the points 1 V + 23 Z n , and 13 V i + 13 3 i t cm since 6,3 = 0, and b8

R i ,k +

1 3

Z n with respect to

1 V 3 i

+ 13 R i , j + 13 Z n ,

Qti,n

=

1 V 3 i

+ 23 R i , j with respect to Qe =  S ei, j , Q e . Otherwise, the C 1 condi-

 S ti,n ,

Q

t m,

Q

t . m

As a matter of fact, b6 only depends on b3 and

t only depends on b3 and cm  since τ8,2 = 0. The calculation of b 5 depends on the type of the edge ε =  V i , R i, j . If ε is a boundary edge, then (18a) and (12b) show that

τ

b5 = τ5,1 b2 + τ5,2 c e , where (τ5,1 , τ5,2 ) are the barycentric coordinates of tions across the edge  V i , R i , j  imply that

b5 = μn,n b6 + μn ,n b8 , where n is the index of the triangle Tn ∈ PS sharing the edge  V i , R i , j  with Tn , and b8 is the Bézier ordinate of s corresponding to the domain point

1 V + 13 3 i

R i , j + 13 Z n . The value of b9 is determined analogously. The values of b10 , . . . , b18

and b19 , . . . , b27 can be calculated in a similar way. The remaining Bézier ordinates can be deduced easily from the C 1 conditions across the edges  R i , j , Z n ,  R j ,k , Z n ,  R i ,k , Z n ,

b28 = λi , j b5 + λ j ,i b18 ,

b29 = λi , j b6 + λ j ,i b17 ,

b30 = λi , j b7 + λ j ,i b16 ,

b31 = λ j ,k b14 + λk, j b27 ,

b32 = λ j ,k b15 + λk, j b26 ,

b33 = λ j ,k b16 + λk, j b25 ,

b34 = λk,i b23 + λi ,k b9 ,

b35 = λk,i b24 + λi ,k b8 ,

b36 = λk,i b25 + λi ,k b7 ,

and at the triangle split point Z n ,

b37 = ξn,i b7 + ξn, j b16 + ξn,k b25 . Clearly, by Properties 1–3, only convex combinations are used in the above calculation of all Bézier ordinates of the cubic spline s in the form (8) starting from its B-spline coefficients. 4.6. Stability of the B-spline basis The stability of a basis ensures that a relatively small perturbation of coefficients does not have a large impact on a relative error of the spline value. The B-spline representation (8) is stable if there exist constants κ and η such that

κ c ∞ ≤ s∞, ≤ η c ∞ S13 (PS ),

(28)

for any s ∈ where c is a vector of all coefficients of s that appear in (8). The condition number η/κ ≥ 1 measures sensitivity of the basis and should be as close to 1 as possible. By (19) and the non-negativity of the B-splines, we can

14

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

set η = 1. The constant κ can be derived based on the Bernstein–Bézier form of s and standard approximation results for polynomials in Bernstein–Bézier form. By (22), the coefficient c iv,r corresponding to the vertex V i ∈ V is bounded by

√  v       c  ≤ |s( V i )| + 2 max | D x s( V i )| ,  D y s( V i )  Q v − V i  . i ,r i ,r

Let TPS be a triangle of PS with a vertex at V i , and consider the Bernstein–Bézier form of s restricted to TPS . Denote the radius of the largest disk contained in TPS by ρPS , and let K be a constant satisfying (4) with d = 3. Then, by Lai and Schumaker (2007, Theorem 2.19) we have







max | D x s( V i )| ,  D y s( V i ) ≤

3K

ρPS

s∞,TPS .

This implies that





κiv c iv,r  ≤ s∞,TPS ≤ s∞, , where 1/κiv = 1 + K iv , and K iv is a constant depending on the ratio between the diameter of Qiv and ρPS . This indicates that the triangle Qiv should be small within the constraints prescribed by Property 1. If Qiv is a triangle with minimal area, then practically the same analysis as in Maes et al. (2004, Theorems 4.3 and 6.1) shows that K iv can be determined so that it only depends on the smallest angle of PS . t that provides a lower bound for the ratio Let Tm =  V i , R i , j , Z n be a triangle of PS . To determine the constant κm t  between s∞, and cm , we consider the Bernstein–Bézier form of s restricted to Tm . From Section 4.5 we know that t cm

=

1

τ6,2

b6 −

1 − τ6,2

τ6,2



b3 =

2σmt

λ j ,i

b6 + 1 −

2σmt



λ j ,i

b3 .

Let again K be a constant satisfying (4) with d = 3. Then,

t  m

 t  c  ≤ K 1 + 4σ m λ j ,i

 ⎞  t   Q m − S ti,n   ⎠ s∞,Tm . = K ⎝1 + 6   R i, j − V i  ⎛

s∞,Tm

t This shows that the line segment  S ti,n , Q m  should be as short as possible. For a fixed parameter

Property 2, the constant estimate

1

λ j ,i

σmt such that Qti,n satisfies

κ can be chosen so that it only depends on the smallest angle θPS of PS . This follows from the t m

  V j − V i  1 ≤ = ,  R i , j − V i  sin(θPS )2

which can be proven between s∞, and κiv , κmt and κe .

e bye using the law of sine twice. For a boundary edge ε of PS , the constant κ bounding the ratio c  can be derived in a similar way. The constant κ in (28) can then be chosen as the minimum of 

5. Control points and control structures Control points are a useful tool for the geometric modelling of spline surfaces. By combining (8) and (21), it is natural to define the control points of the B-spline representation as

c iv,r = ( Q iv,r , c iv,r ),

t t t cm = (Q m , cm ),

c e = ( Q e , c e ),

for all possible indices (i , r ), m, and . As a consequence of (19) and the non-negativity of the B-splines, the graph of s ∈ S13 (PS ) lies in the convex hull of its control points. The control points corresponding to a vertex V i of  form a control triangle Civ = c iv,1 , c iv,2 , c iv,3  whose projection to the domain is the triangle Qiv . By (22), the triangle Civ lies in the plane tangent to s at the point V i , and so it can be used to intuitively handle the graph surface in the neighbourhood of V i . In order for Civ to mimic the surface well, the triangle Qiv should be as small as possible within the constraints of Property 1. This leads to the problem of finding a triangle with minimal area that contains a certain set of points. Such a triangle could be found using the algorithm provided by Klee and Laskowski (1985) or O’Rourke et al. (1986). Alternatively, we could also formulate and solve the following quadratic y y optimization problem. By (11), the area of Qiv is small if αix,1 αi ,2 − αi ,1 αix,2 is large. Thus, we have the objective y

y

max (αix,1 αi ,2 − αi ,1 αix,2 )

(29)

subjected to the constraints (10) and the conditions that the values in (12) are greater than or equal to 0 for every triangle Tn ∈  with a vertex at V i , which ensure that Property 1 is satisfied. A detailed derivation of such an optimization problem can be found in Dierckx (1997); Speleers (2010b).

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

15

The control points corresponding to the triangle B-splines can be organized into control triangles as well. Let Tn be a triangle of  with a vertex at V i . Denote by c ti,n the point on the triangle Civ that lies above the domain point S ti,n . For t t triangles Tm and Tm of PS sharing the vertex V i and lying inside Tn ∈ , the triangle Cit,n = c ti,n , c m , cm   is a lifting

of Qti,n . It does not have a tangent property similar to Civ but it lies in the plane spanned by two triangles of the Bézier control net of s above  V i , Z n . Thus, Cit,n should be small as well, but according to Property 2 and (14a), the conditions

σmt ≥

λ j ,i 2

ξn, j

,

2σmt

+

ξn,k 2σmt 

λk,i

σmt  ≥

≤ 1,

(30)

2

must not be violated. One obvious choice is σmt = σmt  = 12 . It greatly simplifies the construction of the B-splines but the drawback is that for any edge  V i , V j  of  and any Tn =  V i , V j , V k  ∈  the Greville points associated with the triangles  V i , R i , j , Z n  and  V j , R i , j , Z n  coincide with 13 V i + 13 V j + 13 Z n . Therefore, the choice

σmt =

λ j ,i

2

for small values of

ξn, j ξn,k + λ j ,i λk,i



γ , σmt  =

λk,i

2

ξn, j ξn,k + λ j ,i λk,i



γ , γ ≥ 1,

γ is preferable. The conditions (30) are met if γ satisfies

ξn, j ξn,k 1 + ≥ . λ j ,i λk,i γ

(31)

In the typical case that Z n is contained in the triangle  R i , j , R j ,k , R i ,k , the left-hand side of (31) is greater than or equal t t to 1, and we can choose γ = 1. The point 13 V i + 23 Z n then lies on the line segment  Q m ,Qm  . Otherwise, clearly our best option is

γ=

ξn, j ξn,k + λ j ,i λk,i

−1 ,

and so, in summary,

σmt =

λ j ,i 2

 max

 ξn, j ξn,k + ,1 , λ j ,i λk,i

σmt  =

λk,i 2

 max

 ξn, j ξn,k + ,1 . λ j ,i λk,i

(32)

If additionally the left-hand side of (31) is less than 2, or, equivalently, the point Z n lies inside the triangle  V i , V i + t t 2( R i , j − V i ), V i + 2( R i ,k − V i ), this choice ensures that Q m is inside Tm and Q m  is inside Tm . The control point c e corresponding to the boundary edge ε =  V i , R i , j  of PS can be thought of as the endpoint of the control line Ce = c ei, j , c e , where c ei, j is the point on Civ above S ei, j . The line segment Ce is a lifting of Qe . To obtain the shortest line satisfying Property 3, we set

σe =

λ j ,i 2

(33)

,

but the choices σe = 12 or σe = σmt , where Tm is the triangle of PS with the edge ε , are also acceptable. The above selection of control triangles and control lines is directly related to the geometric interpretation of the B-spline construction (see Properties 1–3). This natural choice, however, might not be very practical from a design point of view because not all the vertices of the control triangles Cit,n and the control lines Ce are free to choose as they depend on the other control triangles Civ . Therefore, alternative control structures might be more appealing. For example, one could opt for one control triangle per vertex and one control quadrilateral per edge as follows. 1. For each vertex V i of , set the control triangle Civ = c iv,1 , c iv,2 , c iv,3  collecting the three control points associated with the vertex. t t t t , cm , cm , cm  collecting the 2. For each interior edge εl =  V i , V j  of , set the control quadrilateral Clte = c m 1 2 3 4 control points associated with the triangles Tm1 =  V i , R i , j , Z n , Tm2 =  V j , R i , j , Z n , Tm3 =  V i , R i , j , Z n , Tm4 =  V j , R i , j , Z n . The control quadrilateral is chosen to be the bilinear interpolant parametrized by



1−u

u

 t  cm 1

t cm 2

t cm 3

t cm 4



1−v v

 ,

(u , v ) ∈ [0, 1] × [0, 1].

t e e t 3. For each boundary edge εl =  V i , V j  of , set the control quadrilateral Clte = c m , cm  , c  , c   collecting the control points associated with the triangles Tm =  V i , R i , j , Z n , Tm =  V j , R i , j , Z n , and the boundary edges ε =  V i , R i , j , ε =  V j , R i, j .

16

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

Fig. 9. A spline surface that is a linear combination of six vertex B-splines, six triangle B-splines, and two edge B-splines. (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

As mentioned earlier, each control triangle Civ is tangent to the spline surface s at the point V i . There is no similar tangent property for the control quadrilaterals Clte , but collapsing these quadrilaterals to lines gives rise to some additional smoothness properties. This will be discussed in Section 6.2 in more detail. Fig. 9 (right) shows a spline surface in the B-spline representation (8) together with its control points arranged in red triangles Civ and red quadrilaterals Clte . The free parameters of the B-splines are set as in (29), (32), and (33) in order to satisfy Property 1, Property 2, and Property 3, respectively. The corresponding Powell–Sabin refined triangulation and the Greville points are shown in Fig. 9 (left). In addition, the projected triangles Civ and quadrilaterals Clte are coloured in red; the projected triangles Cit,n and lines Ce are coloured in blue. 6. Relations to other macro-element spaces The spline space S13 (PS ) contains some well-known macro-element subspaces, including the classical quadratic Powell– Sabin splines and cubic Clough–Tocher splines. In this section, we explore how splines in these subspaces can be expressed by means of the B-spline representation (8). Of course, we could use the general formulas in (27) to calculate the corresponding B-spline coefficients, but we will provide specific (easier) formulas fine-tuned for the different subspaces. As already mentioned in Section 5, the construction of the B-splines in Section 3 can be simplified by choosing σmt = 12 for every triangle Tm of PS and σe = 12 for every boundary edge ε of PS . In the following, we make this choice to t simplify the calculations. This has no effect on the generality of the results since the coefficients c˜ iv,r , c˜ m , c˜ e with respect to t e the B-splines  B iv,r ,  Bm , B  obtained with the same triangles Qiv but different parameters σ˜ mt , σ˜ e can be expressed easily in e v t terms of c i ,r , cm , c  . By comparing the Bernstein–Bézier forms of the B-splines, one can verify that

c˜ iv,r = c iv,r ,

t t c˜m = 2σ˜ mt cm + 1 − 2σ˜ mt

3 !

t c iv,r bm ,r ,

c˜ e = 2σ˜ e c e + 1 − 2σ˜ e

r =1

where

t bm ,r

denotes the Bézier ordinate of

B iv,r

3 !

c iv,r be,r ,

r =1

corresponding to the domain point

 V i , R i , j , Z n  ∈ PS , and b,r denotes the Bézier ordinate of boundary edge ε =  V i , R i , j  ∈ ∂ EPS . e

B iv,r

2 V 3 i

+

1 3

Z n in the triangle Tm =

corresponding to the domain point

2 V 3 i

+ 13 R i , j on the

6.1. C 1 quadratic splines on Powell–Sabin triangulations As developed by Dierckx (1997), an element of s ∈ S12 (PS ) can be expressed in a B-spline representation of the form

s(x, y ) =

|V |  3 

PS2 c iPS2 ,r B i ,r (x, y ).

(34)

i =1 r =1

The construction of the quadratic B-splines in (34) is based on a set of triangles associated with the vertices of . Nonnegativity of the basis requires certain geometric constraints on these triangles to be satisfied. Such constrained triangles

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

17

Fig. 10. Schematic representation of the Bézier ordinates of a quadratic Powell–Sabin B-spline (left) and the Bernstein–Bézier form obtained by raising the degree of polynomial patches by one (right).

satisfy Property 1 and can be used as the triangles Qiv in the construction of the cubic B-splines B iv,r . In the following we show how s can be represented in the form (8), which can be regarded as the degree elevation of s. We assume that the 1 e t t B-splines B iv,r , B m , B e are determined by the same triangles Qiv as B PS2 i ,r and the parameters σm = σ = 2 . We start by expressing the single quadratic B-spline B PS2 in terms of the cubic B-splines in (8). Comparing the B-spline i ,r supports, the B-spline B PS2 can be written as a linear combination of the cubic B-splines associated with the vertex V i and i ,r t the B-splines B m and B e associated with the triangles Tm of PS and the boundary edges ε of PS inside the triangles of  sharing the vertex V i . on a triangle Tn =  V i , V j , V k  ∈  depicted in Fig. 10 (left). According to Consider the Bernstein–Bézier form of B PS2 i ,r Dierckx (1997),

bPS2 1,r = αi ,r , bPS2 2,r = αi ,r + bPS2 3,r = αi ,r + bPS2 4,r = αi ,r +

(35a) 1 2 1 2 1 2

y

(35b)

y

(35c)

y

(35d)

(αix,r , αi ,r ) · ( R i , j − V i ), (αix,r , αi ,r ) · ( Z n − V i ), (αix,r , αi ,r ) · ( R i ,k − V i ),

and PS2 bPS2 5,r = λi , j b 2,r ,

PS2 bPS2 6,r = λi , j b 3,r ,

PS2 bPS2 7,r = ξn,i b 3,r ,

PS2 bPS2 8,r = λi ,k b 3,r ,

PS2 bPS2 9,r = λi ,k b 4,r .

(35e)

The Bernstein–Bézier form obtained by raising the degree to three is shown in Fig. 10 (right). After applying standard formulas for degree elevation of polynomials in Bernstein–Bézier form, we get PS2 bˆ PS2 1,r = b 1,r ,

bˆ PS2 2,r =

1 3

bPS2 1,r +

2 3

bPS2 2,r ,

bˆ PS2 3,r =

1 3

bPS2 1,r +

2 3

bPS2 3,r ,

bˆ PS2 4,r =

1 3

bPS2 1,r +

2 3

bPS2 4,r .

(36)

A short calculation using (12) shows that v bˆ PS2 1,r = b 1,r ,

v bˆ PS2 2,r = b 2,r ,

v bˆ PS2 3,r = b 3,r ,

v bˆ PS2 4,r = b 4,r .

(37)

This implies that the coefficient corresponding to B iv,r in the cubic B-spline representation of B PS2 i ,r is equal to 1, and the coefficients of the other two vertex B-splines associated with V i are equal to 0. We now proceed with the determination of the coefficients in the linear combination representing B PS2 that cori ,r t t the coefficient that corresponds to B m associated with the triangle respond to the triangle B-splines. Denote by am

t Tm =  V i , R i , j , Z n . With exception of the vertex B-splines associated with V i , B m is the only B-spline with non-zero Bézier ordinate corresponding to the domain point 13 V i + 13 R i , j + 13 Z n . Thus,

1 PS2 1 PS2 1 PS2 t λ j ,i am + b6v,r = bˆ PS2 6,r = b 2,r + b 3,r + b 6,r , 3

and so by (14a), (36)–(37), and (35e),

3

3

18

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

t am



1

=

1 3

λ j ,i

bPS2 2,r

+

1 3

bPS2 3,r

+

1 3

 bPS2 6,r

− λi , j

By (35a)–(35c), this can be expressed as

1

t am =

3



1

αi,r + (αix,r , αiy,r ) · ( Z n − V i ) + 2

1 3

1 3

bPS2 1,r

+

2 3

 bPS2 3,r

1

1

3

3λ j , i

= bPS2 3,r +

"

#

PS2 bPS2 2,r − λi , j b 1,r .



1

αi,r + (αix,r , αiy,r ) · ( V j − V i ) .

(38)

2

t t On the other hand, the coefficient am  corresponding to the B-spline B m associated with the triangle Tm =  V j , R i , j , Z n  satisfies

1

1

3

3

t PS2 ˆ PS2 λi , j am bPS2  = b 18,r = 6,r = λi , j b 3,r ,

which by (35c) implies t am  =

1

3



1

αi,r + (αix,r , αiy,r ) · ( Z n − V i ) .

(39)

2

The coefficients corresponding to the B-splines that are associated with the triangles  V i , R i ,k , Z n  and  V k , R i ,k , Z n  can be computed by symmetry. Finally, the coefficients corresponding to the B-splines that are associated with the triangles  V j , R j ,k , Z n  and  V k , R j ,k , Z n  are equal to 0. The derivation of the coefficients in the linear combination representing B PS2 i ,r that correspond to the boundary edge B-splines is very similar. Suppose, for example, that  V i , V j  is a boundary edge of . Then, the coefficient ae corresponding to B e that is associated with ε =  V i , R i , j  can be expressed as

ae =

1 3



1

1

αi,r + (αix,r , αiy,r ) · ( R i, j − V i ) + 2

3



1

αi,r + (αix,r , αiy,r ) · ( V j − V i ) ,

(40)

2

and the coefficient ae corresponding to B e that is associated with ε =  V j , R i , j  as

ae =

1

3



1

αi,r + (αix,r , αiy,r ) · ( R i, j − V i ) .

(41)

2

Symmetric expressions hold for the coefficients corresponding to the B-splines associated with the edges  V i , R i ,k  and  V k , R i ,k  if they are boundary edges of PS . On the contrary, the coefficients corresponding to the B-splines associated with  V j , R j ,k  and  V k , R j ,k  are equal to 0. in the cubic B-spline form (8) is a useful ingredient to prove the following The expression of the quadratic B-spline B PS2 i ,r theorem. Theorem 1. Suppose that s ∈ S13 (PS ) is represented in the form (8) with the B-splines determined by choosing 1, . . . , |PS |, and σ = e

t cm =

1

 = 1, . . . , |∂ EPS |. For each triangle Tm =  V i , R i , j , Z n  of PS , suppose 1

1 3

1 3

1

1

3

3

(ρi ,1 c iv,1 + ρi ,2 c iv,2 + ρi ,3 c iv,3 ) + (υi ,1 c iv,1 + υi ,2 c iv,2 + υi ,3 c iv,3 ) + (ρ j ,1 c vj,1 + ρ j ,2 c vj,2 + ρ j ,3 c vj,3 ),

where (υh,1 , υh,2 , υh,3 ), (ρh,1 , ρh,2 , ρh,3 ), and (ζh,1 , ζh,2 , ζh,3 ) are the barycentric coordinates of the points 1 2

R i , j , and

σmt = 12 , m =

(ζi ,1 c iv,1 + ζi ,2 c iv,2 + ζi ,3 c iv,3 ) + (υi ,1 c iv,1 + υi ,2 c iv,2 + υi ,3 c iv,3 ) + (ζ j ,1 c vj,1 + ζ j ,2 c vj,2 + ζ j ,3 c vj,3 ),

3 and for each boundary edge ε =  V i , R i , j  of PS , suppose c e =

3

1 , 2

1 V 2 h

1 2

+ Z n with respect to the triangle

Qhv

= Q

v , h ,1

Q

v , h ,2

Q

v  h ,3

for h = i , j. Then s

∈ S12 (PS ).

1 V 2 i

+ 12 V j ,

1 V 2 h

+

Proof. Using the above expression of the quadratic B-spline B PS2 i ,r in terms of the cubic B-splines (in particular the formulas (38)–(41)), we deduce that the B-spline form (8) of s ∈ S12 (PS ) has the following coefficients: for each vertex V i ∈ V ,

c iv,r = c iPS2 ,r ,

r = 1, 2, 3;

for each triangle Tm =  V i , R i , j , Z n  of PS , t cm =

1

s( V i ) +

3

1 2

 ∇ s( V i ) · ( Z n − V i )

 ∇ s( V i ) · ( V j − V i ) 3 2

 1 1 s( V j ) + ∇ s( V j ) · ( Z n − V j ) ; +

+

1

3

s( V i ) +

1

2

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

19

and for each boundary edge ε =  V i , R i , j  of PS ,

c e =

 ∇ s( V i ) · ( R i , j − V i ) 3 2

 1 1 s( V i ) + ∇ s( V i ) · ( V j − V i ) + 3 2

 1 1 s( V j ) + ∇ s( V j ) · ( R i , j − V j ) . + 1

s( V i ) +

3

1

2

The expression in (22) for the coefficients c iv,r completes the proof.

2

Theorem 1 immediately leads to a degree elevation formula for C 1 quadratic Powell–Sabin splines. Corollary 1. Let s ∈ S12 (PS ) be represented in the form (34). For each vertex V i ∈ V , we set

c iv,r = c iPS2 ,r ,

r = 1, 2, 3,

t and the coefficients cm and c e associated with any triangle Tm of PS and any boundary edge ε of PS , respectively, are determined t as in Theorem 1. Then, s satisfies (8) with the B-splines B iv,r , B m , B e determined by the same triangles Qiv as B PS2 i ,r and the parameters

σmt = σe = 12 . Theorem 1 and Corollary 1 reveal that degree elevation can be nicely described in terms of control points. The spline s ∈ S12 (PS ) in the B-spline representation (34) is determined by the control triangles associated with the vertices V i ∈ V . t To represent s in the B-spline representation (8), we consider these triangles to be Civ . The control point c m associated with the triangle Tm =  V i , R i , j , Z n  of PS is then obtained as the barycentre of a triangle with vertices that lie in the planes spanned by Civ and C vj . More precisely, two vertices are in the plane of Civ above the points 12 V i + 12 Z n and 12 V i + 12 V j , and the third vertex is in the plane of C vj above the point 12 V j + 12 Z n . Similarly, the control point c e associated with the boundary edge ε =  V i , R i , j  of PS is the barycentre of the triangle with two vertices in the plane of Civ above the points 1 V 2 i

+ 12 R i , j and

1 V 2 i

+ 12 V j and one vertex in the plane of C vj above the point

1 Vj 2

+ 12 R i , j .

6.2. C 1 cubic super-splines on Powell–Sabin triangulations The B-splines in (8) obtained by choosing σmt = 12 for every triangle Tm of PS and σe = 12 for every boundary edge ε of PS have a few interesting smoothness properties. It can be rather easily verified that the B-spline B iv,r , r ∈ {1, 2, 3}, associated with V i ∈ V is globally C 1 and locally even C 2 at every triangle split point of PS and across every edge of PS connecting a triangle split point and an edge split point. Similarly, for any two triangles Tm =  V i , R i , j , Z n  and t t Tm =  V j , R i , j , Z n  of PS inside the same macro-triangle Tn ∈  sharing the edge  R i , j , Z n , the spline B m + Bm  is locally 2 C at every triangle split point of PS and across every edge of PS connecting a triangle split point and an edge split point. If  V i , V j  is a boundary edge, this also holds true for the sum B e + B e of the B-splines B e and B e associated with the t t boundary edges ε =  V i , R i , j  and ε =  V j , R i , j , respectively. A careful inspection reveals that the splines B iv,r , B m + Bm , e e and B  + B  coincide with the B-splines provided by Speleers (2015) that form a basis for the super-spline space

% & $ S13 (PS ) = s ∈ S13 (PS ) : s ∈ C 2 ( Z ), Z ∈ ZPS ; s ∈ C 2 (ε ), ε ∈ E$PS ,

$PS is the set of all edges of PS connecting a triangle split point where ZPS is the set of all triangle split points of PS and E and an edge split point. This observation proves the following result.

σmt = 12 , m = t t  = 1, . . . , |∂ EPS |. For each edge ε =  V i , V j  ∈ E and each triangle Tn ∈  with the edge ε , let cm = cm ,

Theorem 2. Suppose that s ∈ S13 (PS ) is represented in the form (8) with the B-splines determined by choosing 1 , 2

1, . . . , |PS |, and σ = t t where cm and cm  are the coefficients associated with the triangles Tm =  V i , R i , j , Z n  and Tm =  V j , R i , j , Z n  of PS , respectively. Additionally, if ε is a boundary edge, let c e = c e , where c e and c e are the coefficients associated with the boundary edges ε =  V i , R i, j  and ε =  V j , R i, j , respectively. Then s ∈ $ S13 (PS ). e

Theorem 2 specifies 2|E | additional conditions on a spline s ∈ S13 (PS ), reducing the number of its degrees of freedom to 3|V | + 2|E |. This number is exactly the dimension of the space $ S13 (PS ), see Speleers (2015). Geometrically, it leads to one control triangle per vertex of  (connecting three control points) and one control line per edge of  (connecting two control points transverse to the edge).

20

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

Another cubic subspace of S13 (PS ) of dimension 6|V | with local C 2 super-smoothness at the vertices has been considered by Chen and Liu (2008); Lamnii et al. (2014):

% &  S13 (PS ) = s ∈ S13 (PS ) : s ∈ C 2 ( V ), V ∈ V .

Thanks to the blossoming expressions of the B-spline coefficients in Section 4.4, we arrive at the following result. Theorem 3. Suppose that s ∈ S13 (PS ) is represented in the form (8) with the B-splines determined by choosing

σmt = 12 , m =

 = 1, . . . , |∂ EPS |. For each vertex V i ∈ V and any given polynomial qi ∈ P3 , let = Bqi ( V i , V i , 3 Q iv,r − t 2V i ). Moreover, for each triangle Tm =  V i , R i , j , Z n  of PS , let cm = Bqi ( V i , V j , Z n ), and for each boundary edge ε =  V i , R i , j  of PS , let c e = Bqi ( V i , V j , R i , j ). Then s ∈  S13 (PS ). 1, . . . , |PS |, and σ = e

1 , 2

c iv,r

Proof. Since  S13 (PS ) ⊂ S13 (PS ), we can use the expressions of the B-spline coefficients in (27) to determine the spline s. Let us first have a closer look at the formula

Sm ( P 1 , P 2 , P 3 ) = sm ( P 1 ) + ∇ sm ( P 1 ) · ( P 3 − P 1 ) +

1 6

( P 2 − P 1 ) · H sm ( P 1 ) · (3 P 3 − 2 P 1 − P 2 ),

where sm is the restriction of s to the triangle Tm =  V i , R i , j , Z n  ∈ PS . If s ∈  S13 (PS ) then s ∈ C 2 ( V i ), and so

Sm ( V i , P 2 , P 3 ) = s( V i ) + ∇ s( V i ) · ( P 3 − V i ) +

1 6

( P 2 − V i ) · H s( V i ) · (3 P 3 − 2V i − P 2 ).

From Chen and Liu (2008) we know that any spline s ∈  S13 (PS ) can be uniquely characterized by specifying its value and all derivatives up to order 2 at each vertex V i . Hence, in particular, we can sample them from a polynomial q i ∈ P3 at vertex V i . In such a case,

Sm ( V i , P 2 , P 3 ) = qi ( V i ) + ∇ qi ( V i ) · ( P 3 − V i ) +

1

( P 2 − V i ) · Hqi ( V i ) · (3 P 3 − 2V i − P 2 ) 6 = Bqi ( V i , P 2 , P 3 + ( P 3 − V i ) + ( P 3 − P 2 )).

Then, taking into account

σmt = σe = 12 , the expressions in (27) give

c iv,r = Bqi ( V i , V i , Q iv,r + ( Q iv,r − V i ) + ( Q iv,r − V i )) = Bqi ( V i , V i , 3 Q iv,r − 2V i ), t t t t cm = Bq i ( V i , Z n , Q m + (Q m − V i) + ( Q m − Z n )) = Bqi ( V i , Z n , V j ),

c e = Bqi ( V i , R i , j , Q e + ( Q e − V i ) + ( Q e − R i , j )) = Bqi ( V i , R i , j , V j ), which completes the proof.

2

Note that every spline from the intersection of $ S13 (PS ) and  S13 (PS ) is globally C 2 . This follows from the following three 1 2 observations. First, every spline of  S3 (PS ) is C at each vertex of  and consequently C 2 across each edge of . Second, 1 2 every spline of $ S3 (PS ) is C at each triangle split point of PS and across each edge of PS connecting a triangle split point and an edge split point. Third, every spline of $ S13 (PS ) is C 2 across each edge of PS connecting a vertex S13 (PS ) ∩  2 of  and a triangle split point of PS because the spline is C at the vertex and at the triangle split point. 6.3. C 1 cubic splines on Clough–Tocher triangulations The splines of $ S13 (PS ) considered in Section 6.2 are interesting because they are closely related to the so-called Clough– Tocher splines, named after Clough and Tocher (1965). These are C 1 cubic splines defined on a Clough–Tocher refinement of  that splits every triangle of  into three smaller triangles. It is obtained by choosing a triangle split point inside every triangle of  and by connecting it to the triangle’s vertices. Suppose that CT is the Clough–Tocher refinement of  that uses the same triangle split points as the Powell–Sabin refinement PS of . Then, S13 (CT ) ⊂ $ S13 (PS ) as proven by Speleers (2015). Among all Clough–Tocher splines, a special class is the one containing all splines that have a linear derivative along every edge of  in a direction not parallel to the edge. They form the subspace

% & $ S13 (CT ) = s ∈ S13 (CT ) : D v ε |ε ∈ P1 , ε ∈ E

of dimension 3|V |, which is |E | less than the dimension of S13 (CT ). This allows us to describe the elements of $ S13 (CT ) by a Hermite interpolation problem with values and derivatives provided only at the vertices of . The vector v ε can be any unit vector not parallel to the edge ε ∈ E . In the following, however, we will focus on a particular case of interest

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

21

(see Speleers, 2010b, Example 2.2) where v ε is parallel to the edge of PS connecting the edge split point on ε and the split point inside a triangle with the edge ε . By combining the results of Speleers (2015, Section 5.1) with Theorem 2, we can formulate the following theorem. Theorem 4. Suppose that s ∈ S13 (PS ) is represented in the form (8) with the B-splines determined by choosing 1, . . . , |PS |, and σ = e

1 , 2

 = 1, . . . , |∂ EPS |. For any edge ε =  V i , V j  ∈ E , let

σmt = 12 , m =

βε = λi , j (νi ,1 c iv,1 + νi ,2 c iv,2 + νi ,3 c iv,3 ) + λ j ,i (ν j ,1 c vj,1 + ν j ,2 c vj,2 + ν j ,3 c vj,3 ), where (νh,1 , νh,2 , νh,3 ) are the barycentric coordinates of the point 13 V i + 13 V j + 13 V h with respect to the triangle Qhv =  Q hv,1 , Q hv,2 , Q hv,3  for h = i , j. Denote by CT the Clough–Tocher refinement of  obtained by choosing the same triangle split points as for PS . 1. For each interior edge ε =  V i , V j  of , suppose t t cm = cm , 1 2

t t cm = cm , 3 4

t t t t μn,n (cm + cm ) + μn ,n (cm + cm ) = 2 βε , 1 2 3 4

t t t t where n and n are the indices of the triangles Tn , Tn ∈  sharing the edge ε and cm , cm , cm , cm are the coefficients associated 1 2 3 4 with the triangles Tm1 =  V i , R i , j , Z n , Tm2 =  V j , R i , j , Z n , Tm3 =  V i , R i , j , Z n , Tm4 =  V j , R i , j , Z n , respectively. For each boundary edge ε =  V i , V j  of , suppose

t t cm = cm ,

c e = c e ,

c e + c e = 2βε ,

t t where cm and cm  are the coefficients associated with the triangles Tm =  V i , R i , j , Z n  and Tm =  V j , R i , j , Z n  inside the only triangle Tn ∈  with the edge ε , and c e and c e are the coefficients associated with the boundary edges ε =  V i , R i , j  and

ε =  V j , R i, j , respectively. Then s ∈ S13 (CT ).

2. For each Tm =  V i , R i , j , Z n  of PS , suppose t cm = βε +

1 2

1

(δi ,1 c iv,1 + δi ,2 c iv,2 + δi ,3 c iv,3 ) + (δ j ,1 c vj,1 + δ j ,2 c vj,2 + δ j ,3 c vj,3 ), 2

where (δh,1 , δh,2 , δh,3 ) are the barycentric coordinates of the direction vector 31 Z n − 13 R i , j with respect to the triangle Qhv =  Q hv,1 , Q hv,2 , Q hv,3  for h = i , j. For each boundary edge ε =  V i , R i , j  of PS , suppose c e = βε . Then s ∈ $ S13 (CT ). Proof. From the representation of a Clough–Tocher spline s ∈ S13 (CT ) in terms of the cubic B-spline basis developed by Speleers (2015) and its relation with the cubic B-spline representation (8) given by Theorem 2, we immediately obtain the equalities in item 1 with

  1 1 βε = λ i , j s ( V i ) + ∇ s ( V i ) · ( V j − V i ) + λ j , i s ( V j ) + ∇ s ( V j ) · ( V i − V j ) . 3

3

Furthermore, from the representation of a reduced Clough–Tocher spline s ∈ $ S13 (CT ) in terms of the cubic B-spline basis developed by Speleers (2015), we get t cm = βε +

1 6

1

∇ s( V i ) · ( Z n − R i , j ) + ∇ s( V j ) · ( Z n − R i , j ), 6

The expression in (22) for the coefficients c iv,r completes the proof.

c e = βε .

2

The conditions in Theorem 4 that ensure that s ∈ S13 (PS ) is an element of S13 (CT ) or $ S13 (CT ) can be easily expressed geometrically in terms of the control points. To get a Clough–Tocher spline, one first needs to satisfy the conditions of Theorem 2. This reduces the degrees of freedom of s to one control triangle per vertex of  and one control line per edge of . Then, for every control line the point above the edge must be fixed at a certain height that is determined by the control triangles associated with the endpoints of the edge. This means that the slopes of the control lines can still be freely modified. By fixing the slopes accordingly, one obtains an element of $ S13 (CT ). More details about such geometric control structures can be found in Speleers (2015). 7. Conclusions We have presented a new B-spline representation for the space S13 (PS ) of C 1 cubic splines defined on a triangulation with a Powell–Sabin refinement. We have proven that the cubic B-splines in this representation enjoy interesting properties such as local support, non-negativity, and partition of unity. Their construction is based on lifting particular triangles

22

J. Grošelj, H. Speleers / Computer Aided Geometric Design 57 (2017) 1–22

and line segments from the domain. These triangles and line segments are determined by few parameters which are constrained to enforce non-negativity of the B-splines (Properties 1–3). It is beneficial to choose these parameters such that the corresponding triangles and line segments are as small as possible, both from a stability and a geometric point of view. Furthermore, explicit expressions for the B-spline coefficients of any element of S13 (PS ) have been provided, and we have shown how to compute the Bernstein–Bézier form of such a spline in a stable way. The proposed B-spline representation has great potential in several applications, especially computer aided geometric design and approximation theory. As the B-spline coefficients naturally support a control structure, the representation might be a useful tool for geometric modelling. Thanks to the absence of local super-smoothness, it is possible to construct sequences of nested spline spaces and this allows the development of subdivision schemes in the same spirit as Vanraes et al. (2004); Windmolders and Dierckx (1999). The representation is also suited for the construction of quasi-interpolation schemes, following an approach similar to Manni and Sablonnière (2007); Sbibih et al. (2009). Finally, the cubic B-splines are attractive because they are able to represent both C 1 quadratic Powell–Sabin splines and C 1 cubic Clough–Tocher splines in a unified context. In addition, they are likely to supersede all known cubic B-spline representations developed on Powell–Sabin triangulations (Grošelj and Krajnc, 2016; Lamnii et al., 2014; Speleers, 2015). Acknowledgements This work was partially supported by the MIUR ‘Futuro in Ricerca 2013’ Programme through the project DREAMS (RBFR13FBI3). References Chen, S.K., Liu, H.W., 2008. A bivariate C 1 cubic super spline space on Powell–Sabin triangulation. Comput. Math. Appl. 56, 1395–1401. Clough, R.W., Tocher, J.L., 1965. Finite element stiffness matrices for analysis of plates in bending. In: Conf. on Matrix Methods in Structural Mechanics. Wright–Patterson Air Force Base, Ohio, pp. 515–545. Cohen, E., Lyche, T., Riesenfeld, R.F., 2013. A B-spline-like basis for the Powell–Sabin 12-split based on simplex splines. Math. Comput. 82, 1667–1707. Dierckx, P., 1997. On calculating normalized Powell–Sabin B-splines. Comput. Aided Geom. Des. 15, 61–78. Farin, G., 1985. A modified Clough–Tocher interpolant. Comput. Aided Geom. Des. 2, 19–27. Farin, G., 1986. Triangular Bernstein–Bézier patches. Comput. Aided Geom. Des. 3, 83–127. Grošelj, J., 2016. A normalized representation of super splines of arbitrary degree on Powell–Sabin triangulations. BIT Numer. Math. 56, 1257–1280. Grošelj, J., Krajnc, M., 2016. C 1 cubic splines on Powell–Sabin triangulations. Appl. Math. Comput. 272, 114–126. Klee, V., Laskowski, M.C., 1985. Finding the smallest triangles containing a given convex polygon. J. Algorithms 6, 359–375. Lai, M.J., Schumaker, L.L., 2007. Spline Functions on Triangulations. Cambridge University Press. Lamnii, A., Lamnii, M., Mraoui, H., 2015. A normalized basis for condensed C 1 Powell–Sabin-12 splines. Comput. Aided Geom. Des. 34, 5–20. Lamnii, M., Mraoui, H., Tijini, A., Zidna, A., 2014. A normalized basis for C 1 cubic super spline space on Powell–Sabin triangulation. Math. Comput. Simul. 99, 108–124. Maes, J., Bultheel, A., 2006. Stable multiresolution analysis on triangles for surface compression. Electron. Trans. Numer. Anal. 25, 224–258. Maes, J., Vanraes, E., Dierckx, P., Bultheel, A., 2004. On the stability of normalized Powell–Sabin B-splines. J. Comput. Appl. Math. 170, 181–196. Mann, S., 1999. Cubic precision Clough–Tocher interpolation. Comput. Aided Geom. Des. 16, 85–88. Manni, C., Sablonnière, P., 2007. Quadratic spline quasi-interpolants on Powell–Sabin partitions. Adv. Comput. Math. 26, 283–304. O’Rourke, J., Aggarwal, A., Maddila, S., Baldwin, M., 1986. An optimal algorithm for finding minimal enclosing triangles. J. Algorithms 7, 258–269. Powell, M.J.D., Sabin, M.A., 1977. Piecewise quadratic approximations on triangles. ACM Trans. Math. Softw. 3, 316–325. Sbibih, D., Serghini, A., Tijini, A., 2009. Polar forms and quadratic spline quasi-interpolants on Powell–Sabin partitions. Appl. Numer. Math. 59, 938–958. Sbibih, D., Serghini, A., Tijini, A., 2012. Normalized trivariate B-splines on Worsey–Piper split and quasi-interpolants. BIT Numer. Math. 52, 221–249. Seidel, H.P., 1993. An introduction to polar forms. IEEE Comput. Graph. Appl. 13, 38–46. Speleers, H., 2010a. A normalized basis for quintic Powell–Sabin splines. Comput. Aided Geom. Des. 27, 438–457. Speleers, H., 2010b. A normalized basis for reduced Clough–Tocher splines. Comput. Aided Geom. Des. 27, 700–712. Speleers, H., 2011. On multivariate polynomials in Bernstein–Bézier form and tensor algebra. J. Comput. Appl. Math. 236, 589–599. Speleers, H., 2013a. Construction of normalized B-splines for a family of smooth spline spaces over Powell–Sabin triangulations. Constr. Approx. 37, 41–72. Speleers, H., 2013b. Multivariate normalized Powell–Sabin B-splines and quasi-interpolants. Comput. Aided Geom. Des. 30, 2–19. Speleers, H., 2015. A new B-spline representation for cubic splines over Powell–Sabin triangulations. Comput. Aided Geom. Des. 37, 42–56. Speleers, H., Dierckx, P., Vandewalle, S., 2006. Numerical solution of partial differential equations with Powell–Sabin splines. J. Comput. Appl. Math. 189, 643–659. Speleers, H., Manni, C., Pelosi, F., Sampoli, M.L., 2012. Isogeometric analysis with Powell–Sabin splines for advection–diffusion–reaction problems. Comput. Methods Appl. Mech. Eng. 221–222, 132–148. Vanraes, E., Windmolders, J., Bultheel, A., Dierckx, P., 2004. Automatic construction of control triangles for subdivided Powel–Sabin splines. Comput. Aided Geom. Des. 21, 671–682. Windmolders, J., Dierckx, P., 1999. Subdivision of uniform Powell–Sabin splines. Comput. Aided Geom. Des. 16, 301–315.