- Email: [email protected]
C1 cubic splines on Powell–Sabin triangulations
Applied Mathematics and Computation 272 (2016) 114–126
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
C 1 cubic splines on Powell–Sabin triangulations Jan Grošelj a,∗, Marjeta Krajnc a,b a b
FMF, University of Ljubljana, Jadranska 19, Ljubljana, Slovenia IMFM, Jadranska 19, Ljubljana, Slovenia
a r t i c l e
i n f o
Article history: Available online 27 July 2015 MSC: 41A05 41A15 65D07 65D17
a b s t r a c t A bivariate C 1 cubic spline space on a triangulation with Powell–Sabin refinement which extends the well-known C 1 quadratic spline space and has a nested structure is introduced. A construction of a locally supported basis that forms a partition of unity is presented based on choosing particular triangles and line segments in the domain. Further, it is shown how these objects can be determined in order to obtain nonnegative basis functions under a natural restriction on the Powell–Sabin refinement. Geometrically intuitive B-spline representation is proposed which makes these splines a useful tool for CAGD applications.
Keywords: Nested spline spaces Powell–Sabin triangulation Normalized B-splines Control structure
© 2015 Elsevier Inc. All rights reserved.
1. Introduction Spline functions are one of the most efficient and widely used tools in approximation theory and computer aided geometric design. Univariate spline spaces are very well examined, but a step to a bivariate case brings many difficulties and many problems are still unsolved. A bivariate spline is a function of a particular smoothness, composed of polynomial pieces defined on some partition that can be rectangular or triangular. Although tensor product splines on rectangular domains have recently attract much attention in isogeometric analysis, splines on triangulations are still widely used in function approximation, finite element method, for design purposes, etc. The bivariate spline spaces on triangulations have been intensively studied by many researchers (see [1] and the references therein), and it turns out that even the dimension determination is an extremely difficult problem since it heavily depends on the geometry of the triangulation. To be able to determine the dimension and to construct locally supported basis functions, one can use the so called macro elements [1]. Namely, every triangle of the initial triangulation is split into several smaller triangles, and a spline space on the refined triangulation is considered. The two refinements that are most commonly used are the Powell–Sabin and the Clough–Tocher split. A Powell–Sabin refinement splits every triangle of a triangulation into six smaller triangles. Splines on a Powell–Sabin refinement have first been studied in [2]. It has been shown therein that a C 1 quadratic spline on the refined triangulation is uniquely specified by the values and the first order derivatives at the vertices of the initial triangulation. The dual basis functions of the associated spline space have local supports, and it was proven in [3] that by slightly changing their interpolation values, nonnegative basis functions that form a partition of unity may be derived. A generalization of quadratic Powell–Sabin splines to higher degrees and higher orders of smoothness led to the so called Powell–Sabin super spline spaces, where the smoothness at some ∗
Corresponding author. E-mail address: [email protected] (J. Grošelj).
http://dx.doi.org/10.1016/j.amc.2015.07.013 0096-3003/© 2015 Elsevier Inc. All rights reserved.
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
115
vertices and across some edges is increased. In particular, C r super splines of degree 3r − 1 were presented in [4] together with a normalized representation. Cubic splines have also been considered in the literature, especially in terms of their application to the theory of quasiinterpolants. Based on the representation of C 2 cubic splines on the uniform Powell–Sabin refinement of a rectangular domain with pairs of translated multi-box splines, quasi-interpolants of differential and discrete type were derived in [5,6]. For arbitrary polygonal domains and general Powell–Sabin refinements, a C 1 cubic super spline space with C 2 smoothness conditions at the vertices of the initial triangulation was proposed in [7]. Its normalized representation with applications to quasi-interpolants was studied extensively in [8–10], thus generalizing the results for the quadratic case (see [11–13]). Increasing the smoothness at some points and across some edges of the triangulation makes the dimension counting problem easier and also simplifies the construction of the basis functions. But the super spline formulation has at least one important drawback, namely the super spline spaces cannot be nested. The sequence of spaces S0 , S1 , S2 , . . . defined on the triangulations 0 , 1 , 2 , . . . , such that n is a refinement of n−1 for each n, is a nested sequence of spaces if Sn−1 ⊂ Sn . Clearly this cannot hold for super spline spaces. The nested structure enables the use of (local) subdivision and allows the construction of hierarchical spline spaces, which are of a great importance in e.g. solving partial differential equations. For these applications with quadratic Powell–Sabin splines see [14–18]. In order to extend the number of Powell–Sabin spline spaces with the nested structure, a space of C 1 cubic splines on a Powell–Sabin refinement is introduced in this paper, which in comparison to the quadratic case allows more flexibility for the approximation and design purposes since more degrees of freedom are provided. Each element of the presented space is uniquely defined by prescribing the values and the first order derivatives at the vertices of the initial triangulation and at the triangle split points inside the initial triangles. Additionally, the directional derivatives at three chosen points in each initial triangle, and the values and the directional derivatives at the split point on each initial boundary edge are interpolated. Particularly chosen interpolation problems give locally supported basis functions which form a partition of unity. It is shown also that the basis functions are nonnegative under some natural restriction on the geometry of the Powell–Sabin refinement. In order to make the presented splines a useful and a user friendly design tool, a B-spline representation and a control structure that intuitively affects the shape of the spline are derived, too. The paper is organized as follows. In Section 2 some preliminaries about polynomials on triangles and splines on triangulations are presented. In Section 3 the C 1 cubic spline space on a Powell–Sabin triangulation is introduced, and the Bézier representation of its elements is derived in dependence of an interpolation problem. Section 4 gives a construction of a locally supported basis which forms a partition of unity. Geometric conditions on the construction, which ensure that the basis functions are nonnegative, are discussed. In the end, the B-spline representation is derived, and its application to geometric design is demonstrated. 2. Preliminaries 2.1. Bivariate polynomials on triangles Consider a triangle T V 1 , V 2 , V 3 with vertices Vi , i = 1, 2, 3. Each point P ∈ R2 can be uniquely expressed as an affine combination
P = τ1 V 1 + τ2 V 2 + τ3 V 3 ,
τ1 + τ2 + τ3 = 1.
The coefficients τ = (τ1 , τ2 , τ3 ) are called the barycentric coordinates of the point P with respect to the triangle T . If the point P lies inside the triangle, its barycentric coordinates are all positive. Let Pd denote the space of polynomials
p(x, y) =
ai j xi y j ,
ai j ∈ R,
0≤i+ j≤d
in two variables of total degree less than or equal to d ∈ N. Each bivariate polynomial p ∈ Pd can be expressed in the Bernstein– Bézier representation with respect to the triangle T as
p(P ) ≡ p(τ) =
bi Bdi (τ),
i = (i1 , i2 , i3 ),
|i| = i1 + i2 + i3 ,
(1)
|i|=d
where Bdi are the Bernstein basis polynomials of degree d, defined as
Bdi (τ) =
d! i1 ! i2 ! i3 !
τ1i1 τ2i2 τ3i3 , |i| = d,
and bi , |i| = d, are called the Bézier ordinates of p. The domain point Di that corresponds to the Bézier ordinate bi is determined
by the barycentric coordinates ( d1 , d2 , d3 ), and a pair (Di , bi ) is called a control point of p. The de Casteljau algorithm allows us to compute the values of p from (1) in a stable and efficient way. = ( 2 , η 3 ) be the barycentric coordinates of points Q = (xQ , yQ ) and Q 1 , η η = (η yQ ) with respect xQ , Let η = (η1 , η2 , η3 ) and with respect to T are equal to to T . The barycentric directions of the vector v := Q − Q i
i
i
1 , η2 − η 2 , η3 − η 3 ) ν = (ν1 , ν2 , ν3 ) = (η1 − η
116
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
and satisfy ν1 + ν2 + ν3 = 0. The directional derivative of p at the point P in the direction of v is defined as
Dv p(P ) :=
d p(P + t v) = Dx p(P )(xQ − xQ ) + Dy p(P )(yQ − yQ ) dt
and equals
Dv p(P ) = d
|i|=d−1
(ν1 bi+(1,0,0) + ν2 bi+(0,1,0) + ν3 bi+(0,0,1) )Bd−1 (τ). i
(2)
(3)
2.2. Spline spaces on triangulations Let ⊂ R2 be a simply connected domain with a polygonal boundary, and let = {Tk }k be its arbitrary regular triangulation. The sets of vertices and edges of are denoted by V = {V i }i and E = {ε } , respectively. The sets are assumed to be finite, and their elements are indexed with the indices from 1 to the cardinality of a set, which is denoted by | · |. The edge ε is sometimes written as ε Vi , Vj to indicate that its endpoints are the vertices Vi and Vj . The edge lying on the boundary of the domain is a boundary edge. The set of all boundary edges is denoted by ∂ E. The subset of all triangles in sharing the vertex Vi is called the molecule of Vi and is denoted by M (Vi ). A space of polynomial splines of degree d ∈ N and smoothness r ∈ N0 on a triangulation is defined as
Sdr () := {s ∈ C r (); s|T ∈ Pd , T ∈ }. It consists of splines composed of polynomials of degree d, defined on triangles of , that join together with C r smoothness. By using (1) each of the polynomial segments can be locally represented in the Bernstein–Bézier form and is completely defined by prescribing its control points. The union of the Bézier representations (1) of all the polynomial segments gives the Bézier spline representation. Smoothness constraints imply certain relations between the control points of the neighbouring polynomials. In particular, let p and p be two polynomials of total degree d written in the form (1) on T V 1 , V 2 , V 3 and T V 4 , V 3 , V 2 with Bézier ordinates bi and b i , respectively. Suppose that V4 lies outside of T and has barycentric coordinates σ = (σ1 , σ2 , σ3 ) with respect to T . Then p and p are C 1 continuous across the common edge V2 , V3 if and only if
b 0,i2 ,i3 = b0,i3 ,i2 , b
1,i2 ,i3
i2 + i3 = d,
(4a)
= σ1 b1,i3 ,i2 + σ2 b0,i3 +1,i2 + σ3 b0,i3 ,i2 +1 ,
i2 + i3 = d − 1.
(4b)
The proof and the higher order smoothness conditions can be found in [1]. 2.3. Powell–Sabin refinement A Powell–Sabin refinement ∗ of the initial triangulation , introduced in [2], splits every triangle into six smaller triangles in the following way. First, for every triangle Tk ∈ , k = 1, 2, . . . , n , a point Zk in the interior of Tk is chosen and is connected to the vertices of the triangle. The newly inserted points can be chosen as the incenters of the corresponding triangles or, more generally, can be chosen in such a way that for every two triangles Tk and Tk with a common edge Vi , Vj the intersection point Rij of the line segments Z k , Z k and Vi , Vj lies in the interior of Vi , Vj . The point Rij is then connected to Zk and Z k . Additional points are chosen on the boundary edges of the triangulation and are connected to the interior split points of the associated triangles. An example of a Powell–Sabin refinement of two neighbouring triangles is shown in Fig. 1. The following notation will further be needed. For every triangle split point Z k ∈ Tk V 1 , V 2 , V 3 , let ξ k, 1 , ξ k, 2 , ξ k, 3 be such that
Z k = ξk,1V 1 + ξk,2V 2 + ξk,3V 3 ,
ξk,1 + ξk,2 + ξk,3 = 1.
For every edge split point Rij ∈ Vi , Vj , let λij and λji be such that
Ri j = λi jV i + λ jiV j ,
λi j + λ ji = 1.
Furthermore, the cartesian coordinates of the points Vi , Rij , and Zk are denoted by (xi , yi ), (x∗i j , y∗i j ), and (x∗k , y∗k ), respectively. Note that the points Rij and Rji coincide. 3. A space of C 1 cubic splines Let be an arbitrary triangulation and ∗ its Powell–Sabin refinement. Consider the spline space S31 (∗ ) of C 1 cubic splines on ∗ . The aim of this section is to provide an interpolation problem that uniquely describes an element of S31 (∗ ). Suppose that s ∈ S31 (∗ ). For a vertex V i ∈ V, let fi , fi, x , and fi, y be the interpolation values such that
s(V i ) = fi ,
Dx s(V i ) = fi,x ,
Dy s(V i ) = fi,y .
(5a)
For a triangle Tk V 1 , V 2 , V 3 ∈ , let gk , gk, x , and gk, y be the interpolation values such that
s(Z k ) = gk ,
Dx s(Z k ) = gk,x ,
Dy s(Z k ) = gk,y .
(5b)
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
117
Fig. 1. A Powell–Sabin refinement of two neighbouring triangles.
Fig. 2. Domain points of a cubic spline on a macro-triangle Tk V 1 , V 2 , V 3 .
Moreover, for a triangle Tk ∈ and each of its vertices Vi , choose a point Wk, i lying on the edge Vi , Zk , i.e.,
W k,i = ωk,iV i + (1 − ωk,i )Z k ,
0 < ωk,i < 1.
Additionally, choose any unit vector wk, i nonparallel to Z k − V i , and let gwk,i be the interpolation values such that
Dwk,i s(W k,i ) = gwk,i .
(5c)
Finally, for a boundary edge V i , V j ∈ ∂ E, let rij be a unit vector parallel to V j − V i , and let hij and hri j be the interpolation values such that
s(Ri j ) = hi j ,
Dri j s(Ri j ) = hri j .
(5d) S31 (∗ ).
The conditions in (5) specify an interpolation problem for In what follows, it is proven that the interpolation problem admits a unique solution. Let Tk V 1 , V 2 , V 3 be an arbitrary triangle. It will be assumed hereafter that the domain points of a spline s ∈ S31 ({Tk }∗ ) on Tk are denoted as in Fig. 2. The points P i j = P ji and T i j = T ji , i, j ∈ {1, 2, 3}, will be identified to simplify the notation. Moreover, the Bézier ordinate of s associated with a domain point will be denoted by the lower case letter. For example, the Bézier ordinate associated with the domain point O12 will be denoted by o12 .
118
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
Lemma 1. Let Tk V 1 , V 2 , V 3 be an arbitrary triangle. There exists a unique spline s ∈ S31 ({Tk }∗ ) satisfying the conditions in (5). Proof. Suppose that an arbitrary set of interpolation values prescribing the interpolation conditions in (5) is given. In what follows, the Bézier ordinates of the spline s from S31 ({Tk }∗ ) that satisfies these conditions are derived, which proves both, the existence and the uniqueness of s. By (2) and (3) the conditions (5a) imply vi = fi and
1 ( fi,x (x∗k − xi ) + fi,y (y∗k − yi )), 3 1 ui j = fi + ( fi,x (x∗i j − xi ) + fi,y (y∗i j − yi )) 3 mi = fi +
(6)
for every {i, j}⊆{1, 2, 3}. Analogously by (5b), zk = gk and
1 (g (xi − x∗k ) + gk,y (yi − y∗k )), 3 k,x 1 ti j = gk + (gk,x (x∗i j − x∗k ) + gk,y (y∗i j − y∗k )). 3 ni = gk +
(7)
Furthermore, let (wij, 1 , wij, 2 , wij, 3 ) denote the barycentric directions of wk, i with respect to the triangle Vi , Zk , Rij . By (3) and (5c),
wi j,2 wi j,2 wi j,1 wi j,1 ωk,i oi j = − vi + mi + ui j − mi + ni wi j,3 wi j,3 wi j,3 2(1 − ωk,i ) wi j,3 wi j,2 1 − ωk,i wi j,1 1 1 − ni + z + ti j + gw . 2ωk,i wi j,3 wi j,3 k 6ωk,i (1 − ωk,i ) wi j,3 k,i
(8)
By the C 1 smoothness conditions (4b) across the edges Rij , Zk , it follows:
pi j = λi j oi j + λ ji o ji .
(9)
Finally, by (2) and (3) the conditions (5d) imply ri j = hi j and
li j = hi j +
1 hr (V i − Ri j )·ri j , 3 ij
(10)
where · denotes the standard scalar product in R2 . The conditions in (4a) are clearly satisfied everywhere by the notation. The conditions (4b) across the edges Rij , Zk are satisfied by (7), (9), (10), and the fulfillment of the same conditions across Vi , Zk can be easily verified by following [1, p. 35]. Since all the Bézier ordinates of s are now determined and all smoothness conditions are satisfied, the proof is completed. Lemma 2. Let Tk V 1 , V 2 , V 3 and Tk V 4 , V 3 , V 2 be two triangles with the common edge V2 , V3 . There exists a unique spline s ∈ S31 ({Tk , Tk }∗ ) satisfying the conditions in (5). Proof. All the Bézier ordinates of the spline s may be expressed in the same way as in Lemma 1, except for the values r23 , l23 , and l32 . It needs to be checked that they can be uniquely determined, and that s satisfies the smoothness conditions along and across the edge V2 , V3 .
Suppose that R23 = μ23,k Z k + μ23,k Z k , μ23,k + μ23,k = 1. Using (6), it may be easily verified that ui j = μ23,k mki + μ32,k mki , 1 {i, j} = {2, 3}, where the upper index denotes the triangle that the ordinate belongs to. The other three C smoothness conditions (4b) across V2 , V3 imply
r23 = μ23,k pk23 + μ23,k pk23 ,
l23 = μ23,k ok23 + μ23,k ok23 , l32 = μ
k 23,k o32
+μ
(11)
k
23,k o32 .
It follows from (9) and (11) that
r23 = μ23,k pk23 + μ23,k pk23
= μ23,k (λ23 ok23 + λ32 ok32 ) + μ23,k (λ23 ok23 + λ32 ok32 )
= λ23 (μ23,k ok23 + μ23,k ok23 ) + λ32 (μ23,k ok32 + μ23,k ok32 ) = λ23 l23 + λ32 l32 , which means that the C 1 smoothness condition at R23 in the direction of V 3 − V 2 is satisfied, and hence the lemma is proven. Theorem 1. Let be an arbitrary triangulation. There exists a unique spline s ∈ S31 (∗ ) satisfying the conditions in (5). Proof. The proof follows directly from Lemmas 1 and 2.
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
119
Fig. 3. A choice of a triangle Qvi Q vi,1 , Q vi,2 , Q vi,3 (left) and the associated basis function Bvi,1 (right).
According to Theorem 1, the dimension of S31 (∗ ) can be determined by a simple count of the interpolation conditions in (5). Namely, (5a) gives three conditions per vertex, each of (5b) and (5c) three conditions per triangle, and (5d) two conditions per boundary edge. Corollary 1. The dimension of S31 (∗ ) is equal to 3|V | + 6|| + 2|∂ E |. In order to show that the considered spline space has a nested structure, a suitable refinement R of has to be found such that a Powell–Sabin refinement ∗R of R may be constructed from the existing Powell–Sabin refinement ∗ . When this is the case, it follows trivially that S31 (∗ ) ⊆ S31 (∗R ) since the restriction of s ∈ S31 (∗ ) to each T ∗ ∈ ∗ is a cubic polynomial, which can be subdivided according to the refinement of T ∗ implied by ∗R , and all smoothness conditions of S31 (∗R ) are fulfilled for the spline obtained in this way. One obvious choice for R is ∗ , and a Powell–Sabin refinement of R can be chosen freely. However, with such a refinement the number of triangles grows very fast, which may not be desirable. A common refinement of a triangulation is a dyadic split that divides every triangle into four smaller triangles and is obtained by adding an additional point on every edge of the triangulation. For a triangle Tk V 1 , V 2 , V 3 ∈ , the candidates for additional points would be R12 , this approach yields an incompatible refinement when the triangle R12 , R23 , R31 does not R23 , R31 , but as illustrated in [18], √ contain Zk . As a better√practice the 3-subdivision rule is proposed therein, which results in a triadic split if one applies it twice. The advantage of the 3-rule is also that it allows a local refinement as proposed in [15] for the quadratic case. 4. A normalized basis In this section a construction of a basis for a space S31 (∗ ) is presented. It introduces basis functions of three kinds, each of them is associated with a vertex, a triangle or a boundary edge of the initial triangulation. More precisely, the aim is to represent an arbitrary spline s ∈ S31 (∗ ) as a linear combination
s=
|V | 3
v Bv + ci,d i,d
i=1 d=1
|| 3 2
ctk,i,d Btk,i,d +
k=1 i=1 d=1
|∂ E | 2 =1 d=1
e c,d Be,d ,
(12)
where the basis functions Bvi,d are associated with the vertex V i ∈ V, the basis functions Btk,i,d are associated with the triangle Tk V 1 , V 2 , V 3 ∈ , and the basis functions Be,d are associated with the boundary edge ε V i , V j ∈ ∂ E. Beside the construction it will be shown that it is possible to define the basis functions in a way that they have local supports and form a partition of unity. It will also be proven that they are nonnegative under a certain assumption on the geometry of the Powell–Sabin refinement. 4.1. Basis functions Inspired by the ideas of a geometric construction of a normalized basis for S21 (∗ ) presented in [3], a basis for S31 (∗ ) may be constructed as follows. For a demonstration, consider Figs. 3–5. v , α v , α v ) be the barycentric coordinates of V 1. For every V i ∈ V, choose a triangle Qvi Q vi,1 , Q vi,2 , Q vi,3 in the plane. Let (αi,1 i i,2 i,3 v , β v , β v ) and (γ v , γ v , γ v ) be the barycentric directions of (1, 0) and (0, 1) with respect to with respect to Qvi , and let (βi,1 i,2 i,3 i,1 i,2 i,3
Qvi . Then Bvi,d , d = 1, 2, 3, are the splines from S31 (∗ ) uniquely defined by setting the interpolation values (5a) to v , Bvi,d (V i ) = αi,d
v , Dx Bvi,d (V i ) = βi,d
v , Dy Bvi,d (V i ) = γi,d
and all the other interpolation values in (5) to zero.
120
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
Fig. 4. A choice of a triangle Qtk Q tk,1 , Q tk,2 , Q tk,3 and line segments Ltk,i Q tk,i,1 , Q tk,i,2 , i = 1, 2, 3 (left), and the associated basis function Btk,3,2 (right).
Fig. 5. A choice of a line segment Le Q e,1 , Q e,2 (left) and the associated basis function Be,2 (right).
t , α t , α t ) be the barycentric co2. For every Tk V 1 , V 2 , V 3 ∈ , choose a triangle Qtk Q tk,1 , Q tk,2 , Q tk,3 in the plane. Let (αk,1 k,2 k,3
t , β t , β t ) and (γ t , γ t , γ t ) be the barycentric directions of (1, 0) and ordinates of Zk with respect to Qtk , and let (βk,1 k,2 k,3 k,1 k,2 k,3
(0, 1) with respect to Qtk . Additionally, choose line segments Ltk,i Q tk,i,1 , Q tk,i,2 , i = 1, 2, 3, on the lines through Qk, i with the
t t direction wk, i , respectively. Let (σk,i,1 , σk,i,2 ) be the barycentric directions of wk, i with respect to Ltk,i . Then Btk,i,d , i = 1, 2, 3,
d = 1, 2, are the splines from S31 (∗ ) uniquely defined by setting the interpolation values (5b) to
Btk,i,d (Z k ) =
1 2
t t t αk,i , Dx Btk,i,d (Z k ) = 12 βk,i , Dy Btk,i,d (Z k ) = 12 γk,i ,
and (5c) to
Dwk,i Btk,i,d (W k, ) =
t δ σk,i,d , = 1, 2, 3,
1 2 i,
with all the other interpolation values in (5) being equal to zero. Note that δ i, denotes the Kronecker delta. e , α e ) be the 3. For every ε V i , V j ∈ ∂ E, choose a line segment Le Q e,1 , Q e,2 on the line containing the edge ε . Let (α,1 ,2 e , β e ) be the barycentric directions of r with respect to Le . barycentric coordinates of Rij with respect to Le , and let (β,1 ij ,2 Then Be,d , d = 1, 2, are the splines from S31 (∗ ) uniquely defined by setting the interpolation values (5d) to e Be,d (Ri j ) = α,d ,
e Dri j Be,d (Ri j ) = β,d ,
and all the other interpolation values in (5) to zero.
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
121
It needs to be proven that the presented splines form a basis for S31 (∗ ) which admits the properties of a local support and a v , β v , γ v , d = 1, 2, 3, satisfy the equality partition of unity. Note first that for every i = 1, 2, . . . , |V | the parameters αi,d i,d i,d
⎡
1 v ⎣Xi,1 v Yi,1
1 v Xi,2 v Yi,2
⎤⎡
v αi,1 1 v ⎦⎣α v Xi,3 i,2 v v Yi,3 αi,3
v βi,1 v βi,2 v βi,3
⎤
v γi,1 1 v ⎦= x γi,2 i v yi γi,3
0 1 0
0 0 , 1
(13)
v , Y v ) denote the cartesian coordinates of Q v , d = 1, 2, 3. Similarly, for every k = 1, 2, . . . , || the parameters α t , where (Xi,d i,d i,d k,i
t , γ t , i = 1, 2, 3, satisfy βk,i k,i
⎡
1 t ⎣Xk,1 t Yk,1
1 t Xk,2 t Yk,2
⎤⎡
t βk,1 t βk,2 t βk,3
t 1 αk,1 t ⎦⎣ t Xk,3 αk,2 t t Yk,3 αk,3
⎤ ⎡ t γk,1 1 t ⎦ γk,2 = ⎣x∗k t γk,3 y∗k
0 1 0
⎤
0 0⎦ , 1
(14)
t , Y t ) are the cartesian coordinates of Q t , i = 1, 2, 3. where (Xk,i k,i k,i
Theorem 2. The splines Bvi,d , Btk,i,d , Be,d that appear in (12) form a basis of S31 (∗ ). Proof. Since the number of the splines is equal to the dimension of the spline space S31 (∗ ), it is sufficient to prove that they are linearly independent. Suppose that the expression in (12) is equal to 0. The evaluation of s, Dx s, and Dy s at V i ∈ V results in the system of equations 3
v αv = ci,d i,d
d=1
3
v βv = ci,d i,d
3
d=1
v γ v = 0, ci,d i,d
d=1
v = 0, d = 1, 2, 3. Let T V , V , V ∈ . Similarly as above, the evaluation which has a unique solution by (13). This means that ci,d 1 2 3 k of s, Dx s, and Dy s at Zk gives the system of equations 3
ctk,i,1 + ctk,i,2
i=1
t = αk,i
3
ctk,i,1 + ctk,i,2
i=1
t = βk,i
3
ctk,i,1 + ctk,i,2
t = 0. γk,i
i=1
The evaluation of Dwk,i s, i = 1, 2, 3, at Wk, i extends the system with the equations t t ctk,i,1 σk,i,1 + ctk,i,2 σk,i,2 = 0,
i = 1, 2, 3.
Using (14), one may easily verify that the system of six equations has a unique solution. This proves that ctk,i,d = 0, i = 1, 2, 3, d = 1, 2. Finally, let ε V i , V j ∈ ∂ E. Evaluation of s and Dri j s at Rij results in 2 d=1
e e c,d α,d =
2
e e c,d β,d = 0,
d=1
e = 0, d = 1, 2. from where it follows that c,d
Theorem 3. The splines Bvi,d , Btk,i,d , Be,d that appear in (12) have local supports. Proof. Consider the Bézier ordinates of Bvi,d on a triangle of outside M (Vi ). Since the spline on this triangle is determined with zero interpolation values, all ordinates are zero as can be verified by following the proofs of Lemmas 1 and 2. This shows that the support of Bvi,d is contained in M (Vi ). A similar deduction shows that on a triangle of that does not share any edge with Tk ∈ all the Bézier ordinates of Btk,i,d are zero, and the support of this basis function is contained in the area spanned by Tk and its neighbours. The only nonzero Bézier ordinates of Be,d correspond to the domain points that are located on the edge ε V i , V j ∈ ∂ E, and so the support of Be,d is contained in a single triangle of . Theorem 4. The splines Bvi,d , Btk,i,d , Be,d that appear in (12) form a partition of unity. Proof. It is easy to verify that the sum of the basis functions satisfies the interpolation conditions for the constant 1. By Theorem 1, the proof follows. 4.2. Nonnegativity of the basis functions v , ct e , In order to view the spline s ∈ S31 (∗ ) expressed in the form (12) as a convex combination of the coefficients ci,d , c,d k,i,d
the basis functions Bvi,d , Btk,i,d , Be,d have to be nonnegative. This property imposes additional constraints on the parameters that determine the basis functions. Unfortunately, these constraints cannot be satisfied in general, and some additional assumptions on the geometry of the Powell–Sabin refinement are needed. In what follows, it will be assumed that for every triangle
122
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
Tk V 1 , V 2 , V 3 ∈ the interior split point Zk lies in the interior of the triangle R12 , R23 , R31 . This requirement is sufficient but not necessary for the construction of nonnegative basis functions. Let Tk V 1 , V 2 , V 3 be a triangle in . Suppose hereafter that the vector wk, i , i = 1, 2, 3, specified in the interpolation problem (5), is a unit vector parallel to Ri, j1 − Ri, j2 , { j1 , j2 } = {1, 2, 3}\{i}. Denote by
Ri = (1 − ρi )V i + ρi Z k the intersection point of Vi , Zk and Ri, j1 , Ri, j2 . Note that by the assumption, 0 < ρ i < 1. The construction of nonnegative basis functions is based on imposing geometric conditions, expressed as the requirements that the triangles Qvi and Qtk and the line segments Ltk,i and Le contain certain sets of points. Lemma 3. Let V1 be an arbitrary vertex of . If the triangle Qv1 contains the points V1 , U12 , U13 , M1 ,
M 12 = M 1 + M 13 = M 1 +
ωk,1
6(1 − ωk,1 )(1 − ρ1 )
ωk,1
6(1 − ωk,1 )(1 − ρ1 )
(R1 − R12 ), (R1 − R13 )
for every Tk V 1 , V 2 , V 3 ∈ M (V 1 ), then the basis functions Bv1,d , d = 1, 2, 3, are nonnegative. Proof. By the nonnegativity of the Bernstein basis polynomials, it is sufficient to check that the Bézier ordinates of Bv1,d , d = 1, 2, 3, are nonnegative. By Theorem 3, only the Bézier ordinates corresponding to the domain points in M (V1 ) are nonzero. Furthermore, as follows from the proofs of Lemmas 1 and 2, the nonnegativity of the ordinates needs to be verified only for v1 , u12 , u13 , m1 , o12 , and o13 associated with an arbitrary triangle Tk V 1 , V 2 , V 3 in M (V1 ). To show this, recall the property that the barycentric coordinates of a point with respect to a triangle are nonnegative if and only if the point lies inside the triangle. Consider the transformation
⎡ 1 π1 v π2 = ⎣X1,1 v π3 Y1,1
⎤−1 ⎡
1 v X1,2 v Y1,2
1 v ⎦ X1,3 v Y1,3
1 ⎣ x1 y1
1 x2 y2
⎤ τ1 τ2 τ3
1 x3 ⎦ y3
between the barycentric coordinates (τ 1 , τ 2 , τ 3 ) with respect to Tk and the barycentric coordinates (π 1 , π 2 , π 3 ) with respect to Qv1 . By (13), it may be rewritten as
⎡α v π1 1,1 v π2 = ⎣α1,2 v π3 α1,3
v + ηv α1,1 1,1,2 v + ηv α1,2 1,2,2 v + ηv α1,3 1,3,2
⎤ v + ηv α1,1 τ1 1,1,3 v + ηv ⎦ τ2 , α1,2 1,2,3 v + ηv τ3 α1,3 1,3,3
(15)
v v (x − x ) + γ v (y − y ), d = 1, 2, 3, j = 2, 3. The barycentric coordinates of V , U , U , and M with respect = β1,d where η1,d, 1 1 1 12 13 1 j j j 1,d to Tk are given by
(1, 0, 0),
2 + λ12 λ21 , ,0 , 3 3
2 + λ13 λ31 , 0, 3 3
,
2 + ξk,1 ξk,2 ξk,3 , , 3 3 3
.
It is easy to see that by applying the transformation (15) one obtains the ordinates v1 , u12 , u13 , and m1 of Bv1,d , d = 1, 2, 3, corresponding to Tk . It remains to check that the values of o12 and o13 are nonnegative. Using the same approach as above, one can see that these values can be expressed as the barycentric coordinates of M12 and M13 with respect to Qv1 , multiplied by 1 − ρ1 . By the assumption, this factor is positive, and the result follows. Lemma 4. For a triangle Tk V 1 , V 2 , V 3 ∈ , choose Qtk Q tk,1 , Q tk,2 , Q tk,3 to be a triangle containing the points
Ni j = Ni +
1 − ωk,i (Ri − Ri j ), 6ωk,i ρi
i ∈ {1, 2, 3}, j ∈ {1, 2, 3}\{i},
such that for each i = 1, 2, 3 none of the points N i, j1 and N i, j2 lies on the line segment Q tk, j , Q tk, j , { j1 , j2 } = {1, 2, 3}\{i}. Suppose 1
2
that χ ij is the barycentric coordinate of Nij with respect to the triangle Qtk corresponding to Q tk,i , and let χi = χi, j1 + χi, j2 , { j1 , j2 } =
{1, 2, 3}\{i}. For each i = 1, 2, 3, choose Ltk,i Q tk,i,1 , Q tk,i,2 to be a line segment containing the points K i j,d = Q tk,i +
( − 1)d (Ri − Ri j ), 6ωk,i (1 − ωk,i )ρi χi χi, j
j ∈ {1, 2, 3}\{i}, d ∈ {1, 2},
χi, j
with the property that Q tk,i = χ 1 Q tk,i,1 + χ 2 Q tk,i,2 . Then the basis functions Btk,i,d , i = 1, 2, 3, d = 1, 2, are nonnegative. i i Proof. By a similar deduction as in Lemma 3, one finds out that the proof follows if and only if the ordinates zk , n1 , n2 , n3 , o12 , o13 , o21 , o23 , o32 , o31 of Btk,i,d corresponding to the domain points in Tk V 1 , V 2 , V 3 are nonnegative. The transformation between
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
123
the barycentric coordinates (τ 1 , τ 2 , τ 3 ) with respect to Tk and the barycentric coordinates (π 1 , π 2 , π 3 ) with respect to Qtk is now given by
⎡αt + ηt π1 k,1 k,1,1 t t + ηk,2,1 π2 = ⎣αk,2 t t π3 αk,3 + ηk,3,1
t t αk,1 + ηk,1,2 t t αk,2 + ηk,2,2 t t αk,3 + ηk,3,2
⎤ t t αk,1 + ηk,1,3 τ1 t t ⎦ τ2 , αk,2 + ηk,2,3 t t τ3 αk,3 + ηk,3,3
(16)
t t (x − x∗ ) + γ t (y − y∗ ), i, j = 1, 2, 3. By applying the transformation (16) to (ξ = βk,i where ηk,i, j j k, 1 , ξ k, 2 , ξ k, 3 ) and k k j k,i
1 + 2ξk,1 2ξk,2 2ξk,3 , , 3 3 3
,
2ξk,1 1 + 2ξk,2 2ξk,3 , , 3 3 3
,
2ξk,1 2ξk,2 1 + 2ξk,3 , , 3 3 3
,
which are the barycentric coordinates of Zk and N1 , N2 , N3 with respect to Tk , one obtains the ordinates zk and n1 , n2 , n3 of Btk,i,d , i = 1, 2, 3, multiplied by 2. These ordinates are thus nonnegative, since the points Zk and N1 , N2 , N3 are all contained in the triangle Qtk,i by the assumptions of the lemma.
2 , π 3 ) to be the barycentric coordinates 1 , π To see that the ordinate o12 of Btk,i,d is nonnegative, consider ( τ1 , τ2 , τ3 ) and (π of N12 with respect to Tk and Qtk , respectively. By applying the transformation (16) and comparing the result to (8), it may be verified that the ordinate o12 of Btk,i,d is expressed as
o12 =
1 2
i + δi,1 ρ1 π
( − 1)d t σk,1,d (R1 − R12 )·wk,1 . 6ωk,1 (1 − ωk,1 )
If i = 1, o12 is obviously nonnegative. Otherwise notice that 2(ρ1 χ1 )−1 o12 is the barycentric coordinate of K12, d with respect to Ltk,1 corresponding to the point Q tk,1,d , which by assumptions proves that o12 is nonnegative. The values of o13 , o21 , o23 , o32 , and o31 of Btk,i,d may be expressed in an analogue way, and their nonnegativity follows by similar arguments.
Lemma 5. For a boundary edge ε Vi , Vj , let Le Q e,1 , Q e,2 be a line segment containing the points Lij and Lji . Then the basis functions Be,d , d = 1, 2, are nonnegative. Proof. All the Bézier ordinates of Be,d are equal to 0, except for rij , lij , lji . It may be easily verified that the transformation between the barycentric coordinates (τ 1 , τ 2 ) with respect to ε and the barycentric coordinates (π 1 , π 2 ) with respect to Le is given by
e e α,1 + η,1,1 π1 = e e π2 α,2 + η,2,1
e e α,1 + η,1,2 e e α,2 + η,2,2
τ1 , τ2
e e where η,d, ι = β,d (V ι − Ri j )·r i j , ι = i, j. By applying the transformation to the values
(λi j , λ ji ),
1 + 2λi j 2λ ji , 3 3
,
2λi j 1 + 2λ ji , 3 3
,
which are the barycentric coordinates of Rij , Lij , Lji with respect to ε , the values rij , lij , lji of Be,d for d = 1, 2 are obtained. This proves their nonnegativity and concludes the proof. With regard to Lemmas 3–5, the construction of nonnegative basis functions is purely geometric. Still, there remains a considerable freedom in the choice of the proposed objects. A choice of the triangle Qvi for a vertex V i ∈ V as prescribed in Lemma 3 is depicted in Fig. 6 (see also Fig. 3), and a choice of the triangle Qtk and line segments Ltk,i , i = 1, 2, 3, for a triangle
Tk V 1 , V 2 , V 3 ∈ according to Lemma 4 is sketched in Fig. 7 (see also Fig. 4). The triangles Qvi and Qtk are chosen as small as possible, and the line segments Ltk,i are as short as they can be. The first motivation for this comes from [19], where it is argued
that the stability of the considered basis for S21 (∗ ) depends on the area of the triangles used in the definition of basis functions. Since the proposed construction of the basis for S31 (∗ ) is much alike, it is reasonable to expect that similar conclusions may be derived here. The second reason will be indicated in the following subsection, where the domain objects will be extended to a control structure of a spline. The control structure will mimic the shape of the spline better if the domain objects are small. While the determination of the shortest line segments satisfying the given conditions is simple, it is much more difficult to find the smallest triangle containing a certain set of points. Since this problem is equivalent to finding the triangle with the smallest area that contains a convex hull of the point set, an algorithm described in [20,21] is proposed to efficiently solve the problem. 4.3. B-spline representation Suppose that a spline s ∈ S31 (∗ ) is determined with the interpolation values as in (5). To express it in the form (12), one may apply to s the interpolation functionals similarly as in the proof of Theorem 2. This leads to the systems of linear equations that uniquely determine the values of the coefficients. For a vertex Vi of , the corresponding coefficients of s may be expressed as v = f + f (X v − x ) + f (Y v − y ) ci,d i i,x i i,y i,d i i,d
(17a)
124
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
Fig. 6. A choice of a triangle Qv1 containing the points specified in Lemma 3.
Fig. 7. A choice of a triangle Qtk and line segments Ltk,i , i = 1, 2, 3, containing the points specified in Lemma 4.
for every d = 1, 2, 3, which can be easily verified by using (13). Considering (14) and the definition of the basis functions associated with a triangle Tk V 1 , V 2 , V 3 ∈ , one can similarly deduce that for every i = 1, 2, 3 and d = 1, 2,
ctk,i,d = ctk,i + ( − 1)d gwk,i (Q tk,i,2 − Q tk,i,1 )·wk,i ,
(17b)
where t t ctk,i = gk + gk,x (Xk,i − x∗k ) + gk,y (Yk,i − y∗k ).
Finally, the coefficients of s corresponding to a boundary edge ε Vi , Vj of , may be expressed as e c,d = h + hri j (Q e,d − Ri j )·ri j ,
d = 1, 2.
(17c)
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
125
Fig. 8. Modifying the control structure of a C 1 cubic spline on a Powell–Sabin triangulation. Only the control points with nonzero third component are depicted.
In what follows, let us introduce the control points C vi,d , C tk,i,d , C e,d to express the graph of s as |V | || |∂ E | 3 2 3 2 C vi,d Bvi,d + C tk,i,d Btk,i,d + C e,d Be,d . i=1 d=1
k=1 i=1 d=1
(18)
=1 d=1
By applying the interpolation data taken from the splines (x, y)→x and (x, y)→y to (17), one may easily verify that C vi,d and C e,d
v ) and (Q e , ce ), respectively, whereas C t are determined by the pairs (Q vi,d , ci,d ,d ,d k,i,d may be expressed as
C tk,i,d = C tk,i +
( − 1)d Q tk,i,2 − Q tk,i,1 , ctk,i,d − ctk,i
with C tk,i corresponding to (Q tk,i , ctk,i ). By Theorem 4, the graph lies in the affine hull of the control points. Naturally, the hull is convex if the basis functions are nonnegative. In this case, the position of the control points provides boundaries of the spline graph. In what follows, the control points of a spline are organized to a control structure that gives a good local control over the shape of the spline. Fig. 8 demonstrates how the control structure can be used in modelling. For every Vi denote by Civ C vi,1 , C vi,2 , C vi,3 the triangle whose projection to the domain is Qvi . Changing Civ so that the projection v , d = 1, 2, 3, and the graph of s is changed locally by Theorem 3. remains unchanged affects only the values of the coefficients ci,d v Since the triangle Ci is by (17a) tangent to the spline graph at Vi , it may be used as a geometric tool to handle the graph surface, and it thus deserves to be called a control triangle of s (Fig. 8a and b). t t C t Let Tk V 1 , V 2 , V 3 ∈ . Consider the triangle Ckt C tk,1 , C tk,2 , C tk,3 and the line segments Ck,i k,i,1 , C k,i,2 , i = 1, 2, 3, which are attached to the vertices of Ckt (Fig. 8c). This object is a generalization of the triangle Civ presented above. By moving the line t in a way that their projections remain unchanged, only the coefficients ct are modified, which by Theorem 3 segments Ck,i k,i,d
t is highly intuitive, and it may be approached has a local effect. Modifying the graph surface by changing the line segments Ck,i t with the point C t being clamped, one modifies in different ways. By fixing the triangle Ckt and changing the slope of a single Ck,i k,i
exclusively the coefficients ctk,i,d , d = 1, 2, and geometrically regulates the differential of s at Wk, i in the direction of wk, i (Fig. 8c
t , i = 1, 2, 3, and modifying C t on the other hand, one intuitively affects the graph surface in the and d). By fixing the slopes of Ck,i k
neighbourhood of s(Zk ) since Ckt is tangent to s at Zk (Fig. 8d and e). Raising and lowering of the fixed object gives the effect of extending and shrinking the spline locally. For a boundary edge ε Vi , Vj of , the control line Ce Q e,1 , Q e,2 provides a geometric tool similar to a control triangle. Through changing Ce above Le one may effectively modify the shape of the graph alongside the boundary of the domain (Fig. 8f).
126
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
The componentwise extension of spline functions to parametric spline surfaces is straightforward. The parametrization of the spline surface is given by (18) with the control points being considered as arbitrary points in the space. In this way, the control triangles and lines become independent of the domain and can be modified freely. The representation of parametric spline surfaces in terms of the introduced basis functions also enables a further generalization to rational spline surfaces, which indicates a rich potential of cubic splines on Powell–Sabin triangulations in computer aided geometric design. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]
M.J. Lai, L.L. Schumaker, Spline Functions on Triangulations, Cambridge University Press, Cambridge, 2007. M.J.D. Powell, M.A. Sabin, Piecewise quadratic approximations on triangles, ACM Trans. Math. Softw. 3 (1977) 316–325. P. Dierckx, On calculating normalized Powell–Sabin B-splines, Comput. Aided Geom. Des. 15 (1997) 61–78. H. Speleers, Construction of normalized B-splines for a family of smooth spline spaces over Powell–Sabin triangulations, Constr. Approx. 37 (2013) 41–72. O. Davydov, P. Sablonnière, C 2 piecewise cubic quasi-interpolants on a 6-direction mesh, J. Approx. Theory 162 (2010) 528–544. S. Remogna, Bivariate C 2 cubic spline quasi-interpolants on uniform Powell–Sabin triangulations, Adv. Comput. Math. 36 (2012) 39–65. S.-K. Chen, H.-W. Liu, A bivariate C 1 cubic super spline space on Powell–Sabin triangulation, Comput. Math. Appl. 56 (2008) 1395–1401. M. Lamnii, H. Mraoui, A. Tijini, A. Zidna, A normalized basis for C 1 cubic super spline space on Powell–Sabin triangulation, Math. Comput. Simul. 99 (2014) 108–124. A. Lamnii, M. Lamnii, H. Mraoui, Cubic spline quasi-interpolants on Powell–Sabin partitions, BIT Numer. Math. 54 (2014) 1099–1118. D. Sbibih, A. Serghini, A. Tijini, A. Zidna, Superconvergent C 1 cubic spline quasi-interpolants on Powell–Sabin partitions, BIT Numer. Math. (2014), doi:10.1007/s10543-014-0523-z. C. Manni, P. Sablonnière, Quadratic spline quasi-interpolants on Powell–Sabin partitions, Adv. Comput. Math. 26 (2007) 283–304. D. Sbibih, A. Seghini, A. Tijini, Polar forms and quadratic spline quasi-interpolants over Powell–Sabin partitions, Appl. Numer. Math. 59 (2009) 938–958. D. Sbibih, A. Serghini, A. Tijini, Superconvergent quadratic spline quasi-interpolants on Powell–Sabin partitions, Appl. Numer. Math. 87 (2015) 74–86. H. Speleers, P. Dierckx, S. Vandewalle, Numerical solution of partial differential equations with Powell–Sabin splines, J. Comput. Appl. Math. 189 (2006) 643–659. H. Speleers, P. Dierckx, S. Vandewalle, Local subdivision of Powell–Sabin splines, Comput. Aided Geom. Des. 23 (2006) 446–462. H. Speleers, P. Dierckx, S. Vandewalle, Quasi-hierarchical Powell-Sabin B-splines, Comput. Aided Geom. Des. 26 (2009) 174–191. H. Speleers, C. Manni, F. Pelosi, M.L. Sampoli, Isogeometric analysis with Powell–Sabin splines for advection-diffusion-reaction problems, Comput. Methods Appl. Mech. Eng. 221–222 (2012) 132–148. E. Vanraes, J. Windmolders, A. Bultheel, P. Dierckx, Automatic construction of control triangles for subdivided Powell–Sabin splines, Comput. Aided Geom. Des. 21 (2004) 671–682. J. Maes, E. Vanraes, P. Dierckx, A. Bultheel, On the stability of normalized Powell–Sabin B-splines, J. Comput. Appl. Math. 170 (2004) 181–196. V. Klee, M.C. Laskowski, Finding the smallest triangles containing a given convex polygon, J. Algorithms 6 (1985) 359–375. J. O’Rourke, A. Aggarwal, S. Maddila, M. Baldwin, An optimal algorithm for finding minimal enclosing triangles, J. Algorithms 7 (1986) 258–269.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
C 1 cubic splines on Powell–Sabin triangulations Jan Grošelj a,∗, Marjeta Krajnc a,b a b
FMF, University of Ljubljana, Jadranska 19, Ljubljana, Slovenia IMFM, Jadranska 19, Ljubljana, Slovenia
a r t i c l e
i n f o
Article history: Available online 27 July 2015 MSC: 41A05 41A15 65D07 65D17
a b s t r a c t A bivariate C 1 cubic spline space on a triangulation with Powell–Sabin refinement which extends the well-known C 1 quadratic spline space and has a nested structure is introduced. A construction of a locally supported basis that forms a partition of unity is presented based on choosing particular triangles and line segments in the domain. Further, it is shown how these objects can be determined in order to obtain nonnegative basis functions under a natural restriction on the Powell–Sabin refinement. Geometrically intuitive B-spline representation is proposed which makes these splines a useful tool for CAGD applications.
Keywords: Nested spline spaces Powell–Sabin triangulation Normalized B-splines Control structure
© 2015 Elsevier Inc. All rights reserved.
1. Introduction Spline functions are one of the most efficient and widely used tools in approximation theory and computer aided geometric design. Univariate spline spaces are very well examined, but a step to a bivariate case brings many difficulties and many problems are still unsolved. A bivariate spline is a function of a particular smoothness, composed of polynomial pieces defined on some partition that can be rectangular or triangular. Although tensor product splines on rectangular domains have recently attract much attention in isogeometric analysis, splines on triangulations are still widely used in function approximation, finite element method, for design purposes, etc. The bivariate spline spaces on triangulations have been intensively studied by many researchers (see [1] and the references therein), and it turns out that even the dimension determination is an extremely difficult problem since it heavily depends on the geometry of the triangulation. To be able to determine the dimension and to construct locally supported basis functions, one can use the so called macro elements [1]. Namely, every triangle of the initial triangulation is split into several smaller triangles, and a spline space on the refined triangulation is considered. The two refinements that are most commonly used are the Powell–Sabin and the Clough–Tocher split. A Powell–Sabin refinement splits every triangle of a triangulation into six smaller triangles. Splines on a Powell–Sabin refinement have first been studied in [2]. It has been shown therein that a C 1 quadratic spline on the refined triangulation is uniquely specified by the values and the first order derivatives at the vertices of the initial triangulation. The dual basis functions of the associated spline space have local supports, and it was proven in [3] that by slightly changing their interpolation values, nonnegative basis functions that form a partition of unity may be derived. A generalization of quadratic Powell–Sabin splines to higher degrees and higher orders of smoothness led to the so called Powell–Sabin super spline spaces, where the smoothness at some ∗
Corresponding author. E-mail address: [email protected] (J. Grošelj).
http://dx.doi.org/10.1016/j.amc.2015.07.013 0096-3003/© 2015 Elsevier Inc. All rights reserved.
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
115
vertices and across some edges is increased. In particular, C r super splines of degree 3r − 1 were presented in [4] together with a normalized representation. Cubic splines have also been considered in the literature, especially in terms of their application to the theory of quasiinterpolants. Based on the representation of C 2 cubic splines on the uniform Powell–Sabin refinement of a rectangular domain with pairs of translated multi-box splines, quasi-interpolants of differential and discrete type were derived in [5,6]. For arbitrary polygonal domains and general Powell–Sabin refinements, a C 1 cubic super spline space with C 2 smoothness conditions at the vertices of the initial triangulation was proposed in [7]. Its normalized representation with applications to quasi-interpolants was studied extensively in [8–10], thus generalizing the results for the quadratic case (see [11–13]). Increasing the smoothness at some points and across some edges of the triangulation makes the dimension counting problem easier and also simplifies the construction of the basis functions. But the super spline formulation has at least one important drawback, namely the super spline spaces cannot be nested. The sequence of spaces S0 , S1 , S2 , . . . defined on the triangulations 0 , 1 , 2 , . . . , such that n is a refinement of n−1 for each n, is a nested sequence of spaces if Sn−1 ⊂ Sn . Clearly this cannot hold for super spline spaces. The nested structure enables the use of (local) subdivision and allows the construction of hierarchical spline spaces, which are of a great importance in e.g. solving partial differential equations. For these applications with quadratic Powell–Sabin splines see [14–18]. In order to extend the number of Powell–Sabin spline spaces with the nested structure, a space of C 1 cubic splines on a Powell–Sabin refinement is introduced in this paper, which in comparison to the quadratic case allows more flexibility for the approximation and design purposes since more degrees of freedom are provided. Each element of the presented space is uniquely defined by prescribing the values and the first order derivatives at the vertices of the initial triangulation and at the triangle split points inside the initial triangles. Additionally, the directional derivatives at three chosen points in each initial triangle, and the values and the directional derivatives at the split point on each initial boundary edge are interpolated. Particularly chosen interpolation problems give locally supported basis functions which form a partition of unity. It is shown also that the basis functions are nonnegative under some natural restriction on the geometry of the Powell–Sabin refinement. In order to make the presented splines a useful and a user friendly design tool, a B-spline representation and a control structure that intuitively affects the shape of the spline are derived, too. The paper is organized as follows. In Section 2 some preliminaries about polynomials on triangles and splines on triangulations are presented. In Section 3 the C 1 cubic spline space on a Powell–Sabin triangulation is introduced, and the Bézier representation of its elements is derived in dependence of an interpolation problem. Section 4 gives a construction of a locally supported basis which forms a partition of unity. Geometric conditions on the construction, which ensure that the basis functions are nonnegative, are discussed. In the end, the B-spline representation is derived, and its application to geometric design is demonstrated. 2. Preliminaries 2.1. Bivariate polynomials on triangles Consider a triangle T V 1 , V 2 , V 3 with vertices Vi , i = 1, 2, 3. Each point P ∈ R2 can be uniquely expressed as an affine combination
P = τ1 V 1 + τ2 V 2 + τ3 V 3 ,
τ1 + τ2 + τ3 = 1.
The coefficients τ = (τ1 , τ2 , τ3 ) are called the barycentric coordinates of the point P with respect to the triangle T . If the point P lies inside the triangle, its barycentric coordinates are all positive. Let Pd denote the space of polynomials
p(x, y) =
ai j xi y j ,
ai j ∈ R,
0≤i+ j≤d
in two variables of total degree less than or equal to d ∈ N. Each bivariate polynomial p ∈ Pd can be expressed in the Bernstein– Bézier representation with respect to the triangle T as
p(P ) ≡ p(τ) =
bi Bdi (τ),
i = (i1 , i2 , i3 ),
|i| = i1 + i2 + i3 ,
(1)
|i|=d
where Bdi are the Bernstein basis polynomials of degree d, defined as
Bdi (τ) =
d! i1 ! i2 ! i3 !
τ1i1 τ2i2 τ3i3 , |i| = d,
and bi , |i| = d, are called the Bézier ordinates of p. The domain point Di that corresponds to the Bézier ordinate bi is determined
by the barycentric coordinates ( d1 , d2 , d3 ), and a pair (Di , bi ) is called a control point of p. The de Casteljau algorithm allows us to compute the values of p from (1) in a stable and efficient way. = ( 2 , η 3 ) be the barycentric coordinates of points Q = (xQ , yQ ) and Q 1 , η η = (η yQ ) with respect xQ , Let η = (η1 , η2 , η3 ) and with respect to T are equal to to T . The barycentric directions of the vector v := Q − Q i
i
i
1 , η2 − η 2 , η3 − η 3 ) ν = (ν1 , ν2 , ν3 ) = (η1 − η
116
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
and satisfy ν1 + ν2 + ν3 = 0. The directional derivative of p at the point P in the direction of v is defined as
Dv p(P ) :=
d p(P + t v) = Dx p(P )(xQ − xQ ) + Dy p(P )(yQ − yQ ) dt
and equals
Dv p(P ) = d
|i|=d−1
(ν1 bi+(1,0,0) + ν2 bi+(0,1,0) + ν3 bi+(0,0,1) )Bd−1 (τ). i
(2)
(3)
2.2. Spline spaces on triangulations Let ⊂ R2 be a simply connected domain with a polygonal boundary, and let = {Tk }k be its arbitrary regular triangulation. The sets of vertices and edges of are denoted by V = {V i }i and E = {ε } , respectively. The sets are assumed to be finite, and their elements are indexed with the indices from 1 to the cardinality of a set, which is denoted by | · |. The edge ε is sometimes written as ε Vi , Vj to indicate that its endpoints are the vertices Vi and Vj . The edge lying on the boundary of the domain is a boundary edge. The set of all boundary edges is denoted by ∂ E. The subset of all triangles in sharing the vertex Vi is called the molecule of Vi and is denoted by M (Vi ). A space of polynomial splines of degree d ∈ N and smoothness r ∈ N0 on a triangulation is defined as
Sdr () := {s ∈ C r (); s|T ∈ Pd , T ∈ }. It consists of splines composed of polynomials of degree d, defined on triangles of , that join together with C r smoothness. By using (1) each of the polynomial segments can be locally represented in the Bernstein–Bézier form and is completely defined by prescribing its control points. The union of the Bézier representations (1) of all the polynomial segments gives the Bézier spline representation. Smoothness constraints imply certain relations between the control points of the neighbouring polynomials. In particular, let p and p be two polynomials of total degree d written in the form (1) on T V 1 , V 2 , V 3 and T V 4 , V 3 , V 2 with Bézier ordinates bi and b i , respectively. Suppose that V4 lies outside of T and has barycentric coordinates σ = (σ1 , σ2 , σ3 ) with respect to T . Then p and p are C 1 continuous across the common edge V2 , V3 if and only if
b 0,i2 ,i3 = b0,i3 ,i2 , b
1,i2 ,i3
i2 + i3 = d,
(4a)
= σ1 b1,i3 ,i2 + σ2 b0,i3 +1,i2 + σ3 b0,i3 ,i2 +1 ,
i2 + i3 = d − 1.
(4b)
The proof and the higher order smoothness conditions can be found in [1]. 2.3. Powell–Sabin refinement A Powell–Sabin refinement ∗ of the initial triangulation , introduced in [2], splits every triangle into six smaller triangles in the following way. First, for every triangle Tk ∈ , k = 1, 2, . . . , n , a point Zk in the interior of Tk is chosen and is connected to the vertices of the triangle. The newly inserted points can be chosen as the incenters of the corresponding triangles or, more generally, can be chosen in such a way that for every two triangles Tk and Tk with a common edge Vi , Vj the intersection point Rij of the line segments Z k , Z k and Vi , Vj lies in the interior of Vi , Vj . The point Rij is then connected to Zk and Z k . Additional points are chosen on the boundary edges of the triangulation and are connected to the interior split points of the associated triangles. An example of a Powell–Sabin refinement of two neighbouring triangles is shown in Fig. 1. The following notation will further be needed. For every triangle split point Z k ∈ Tk V 1 , V 2 , V 3 , let ξ k, 1 , ξ k, 2 , ξ k, 3 be such that
Z k = ξk,1V 1 + ξk,2V 2 + ξk,3V 3 ,
ξk,1 + ξk,2 + ξk,3 = 1.
For every edge split point Rij ∈ Vi , Vj , let λij and λji be such that
Ri j = λi jV i + λ jiV j ,
λi j + λ ji = 1.
Furthermore, the cartesian coordinates of the points Vi , Rij , and Zk are denoted by (xi , yi ), (x∗i j , y∗i j ), and (x∗k , y∗k ), respectively. Note that the points Rij and Rji coincide. 3. A space of C 1 cubic splines Let be an arbitrary triangulation and ∗ its Powell–Sabin refinement. Consider the spline space S31 (∗ ) of C 1 cubic splines on ∗ . The aim of this section is to provide an interpolation problem that uniquely describes an element of S31 (∗ ). Suppose that s ∈ S31 (∗ ). For a vertex V i ∈ V, let fi , fi, x , and fi, y be the interpolation values such that
s(V i ) = fi ,
Dx s(V i ) = fi,x ,
Dy s(V i ) = fi,y .
(5a)
For a triangle Tk V 1 , V 2 , V 3 ∈ , let gk , gk, x , and gk, y be the interpolation values such that
s(Z k ) = gk ,
Dx s(Z k ) = gk,x ,
Dy s(Z k ) = gk,y .
(5b)
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
117
Fig. 1. A Powell–Sabin refinement of two neighbouring triangles.
Fig. 2. Domain points of a cubic spline on a macro-triangle Tk V 1 , V 2 , V 3 .
Moreover, for a triangle Tk ∈ and each of its vertices Vi , choose a point Wk, i lying on the edge Vi , Zk , i.e.,
W k,i = ωk,iV i + (1 − ωk,i )Z k ,
0 < ωk,i < 1.
Additionally, choose any unit vector wk, i nonparallel to Z k − V i , and let gwk,i be the interpolation values such that
Dwk,i s(W k,i ) = gwk,i .
(5c)
Finally, for a boundary edge V i , V j ∈ ∂ E, let rij be a unit vector parallel to V j − V i , and let hij and hri j be the interpolation values such that
s(Ri j ) = hi j ,
Dri j s(Ri j ) = hri j .
(5d) S31 (∗ ).
The conditions in (5) specify an interpolation problem for In what follows, it is proven that the interpolation problem admits a unique solution. Let Tk V 1 , V 2 , V 3 be an arbitrary triangle. It will be assumed hereafter that the domain points of a spline s ∈ S31 ({Tk }∗ ) on Tk are denoted as in Fig. 2. The points P i j = P ji and T i j = T ji , i, j ∈ {1, 2, 3}, will be identified to simplify the notation. Moreover, the Bézier ordinate of s associated with a domain point will be denoted by the lower case letter. For example, the Bézier ordinate associated with the domain point O12 will be denoted by o12 .
118
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
Lemma 1. Let Tk V 1 , V 2 , V 3 be an arbitrary triangle. There exists a unique spline s ∈ S31 ({Tk }∗ ) satisfying the conditions in (5). Proof. Suppose that an arbitrary set of interpolation values prescribing the interpolation conditions in (5) is given. In what follows, the Bézier ordinates of the spline s from S31 ({Tk }∗ ) that satisfies these conditions are derived, which proves both, the existence and the uniqueness of s. By (2) and (3) the conditions (5a) imply vi = fi and
1 ( fi,x (x∗k − xi ) + fi,y (y∗k − yi )), 3 1 ui j = fi + ( fi,x (x∗i j − xi ) + fi,y (y∗i j − yi )) 3 mi = fi +
(6)
for every {i, j}⊆{1, 2, 3}. Analogously by (5b), zk = gk and
1 (g (xi − x∗k ) + gk,y (yi − y∗k )), 3 k,x 1 ti j = gk + (gk,x (x∗i j − x∗k ) + gk,y (y∗i j − y∗k )). 3 ni = gk +
(7)
Furthermore, let (wij, 1 , wij, 2 , wij, 3 ) denote the barycentric directions of wk, i with respect to the triangle Vi , Zk , Rij . By (3) and (5c),
wi j,2 wi j,2 wi j,1 wi j,1 ωk,i oi j = − vi + mi + ui j − mi + ni wi j,3 wi j,3 wi j,3 2(1 − ωk,i ) wi j,3 wi j,2 1 − ωk,i wi j,1 1 1 − ni + z + ti j + gw . 2ωk,i wi j,3 wi j,3 k 6ωk,i (1 − ωk,i ) wi j,3 k,i
(8)
By the C 1 smoothness conditions (4b) across the edges Rij , Zk , it follows:
pi j = λi j oi j + λ ji o ji .
(9)
Finally, by (2) and (3) the conditions (5d) imply ri j = hi j and
li j = hi j +
1 hr (V i − Ri j )·ri j , 3 ij
(10)
where · denotes the standard scalar product in R2 . The conditions in (4a) are clearly satisfied everywhere by the notation. The conditions (4b) across the edges Rij , Zk are satisfied by (7), (9), (10), and the fulfillment of the same conditions across Vi , Zk can be easily verified by following [1, p. 35]. Since all the Bézier ordinates of s are now determined and all smoothness conditions are satisfied, the proof is completed. Lemma 2. Let Tk V 1 , V 2 , V 3 and Tk V 4 , V 3 , V 2 be two triangles with the common edge V2 , V3 . There exists a unique spline s ∈ S31 ({Tk , Tk }∗ ) satisfying the conditions in (5). Proof. All the Bézier ordinates of the spline s may be expressed in the same way as in Lemma 1, except for the values r23 , l23 , and l32 . It needs to be checked that they can be uniquely determined, and that s satisfies the smoothness conditions along and across the edge V2 , V3 .
Suppose that R23 = μ23,k Z k + μ23,k Z k , μ23,k + μ23,k = 1. Using (6), it may be easily verified that ui j = μ23,k mki + μ32,k mki , 1 {i, j} = {2, 3}, where the upper index denotes the triangle that the ordinate belongs to. The other three C smoothness conditions (4b) across V2 , V3 imply
r23 = μ23,k pk23 + μ23,k pk23 ,
l23 = μ23,k ok23 + μ23,k ok23 , l32 = μ
k 23,k o32
+μ
(11)
k
23,k o32 .
It follows from (9) and (11) that
r23 = μ23,k pk23 + μ23,k pk23
= μ23,k (λ23 ok23 + λ32 ok32 ) + μ23,k (λ23 ok23 + λ32 ok32 )
= λ23 (μ23,k ok23 + μ23,k ok23 ) + λ32 (μ23,k ok32 + μ23,k ok32 ) = λ23 l23 + λ32 l32 , which means that the C 1 smoothness condition at R23 in the direction of V 3 − V 2 is satisfied, and hence the lemma is proven. Theorem 1. Let be an arbitrary triangulation. There exists a unique spline s ∈ S31 (∗ ) satisfying the conditions in (5). Proof. The proof follows directly from Lemmas 1 and 2.
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
119
Fig. 3. A choice of a triangle Qvi Q vi,1 , Q vi,2 , Q vi,3 (left) and the associated basis function Bvi,1 (right).
According to Theorem 1, the dimension of S31 (∗ ) can be determined by a simple count of the interpolation conditions in (5). Namely, (5a) gives three conditions per vertex, each of (5b) and (5c) three conditions per triangle, and (5d) two conditions per boundary edge. Corollary 1. The dimension of S31 (∗ ) is equal to 3|V | + 6|| + 2|∂ E |. In order to show that the considered spline space has a nested structure, a suitable refinement R of has to be found such that a Powell–Sabin refinement ∗R of R may be constructed from the existing Powell–Sabin refinement ∗ . When this is the case, it follows trivially that S31 (∗ ) ⊆ S31 (∗R ) since the restriction of s ∈ S31 (∗ ) to each T ∗ ∈ ∗ is a cubic polynomial, which can be subdivided according to the refinement of T ∗ implied by ∗R , and all smoothness conditions of S31 (∗R ) are fulfilled for the spline obtained in this way. One obvious choice for R is ∗ , and a Powell–Sabin refinement of R can be chosen freely. However, with such a refinement the number of triangles grows very fast, which may not be desirable. A common refinement of a triangulation is a dyadic split that divides every triangle into four smaller triangles and is obtained by adding an additional point on every edge of the triangulation. For a triangle Tk V 1 , V 2 , V 3 ∈ , the candidates for additional points would be R12 , this approach yields an incompatible refinement when the triangle R12 , R23 , R31 does not R23 , R31 , but as illustrated in [18], √ contain Zk . As a better√practice the 3-subdivision rule is proposed therein, which results in a triadic split if one applies it twice. The advantage of the 3-rule is also that it allows a local refinement as proposed in [15] for the quadratic case. 4. A normalized basis In this section a construction of a basis for a space S31 (∗ ) is presented. It introduces basis functions of three kinds, each of them is associated with a vertex, a triangle or a boundary edge of the initial triangulation. More precisely, the aim is to represent an arbitrary spline s ∈ S31 (∗ ) as a linear combination
s=
|V | 3
v Bv + ci,d i,d
i=1 d=1
|| 3 2
ctk,i,d Btk,i,d +
k=1 i=1 d=1
|∂ E | 2 =1 d=1
e c,d Be,d ,
(12)
where the basis functions Bvi,d are associated with the vertex V i ∈ V, the basis functions Btk,i,d are associated with the triangle Tk V 1 , V 2 , V 3 ∈ , and the basis functions Be,d are associated with the boundary edge ε V i , V j ∈ ∂ E. Beside the construction it will be shown that it is possible to define the basis functions in a way that they have local supports and form a partition of unity. It will also be proven that they are nonnegative under a certain assumption on the geometry of the Powell–Sabin refinement. 4.1. Basis functions Inspired by the ideas of a geometric construction of a normalized basis for S21 (∗ ) presented in [3], a basis for S31 (∗ ) may be constructed as follows. For a demonstration, consider Figs. 3–5. v , α v , α v ) be the barycentric coordinates of V 1. For every V i ∈ V, choose a triangle Qvi Q vi,1 , Q vi,2 , Q vi,3 in the plane. Let (αi,1 i i,2 i,3 v , β v , β v ) and (γ v , γ v , γ v ) be the barycentric directions of (1, 0) and (0, 1) with respect to with respect to Qvi , and let (βi,1 i,2 i,3 i,1 i,2 i,3
Qvi . Then Bvi,d , d = 1, 2, 3, are the splines from S31 (∗ ) uniquely defined by setting the interpolation values (5a) to v , Bvi,d (V i ) = αi,d
v , Dx Bvi,d (V i ) = βi,d
v , Dy Bvi,d (V i ) = γi,d
and all the other interpolation values in (5) to zero.
120
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
Fig. 4. A choice of a triangle Qtk Q tk,1 , Q tk,2 , Q tk,3 and line segments Ltk,i Q tk,i,1 , Q tk,i,2 , i = 1, 2, 3 (left), and the associated basis function Btk,3,2 (right).
Fig. 5. A choice of a line segment Le Q e,1 , Q e,2 (left) and the associated basis function Be,2 (right).
t , α t , α t ) be the barycentric co2. For every Tk V 1 , V 2 , V 3 ∈ , choose a triangle Qtk Q tk,1 , Q tk,2 , Q tk,3 in the plane. Let (αk,1 k,2 k,3
t , β t , β t ) and (γ t , γ t , γ t ) be the barycentric directions of (1, 0) and ordinates of Zk with respect to Qtk , and let (βk,1 k,2 k,3 k,1 k,2 k,3
(0, 1) with respect to Qtk . Additionally, choose line segments Ltk,i Q tk,i,1 , Q tk,i,2 , i = 1, 2, 3, on the lines through Qk, i with the
t t direction wk, i , respectively. Let (σk,i,1 , σk,i,2 ) be the barycentric directions of wk, i with respect to Ltk,i . Then Btk,i,d , i = 1, 2, 3,
d = 1, 2, are the splines from S31 (∗ ) uniquely defined by setting the interpolation values (5b) to
Btk,i,d (Z k ) =
1 2
t t t αk,i , Dx Btk,i,d (Z k ) = 12 βk,i , Dy Btk,i,d (Z k ) = 12 γk,i ,
and (5c) to
Dwk,i Btk,i,d (W k, ) =
t δ σk,i,d , = 1, 2, 3,
1 2 i,
with all the other interpolation values in (5) being equal to zero. Note that δ i, denotes the Kronecker delta. e , α e ) be the 3. For every ε V i , V j ∈ ∂ E, choose a line segment Le Q e,1 , Q e,2 on the line containing the edge ε . Let (α,1 ,2 e , β e ) be the barycentric directions of r with respect to Le . barycentric coordinates of Rij with respect to Le , and let (β,1 ij ,2 Then Be,d , d = 1, 2, are the splines from S31 (∗ ) uniquely defined by setting the interpolation values (5d) to e Be,d (Ri j ) = α,d ,
e Dri j Be,d (Ri j ) = β,d ,
and all the other interpolation values in (5) to zero.
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
121
It needs to be proven that the presented splines form a basis for S31 (∗ ) which admits the properties of a local support and a v , β v , γ v , d = 1, 2, 3, satisfy the equality partition of unity. Note first that for every i = 1, 2, . . . , |V | the parameters αi,d i,d i,d
⎡
1 v ⎣Xi,1 v Yi,1
1 v Xi,2 v Yi,2
⎤⎡
v αi,1 1 v ⎦⎣α v Xi,3 i,2 v v Yi,3 αi,3
v βi,1 v βi,2 v βi,3
⎤
v γi,1 1 v ⎦= x γi,2 i v yi γi,3
0 1 0
0 0 , 1
(13)
v , Y v ) denote the cartesian coordinates of Q v , d = 1, 2, 3. Similarly, for every k = 1, 2, . . . , || the parameters α t , where (Xi,d i,d i,d k,i
t , γ t , i = 1, 2, 3, satisfy βk,i k,i
⎡
1 t ⎣Xk,1 t Yk,1
1 t Xk,2 t Yk,2
⎤⎡
t βk,1 t βk,2 t βk,3
t 1 αk,1 t ⎦⎣ t Xk,3 αk,2 t t Yk,3 αk,3
⎤ ⎡ t γk,1 1 t ⎦ γk,2 = ⎣x∗k t γk,3 y∗k
0 1 0
⎤
0 0⎦ , 1
(14)
t , Y t ) are the cartesian coordinates of Q t , i = 1, 2, 3. where (Xk,i k,i k,i
Theorem 2. The splines Bvi,d , Btk,i,d , Be,d that appear in (12) form a basis of S31 (∗ ). Proof. Since the number of the splines is equal to the dimension of the spline space S31 (∗ ), it is sufficient to prove that they are linearly independent. Suppose that the expression in (12) is equal to 0. The evaluation of s, Dx s, and Dy s at V i ∈ V results in the system of equations 3
v αv = ci,d i,d
d=1
3
v βv = ci,d i,d
3
d=1
v γ v = 0, ci,d i,d
d=1
v = 0, d = 1, 2, 3. Let T V , V , V ∈ . Similarly as above, the evaluation which has a unique solution by (13). This means that ci,d 1 2 3 k of s, Dx s, and Dy s at Zk gives the system of equations 3
ctk,i,1 + ctk,i,2
i=1
t = αk,i
3
ctk,i,1 + ctk,i,2
i=1
t = βk,i
3
ctk,i,1 + ctk,i,2
t = 0. γk,i
i=1
The evaluation of Dwk,i s, i = 1, 2, 3, at Wk, i extends the system with the equations t t ctk,i,1 σk,i,1 + ctk,i,2 σk,i,2 = 0,
i = 1, 2, 3.
Using (14), one may easily verify that the system of six equations has a unique solution. This proves that ctk,i,d = 0, i = 1, 2, 3, d = 1, 2. Finally, let ε V i , V j ∈ ∂ E. Evaluation of s and Dri j s at Rij results in 2 d=1
e e c,d α,d =
2
e e c,d β,d = 0,
d=1
e = 0, d = 1, 2. from where it follows that c,d
Theorem 3. The splines Bvi,d , Btk,i,d , Be,d that appear in (12) have local supports. Proof. Consider the Bézier ordinates of Bvi,d on a triangle of outside M (Vi ). Since the spline on this triangle is determined with zero interpolation values, all ordinates are zero as can be verified by following the proofs of Lemmas 1 and 2. This shows that the support of Bvi,d is contained in M (Vi ). A similar deduction shows that on a triangle of that does not share any edge with Tk ∈ all the Bézier ordinates of Btk,i,d are zero, and the support of this basis function is contained in the area spanned by Tk and its neighbours. The only nonzero Bézier ordinates of Be,d correspond to the domain points that are located on the edge ε V i , V j ∈ ∂ E, and so the support of Be,d is contained in a single triangle of . Theorem 4. The splines Bvi,d , Btk,i,d , Be,d that appear in (12) form a partition of unity. Proof. It is easy to verify that the sum of the basis functions satisfies the interpolation conditions for the constant 1. By Theorem 1, the proof follows. 4.2. Nonnegativity of the basis functions v , ct e , In order to view the spline s ∈ S31 (∗ ) expressed in the form (12) as a convex combination of the coefficients ci,d , c,d k,i,d
the basis functions Bvi,d , Btk,i,d , Be,d have to be nonnegative. This property imposes additional constraints on the parameters that determine the basis functions. Unfortunately, these constraints cannot be satisfied in general, and some additional assumptions on the geometry of the Powell–Sabin refinement are needed. In what follows, it will be assumed that for every triangle
122
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
Tk V 1 , V 2 , V 3 ∈ the interior split point Zk lies in the interior of the triangle R12 , R23 , R31 . This requirement is sufficient but not necessary for the construction of nonnegative basis functions. Let Tk V 1 , V 2 , V 3 be a triangle in . Suppose hereafter that the vector wk, i , i = 1, 2, 3, specified in the interpolation problem (5), is a unit vector parallel to Ri, j1 − Ri, j2 , { j1 , j2 } = {1, 2, 3}\{i}. Denote by
Ri = (1 − ρi )V i + ρi Z k the intersection point of Vi , Zk and Ri, j1 , Ri, j2 . Note that by the assumption, 0 < ρ i < 1. The construction of nonnegative basis functions is based on imposing geometric conditions, expressed as the requirements that the triangles Qvi and Qtk and the line segments Ltk,i and Le contain certain sets of points. Lemma 3. Let V1 be an arbitrary vertex of . If the triangle Qv1 contains the points V1 , U12 , U13 , M1 ,
M 12 = M 1 + M 13 = M 1 +
ωk,1
6(1 − ωk,1 )(1 − ρ1 )
ωk,1
6(1 − ωk,1 )(1 − ρ1 )
(R1 − R12 ), (R1 − R13 )
for every Tk V 1 , V 2 , V 3 ∈ M (V 1 ), then the basis functions Bv1,d , d = 1, 2, 3, are nonnegative. Proof. By the nonnegativity of the Bernstein basis polynomials, it is sufficient to check that the Bézier ordinates of Bv1,d , d = 1, 2, 3, are nonnegative. By Theorem 3, only the Bézier ordinates corresponding to the domain points in M (V1 ) are nonzero. Furthermore, as follows from the proofs of Lemmas 1 and 2, the nonnegativity of the ordinates needs to be verified only for v1 , u12 , u13 , m1 , o12 , and o13 associated with an arbitrary triangle Tk V 1 , V 2 , V 3 in M (V1 ). To show this, recall the property that the barycentric coordinates of a point with respect to a triangle are nonnegative if and only if the point lies inside the triangle. Consider the transformation
⎡ 1 π1 v π2 = ⎣X1,1 v π3 Y1,1
⎤−1 ⎡
1 v X1,2 v Y1,2
1 v ⎦ X1,3 v Y1,3
1 ⎣ x1 y1
1 x2 y2
⎤ τ1 τ2 τ3
1 x3 ⎦ y3
between the barycentric coordinates (τ 1 , τ 2 , τ 3 ) with respect to Tk and the barycentric coordinates (π 1 , π 2 , π 3 ) with respect to Qv1 . By (13), it may be rewritten as
⎡α v π1 1,1 v π2 = ⎣α1,2 v π3 α1,3
v + ηv α1,1 1,1,2 v + ηv α1,2 1,2,2 v + ηv α1,3 1,3,2
⎤ v + ηv α1,1 τ1 1,1,3 v + ηv ⎦ τ2 , α1,2 1,2,3 v + ηv τ3 α1,3 1,3,3
(15)
v v (x − x ) + γ v (y − y ), d = 1, 2, 3, j = 2, 3. The barycentric coordinates of V , U , U , and M with respect = β1,d where η1,d, 1 1 1 12 13 1 j j j 1,d to Tk are given by
(1, 0, 0),
2 + λ12 λ21 , ,0 , 3 3
2 + λ13 λ31 , 0, 3 3
,
2 + ξk,1 ξk,2 ξk,3 , , 3 3 3
.
It is easy to see that by applying the transformation (15) one obtains the ordinates v1 , u12 , u13 , and m1 of Bv1,d , d = 1, 2, 3, corresponding to Tk . It remains to check that the values of o12 and o13 are nonnegative. Using the same approach as above, one can see that these values can be expressed as the barycentric coordinates of M12 and M13 with respect to Qv1 , multiplied by 1 − ρ1 . By the assumption, this factor is positive, and the result follows. Lemma 4. For a triangle Tk V 1 , V 2 , V 3 ∈ , choose Qtk Q tk,1 , Q tk,2 , Q tk,3 to be a triangle containing the points
Ni j = Ni +
1 − ωk,i (Ri − Ri j ), 6ωk,i ρi
i ∈ {1, 2, 3}, j ∈ {1, 2, 3}\{i},
such that for each i = 1, 2, 3 none of the points N i, j1 and N i, j2 lies on the line segment Q tk, j , Q tk, j , { j1 , j2 } = {1, 2, 3}\{i}. Suppose 1
2
that χ ij is the barycentric coordinate of Nij with respect to the triangle Qtk corresponding to Q tk,i , and let χi = χi, j1 + χi, j2 , { j1 , j2 } =
{1, 2, 3}\{i}. For each i = 1, 2, 3, choose Ltk,i Q tk,i,1 , Q tk,i,2 to be a line segment containing the points K i j,d = Q tk,i +
( − 1)d (Ri − Ri j ), 6ωk,i (1 − ωk,i )ρi χi χi, j
j ∈ {1, 2, 3}\{i}, d ∈ {1, 2},
χi, j
with the property that Q tk,i = χ 1 Q tk,i,1 + χ 2 Q tk,i,2 . Then the basis functions Btk,i,d , i = 1, 2, 3, d = 1, 2, are nonnegative. i i Proof. By a similar deduction as in Lemma 3, one finds out that the proof follows if and only if the ordinates zk , n1 , n2 , n3 , o12 , o13 , o21 , o23 , o32 , o31 of Btk,i,d corresponding to the domain points in Tk V 1 , V 2 , V 3 are nonnegative. The transformation between
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
123
the barycentric coordinates (τ 1 , τ 2 , τ 3 ) with respect to Tk and the barycentric coordinates (π 1 , π 2 , π 3 ) with respect to Qtk is now given by
⎡αt + ηt π1 k,1 k,1,1 t t + ηk,2,1 π2 = ⎣αk,2 t t π3 αk,3 + ηk,3,1
t t αk,1 + ηk,1,2 t t αk,2 + ηk,2,2 t t αk,3 + ηk,3,2
⎤ t t αk,1 + ηk,1,3 τ1 t t ⎦ τ2 , αk,2 + ηk,2,3 t t τ3 αk,3 + ηk,3,3
(16)
t t (x − x∗ ) + γ t (y − y∗ ), i, j = 1, 2, 3. By applying the transformation (16) to (ξ = βk,i where ηk,i, j j k, 1 , ξ k, 2 , ξ k, 3 ) and k k j k,i
1 + 2ξk,1 2ξk,2 2ξk,3 , , 3 3 3
,
2ξk,1 1 + 2ξk,2 2ξk,3 , , 3 3 3
,
2ξk,1 2ξk,2 1 + 2ξk,3 , , 3 3 3
,
which are the barycentric coordinates of Zk and N1 , N2 , N3 with respect to Tk , one obtains the ordinates zk and n1 , n2 , n3 of Btk,i,d , i = 1, 2, 3, multiplied by 2. These ordinates are thus nonnegative, since the points Zk and N1 , N2 , N3 are all contained in the triangle Qtk,i by the assumptions of the lemma.
2 , π 3 ) to be the barycentric coordinates 1 , π To see that the ordinate o12 of Btk,i,d is nonnegative, consider ( τ1 , τ2 , τ3 ) and (π of N12 with respect to Tk and Qtk , respectively. By applying the transformation (16) and comparing the result to (8), it may be verified that the ordinate o12 of Btk,i,d is expressed as
o12 =
1 2
i + δi,1 ρ1 π
( − 1)d t σk,1,d (R1 − R12 )·wk,1 . 6ωk,1 (1 − ωk,1 )
If i = 1, o12 is obviously nonnegative. Otherwise notice that 2(ρ1 χ1 )−1 o12 is the barycentric coordinate of K12, d with respect to Ltk,1 corresponding to the point Q tk,1,d , which by assumptions proves that o12 is nonnegative. The values of o13 , o21 , o23 , o32 , and o31 of Btk,i,d may be expressed in an analogue way, and their nonnegativity follows by similar arguments.
Lemma 5. For a boundary edge ε Vi , Vj , let Le Q e,1 , Q e,2 be a line segment containing the points Lij and Lji . Then the basis functions Be,d , d = 1, 2, are nonnegative. Proof. All the Bézier ordinates of Be,d are equal to 0, except for rij , lij , lji . It may be easily verified that the transformation between the barycentric coordinates (τ 1 , τ 2 ) with respect to ε and the barycentric coordinates (π 1 , π 2 ) with respect to Le is given by
e e α,1 + η,1,1 π1 = e e π2 α,2 + η,2,1
e e α,1 + η,1,2 e e α,2 + η,2,2
τ1 , τ2
e e where η,d, ι = β,d (V ι − Ri j )·r i j , ι = i, j. By applying the transformation to the values
(λi j , λ ji ),
1 + 2λi j 2λ ji , 3 3
,
2λi j 1 + 2λ ji , 3 3
,
which are the barycentric coordinates of Rij , Lij , Lji with respect to ε , the values rij , lij , lji of Be,d for d = 1, 2 are obtained. This proves their nonnegativity and concludes the proof. With regard to Lemmas 3–5, the construction of nonnegative basis functions is purely geometric. Still, there remains a considerable freedom in the choice of the proposed objects. A choice of the triangle Qvi for a vertex V i ∈ V as prescribed in Lemma 3 is depicted in Fig. 6 (see also Fig. 3), and a choice of the triangle Qtk and line segments Ltk,i , i = 1, 2, 3, for a triangle
Tk V 1 , V 2 , V 3 ∈ according to Lemma 4 is sketched in Fig. 7 (see also Fig. 4). The triangles Qvi and Qtk are chosen as small as possible, and the line segments Ltk,i are as short as they can be. The first motivation for this comes from [19], where it is argued
that the stability of the considered basis for S21 (∗ ) depends on the area of the triangles used in the definition of basis functions. Since the proposed construction of the basis for S31 (∗ ) is much alike, it is reasonable to expect that similar conclusions may be derived here. The second reason will be indicated in the following subsection, where the domain objects will be extended to a control structure of a spline. The control structure will mimic the shape of the spline better if the domain objects are small. While the determination of the shortest line segments satisfying the given conditions is simple, it is much more difficult to find the smallest triangle containing a certain set of points. Since this problem is equivalent to finding the triangle with the smallest area that contains a convex hull of the point set, an algorithm described in [20,21] is proposed to efficiently solve the problem. 4.3. B-spline representation Suppose that a spline s ∈ S31 (∗ ) is determined with the interpolation values as in (5). To express it in the form (12), one may apply to s the interpolation functionals similarly as in the proof of Theorem 2. This leads to the systems of linear equations that uniquely determine the values of the coefficients. For a vertex Vi of , the corresponding coefficients of s may be expressed as v = f + f (X v − x ) + f (Y v − y ) ci,d i i,x i i,y i,d i i,d
(17a)
124
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
Fig. 6. A choice of a triangle Qv1 containing the points specified in Lemma 3.
Fig. 7. A choice of a triangle Qtk and line segments Ltk,i , i = 1, 2, 3, containing the points specified in Lemma 4.
for every d = 1, 2, 3, which can be easily verified by using (13). Considering (14) and the definition of the basis functions associated with a triangle Tk V 1 , V 2 , V 3 ∈ , one can similarly deduce that for every i = 1, 2, 3 and d = 1, 2,
ctk,i,d = ctk,i + ( − 1)d gwk,i (Q tk,i,2 − Q tk,i,1 )·wk,i ,
(17b)
where t t ctk,i = gk + gk,x (Xk,i − x∗k ) + gk,y (Yk,i − y∗k ).
Finally, the coefficients of s corresponding to a boundary edge ε Vi , Vj of , may be expressed as e c,d = h + hri j (Q e,d − Ri j )·ri j ,
d = 1, 2.
(17c)
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
125
Fig. 8. Modifying the control structure of a C 1 cubic spline on a Powell–Sabin triangulation. Only the control points with nonzero third component are depicted.
In what follows, let us introduce the control points C vi,d , C tk,i,d , C e,d to express the graph of s as |V | || |∂ E | 3 2 3 2 C vi,d Bvi,d + C tk,i,d Btk,i,d + C e,d Be,d . i=1 d=1
k=1 i=1 d=1
(18)
=1 d=1
By applying the interpolation data taken from the splines (x, y)→x and (x, y)→y to (17), one may easily verify that C vi,d and C e,d
v ) and (Q e , ce ), respectively, whereas C t are determined by the pairs (Q vi,d , ci,d ,d ,d k,i,d may be expressed as
C tk,i,d = C tk,i +
( − 1)d Q tk,i,2 − Q tk,i,1 , ctk,i,d − ctk,i
with C tk,i corresponding to (Q tk,i , ctk,i ). By Theorem 4, the graph lies in the affine hull of the control points. Naturally, the hull is convex if the basis functions are nonnegative. In this case, the position of the control points provides boundaries of the spline graph. In what follows, the control points of a spline are organized to a control structure that gives a good local control over the shape of the spline. Fig. 8 demonstrates how the control structure can be used in modelling. For every Vi denote by Civ C vi,1 , C vi,2 , C vi,3 the triangle whose projection to the domain is Qvi . Changing Civ so that the projection v , d = 1, 2, 3, and the graph of s is changed locally by Theorem 3. remains unchanged affects only the values of the coefficients ci,d v Since the triangle Ci is by (17a) tangent to the spline graph at Vi , it may be used as a geometric tool to handle the graph surface, and it thus deserves to be called a control triangle of s (Fig. 8a and b). t t C t Let Tk V 1 , V 2 , V 3 ∈ . Consider the triangle Ckt C tk,1 , C tk,2 , C tk,3 and the line segments Ck,i k,i,1 , C k,i,2 , i = 1, 2, 3, which are attached to the vertices of Ckt (Fig. 8c). This object is a generalization of the triangle Civ presented above. By moving the line t in a way that their projections remain unchanged, only the coefficients ct are modified, which by Theorem 3 segments Ck,i k,i,d
t is highly intuitive, and it may be approached has a local effect. Modifying the graph surface by changing the line segments Ck,i t with the point C t being clamped, one modifies in different ways. By fixing the triangle Ckt and changing the slope of a single Ck,i k,i
exclusively the coefficients ctk,i,d , d = 1, 2, and geometrically regulates the differential of s at Wk, i in the direction of wk, i (Fig. 8c
t , i = 1, 2, 3, and modifying C t on the other hand, one intuitively affects the graph surface in the and d). By fixing the slopes of Ck,i k
neighbourhood of s(Zk ) since Ckt is tangent to s at Zk (Fig. 8d and e). Raising and lowering of the fixed object gives the effect of extending and shrinking the spline locally. For a boundary edge ε Vi , Vj of , the control line Ce Q e,1 , Q e,2 provides a geometric tool similar to a control triangle. Through changing Ce above Le one may effectively modify the shape of the graph alongside the boundary of the domain (Fig. 8f).
126
J. Grošelj, M. Krajnc / Applied Mathematics and Computation 272 (2016) 114–126
The componentwise extension of spline functions to parametric spline surfaces is straightforward. The parametrization of the spline surface is given by (18) with the control points being considered as arbitrary points in the space. In this way, the control triangles and lines become independent of the domain and can be modified freely. The representation of parametric spline surfaces in terms of the introduced basis functions also enables a further generalization to rational spline surfaces, which indicates a rich potential of cubic splines on Powell–Sabin triangulations in computer aided geometric design. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]
M.J. Lai, L.L. Schumaker, Spline Functions on Triangulations, Cambridge University Press, Cambridge, 2007. M.J.D. Powell, M.A. Sabin, Piecewise quadratic approximations on triangles, ACM Trans. Math. Softw. 3 (1977) 316–325. P. Dierckx, On calculating normalized Powell–Sabin B-splines, Comput. Aided Geom. Des. 15 (1997) 61–78. H. Speleers, Construction of normalized B-splines for a family of smooth spline spaces over Powell–Sabin triangulations, Constr. Approx. 37 (2013) 41–72. O. Davydov, P. Sablonnière, C 2 piecewise cubic quasi-interpolants on a 6-direction mesh, J. Approx. Theory 162 (2010) 528–544. S. Remogna, Bivariate C 2 cubic spline quasi-interpolants on uniform Powell–Sabin triangulations, Adv. Comput. Math. 36 (2012) 39–65. S.-K. Chen, H.-W. Liu, A bivariate C 1 cubic super spline space on Powell–Sabin triangulation, Comput. Math. Appl. 56 (2008) 1395–1401. M. Lamnii, H. Mraoui, A. Tijini, A. Zidna, A normalized basis for C 1 cubic super spline space on Powell–Sabin triangulation, Math. Comput. Simul. 99 (2014) 108–124. A. Lamnii, M. Lamnii, H. Mraoui, Cubic spline quasi-interpolants on Powell–Sabin partitions, BIT Numer. Math. 54 (2014) 1099–1118. D. Sbibih, A. Serghini, A. Tijini, A. Zidna, Superconvergent C 1 cubic spline quasi-interpolants on Powell–Sabin partitions, BIT Numer. Math. (2014), doi:10.1007/s10543-014-0523-z. C. Manni, P. Sablonnière, Quadratic spline quasi-interpolants on Powell–Sabin partitions, Adv. Comput. Math. 26 (2007) 283–304. D. Sbibih, A. Seghini, A. Tijini, Polar forms and quadratic spline quasi-interpolants over Powell–Sabin partitions, Appl. Numer. Math. 59 (2009) 938–958. D. Sbibih, A. Serghini, A. Tijini, Superconvergent quadratic spline quasi-interpolants on Powell–Sabin partitions, Appl. Numer. Math. 87 (2015) 74–86. H. Speleers, P. Dierckx, S. Vandewalle, Numerical solution of partial differential equations with Powell–Sabin splines, J. Comput. Appl. Math. 189 (2006) 643–659. H. Speleers, P. Dierckx, S. Vandewalle, Local subdivision of Powell–Sabin splines, Comput. Aided Geom. Des. 23 (2006) 446–462. H. Speleers, P. Dierckx, S. Vandewalle, Quasi-hierarchical Powell-Sabin B-splines, Comput. Aided Geom. Des. 26 (2009) 174–191. H. Speleers, C. Manni, F. Pelosi, M.L. Sampoli, Isogeometric analysis with Powell–Sabin splines for advection-diffusion-reaction problems, Comput. Methods Appl. Mech. Eng. 221–222 (2012) 132–148. E. Vanraes, J. Windmolders, A. Bultheel, P. Dierckx, Automatic construction of control triangles for subdivided Powell–Sabin splines, Comput. Aided Geom. Des. 21 (2004) 671–682. J. Maes, E. Vanraes, P. Dierckx, A. Bultheel, On the stability of normalized Powell–Sabin B-splines, J. Comput. Appl. Math. 170 (2004) 181–196. V. Klee, M.C. Laskowski, Finding the smallest triangles containing a given convex polygon, J. Algorithms 6 (1985) 359–375. J. O’Rourke, A. Aggarwal, S. Maddila, M. Baldwin, An optimal algorithm for finding minimal enclosing triangles, J. Algorithms 7 (1986) 258–269.