A hole filling method for explicit and parametric surfaces by using C1 -Powell Sabin splines

Available online at www.sciencedirect.com

Mathematics and Computers in Simulation 99 (2014) 71–81

Original article

A hole filling method for explicit and parametric surfaces by using C1-Powell Sabin splines M.A. Fortes, P. González ∗ , M. Pasadas, M.L. Rodríguez Departamento de Matemática Aplicada, Universidad de Granada, E.T.S.I.C.C.P, Avda. Severo Ochoa, s/n, 18071 Granada, Spain Received 31 October 2011; received in revised form 19 October 2012; accepted 17 April 2013 Available online 13 May 2013

Abstract In this work we develop a method to fill a hole in a surface, either explicit or in parametric form, or just in a set of three dimensional scattered data. We will construct a new surface which is very close to the original one where it is known and that fills the hole in a homogeneous way, in such a way that the final reconstruction is of class C1 . We give results which prove the existence and uniqueness of solution of the proposed method, and we present several graphical examples which show the efficiency of the theory developed. © 2013 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Hole filling; Powell-Sabin; Splines; Parametric surfaces

1. Introduction In the past decades, variational methods in CAD, CAGD and Earth Sciences have received considerable attention (see for example [3,4]. A great effort to develop the theory and improve their efficiency have been also done, due to their usefulness in the fitting and design of curves and surfaces. Special mention may have the works [1,2,13]. The basic idea behind these methods is to minimize a functional which typically contains two terms: the first one indicates how well the curve or surface approximates a given data set, while the second one controls the degree of smoothness or fairness of the curve or surface. A wide range of minimization functionals have been proposed, derived from physical considerations: stretch and/or bending energy, etc.; or geometric entities: curve length, curvature, surface area, etc. (see for example [5,12,6,14]). Discrete smoothing variational splines and discrete smoothing Dm -splines provide specific examples of variational curves and surfaces. In [16] a functional of the above type is minimized in a parametric space of bicubic splines. Moreover, in all cases the obtained splines approximate a Lagrangian or Hermite data set. Other papers related to this matter are [12,14] and references therein. In this paper we treat now the problem of adequately filling the holes of 3D-surfaces, not necessarily explicit, and even closed in parametric form, as spheroid type ones: that is, such obtained by a radial function of the spherical coordinates. Several works related to the field of filling holes have been published in the last years (see [7,10,15,20]).

∗

Corresponding author. Tel.: +34 958246190; fax: +34 958249513. E-mail addresses: [email protected], [email protected] (P. González).

0378-4754/$36.00 © 2013 IMACS. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.matcom.2013.04.013

72

M.A. Fortes et al. / Mathematics and Computers in Simulation 99 (2014) 71–81

Fig. 1. 3D and 2D projected points.

Also, the problem of interpolation or approximation by radial basis functions (RBF) may have some basic points in common with this technique but the final procedure is completely different. The RBF theory have been treated since some time ago, see for example [18] and still continue having a great interest now. For a given data set of 3D-points in the space, belonging or obtained from a given explicit or parametric surface, in any case we treat to find the spherical coordinates of all of these points respect to a specific point, that will be the corresponding centroid of the data set. Afterwards, all these 3D-points are projected over a least squares plane, that is appropriately rotated to find the correct (ϕ, θ) angles and the r (ϕ, θ) associated to every original point in the data set. All these steps can be summarised graphically in the following pictures for the Franke’s function example, see Figs. 1 and 2, where we have taken only 50 points in order to improve the visualization of the process (in black color are the original 3D-points and in green it appears the centroid, that also belong to the least squares plane in this particular case, where they also are the projected points, in red). But the same procedure can be done with all the examples presented in this paper, both in the explicit and parametric cases, where spheroidal surfaces are mainly considered. This particular process is just one of the multiple possibilities of finding a set of 2D points of parameters (u, v) over which we have (x (u, v) , y (u, v) , z (u, v)) values for all of them that finally corresponds to the 3D-points in the data set. We have been trying with other procedures but finally we have chosen this one, mainly for its simplicity and validity both in the case of explicit or parametric functions of spheroidal type. In a more general case, we consider a set of 3D-points that we want to parametrize and then approximate. This set of 3D points may be given previously from empirical measurements or from a given parametric surface that we want

Fig. 2. Spherical coordinates w.r.t. the centroid.

M.A. Fortes et al. / Mathematics and Computers in Simulation 99 (2014) 71–81

73

Fig. 3. Lift-up of the triangulation of the first patch.

to approximate. To this end, we begin by choosing an adequate subset of points sufficiently close to each other, that can be projected over a concrete least squares approximation plane in such a way that we can triangulate the projected points, maintaining their original topological disposition in 3D. This first subset of points are usually chosen to be near the center of gravity of the whole set of points. We triangulate then the corresponding convex-hull of this set of projected points in the plane by a standard DelaunayVoronoi triangulation procedure and, afterwards, it is lifted again into 3D. In this way, each patch can have its proper approximation considering the projection over an appropriate plane. We can see in the Fig. 3 the lifted triangulation in 3D from the corresponding least square plane of the first patch. In this way we have the first patch of the surface by a local parametrization of the first subset of points. Once the first patch of the surface is lifted up. The other patches are being added considering further subsets of points that are near the points of the contour of the previous patches into 3D (see Fig. 4). And we have to continue adding new subsets of points in the original set of data until the total set of points in the original data set is considered. The new patches, are then triangulated and approximated in the same way, but this general procedure is much more complicated that the previous ones and require a lot of programming skills to avoid problems, so the final version of the whole algorithm is still under development. The paper is organized as follows. In Section 2 we recall the main basic concepts and preliminaries to be used through the work. In Section 3 we formulate the problem and we prove the existence and uniqueness of solution for the filling method proposed. Finally, in the last section, we give several graphical examples.

Fig. 4. Adjunction of a new patch.

74

M.A. Fortes et al. / Mathematics and Computers in Simulation 99 (2014) 71–81

2. Basic concepts and preliminaries Let D ⊂ R2 be a polygonal domain (an open polygonal and nonempty connected set) and let us consider the Sobolev space H2 (D; R3 ) = {f : D → R3 : ∂β f ∈ L2 (D; R3 ), ∀|β| ≤ 2}, where β = (β1 , β2 ) ∈ N2 and |β| = β1 + β2 . That is, f ∈ H2 (D; R3 ) if and only if f = (f1 , f2 , f3 ), where fi : D → R satisfies ∂β fi ∈ L2 (D) for |β| ≤ 2 and i = 1, 2, 3 . Let  · , · 3 and  · 3 denote the usual Euclidean inner product and norm in R3 , respectively, and, for any open subset X ⊂ D, consider the inner semi-products   ∂β f(x), ∂β g(x) 3 dx, m = 0, 1, 2, (f, g)m,X = |β|=m X

the semi-norms |f|m,X =

⎛

1/2 (f, f)m,X

=⎝

|β|=m X

and the norm 

f X =

⎛

1/2

2 

⎞1/2

 

|f|2m,X

=⎝

∂β f(x) 23 dx⎠



|β|≤2 X

m=0

m = 0, 1, 2,

,

⎞1/2 ∂β f(x) 23 dx⎠

.

Let T be a triangulation of D and consider the associated Powell-Sabin triangulation T6 of T (see e.g. [17]) obtained by joining an interior split point ΩT (for example the incenter of T) in each interior triangle T ∈ T to the vertices of T and to the split points ΩT  of the neighbouring triangles T  ∈ T. When T has a side lying on the boundary of D, the point ΩT is joined to the mid-point of this side, to the vertices of T and to the split points ΩT  of the neighbouring triangles T  ∈ T. We want to mention that in the particular case in which T be an α-triangulation (see e.g. [17]), a convergence result similar to the one given in [7] for the continuous filling can be stated. Let us consider the splines space 

S12 (D, T6 ) = s ∈ C1 (D) : s|T  ∈ P2 (T  ) ∀T  ∈ T6 , where P2 (T  ) denotes the space of bivariate polynomials of total degree at most two over T . In [19] it is shown that given the values of a function f (defined on D) and the values of the first partial derivatives of f at all the vertices of T, there exists a unique s ∈ S12 (D, T6 ) such that the values of s and the values of its first partial derivatives coincide with those of f at all the vertices of T (Fig. 5). 3. Formulation of the problem Let S ⊂ R3 be a parametric surface with a hole H3D in it. Maybe we know the expression of a parametric function ϒ(u, v) = (x(u, v), y(u, v), z(u, v)) whose graphic is S, or maybe we just know some scattered data point set P = {qi = (xi , yi , zi )}ki=1 in S (see Fig. 6). Let P2D = {(ui , vi )}ki=1 be the parameters values associated to the points of P, that is, ϒ(ui , vi ) = pi (in the final section of this work we will develop a method to obtain such parametrization when only the scattered data point set P is known). T of D, and the Powell-Sabin Let us consider now a polygonal domain D containing the set P2D , a triangulation triangulation T6 associated to T. Let TH = {T ∈ T : T ∩ P2D = ∅} and H ∗ = T ∈TH T. Consider the finite element spaces WD−H ∗ = {v|

◦ D−H

∗

: v ∈ S12 (D, T6 ), v|

◦ ∗ D−H

= / 0}

and

M.A. Fortes et al. / Mathematics and Computers in Simulation 99 (2014) 71–81

75

Fig. 5. Powell-Sabin sub-triangulation.

WH ∗ = {v|H ∗ : v ∈ S12 (D, T6 ), v|H ∗ = / 0}.

be the subset of P consisting of the points such that their parameters values associated lay in D − H* , that is, Let P

= {pi ∈ P : (ui , vi ) ∈ D − H ∗ }, P 2D

and let P

For the sake of simplicity we ⊂ D − H ∗ the set of the parameters values associated to the points in P. 2D k

as {qi = (xi , yi , zi )}

as {(ui , vi )}k . will represent the points in P i=1 and the points in P i=1 The reconstruction process will be done in two steps:

outside the hole H3D in S; Step 1: We will define a parametric surface approximating the data set points P 3D Step 2: We will reconstruct the surface inside the hole H in such a way that the global reconstruction of S will be of class C1 . Let us develop the two stages of the filling method: Step 1. Let Rk,3 be the space of real matrices of k rows and 3 columns equipped with the inner semi-product A, B k,3 =

k  3 

aij bij

i=1 j=1

Fig. 6. Parametric surface and scattered data set with holes.

76

M.A. Fortes et al. / Mathematics and Computers in Simulation 99 (2014) 71–81

and the norm A k,3

⎞1/2 ⎛ k  3  =⎝ aij2 ⎠ , i=1 j=1

where A = (aij ) 1 ≤ i ≤ k and B = (bij ) 1 ≤ i ≤ k . 1≤j≤3 1≤j≤3 Let ρ be the evaluation operator ρ : H2 (D − H ∗ ; R3 ) −→ Rk,3 defined by ⎛ ⎞ f1 (u1 , v1 ) f2 (u1 , v1 ) f3 (u1 , v1 )  ⎜  ⎟ .. .. .. ⎟. ρ(f) = (f(ui , vi ))ki=1 =⎜ ⎝ ⎠ . . . f1 (uk , vk )

f2 (uk , vk )

f3 (uk , vk )

Let us suppose that P1 (D − H ∗ ; R3 ) ∩ Ker(ρ) = {0}.

(1)

Let us consider the functional energy J1 : H (D − H ∗ ; R3 ) −→ R defined by 2

J1 (f) = ρ(f) − ((qi )ki=1 ) k,3 + τ1 |f|21,D−H ∗ + τ2 |f|22,D−H ∗ , where τ 1 ≥ 0 and τ 2 > 0 .

2D approximates the values of P

(in the least Observe that the first summand of J1 measures how well f over P squares sense), while the last two summands represent the “minimal energy condition” over the semi-norms | · |1 and | · |2 , both weighted by the parameters τ 1 and τ 2 , respectively. This functional allows us to consider the following problem: Problem P1 We look for an element σ 1 ∈ (WD−H ∗ )3 such that minimizes the functional J1 , that is, J1 (σ 1 ) ≤ J1 (f) for all f ∈ (WD−H ∗ )3 . Proposition 1. Problem P1 has a unique solution which is also the unique solution to the following variational problem: Find σ 1 ∈ (WD−H ∗ )3 such that for all f ∈ (WD−H ∗ )3 it holds ρ(σ 1 ), ρ(f) k,3 +

2 

τm (σ 1 , f)m,D−H ∗ = ((qi )ki=1 ) , ρ(f) k,3 .

(2)

m=1

Proof. By using (1), it can be shown that the application defined in (WD−H ∗ )3 by 1/2  2  2 [[v]] = ρ(v) k,3 + τm |v|m,D−H ∗ m=1

is a norm in (WD−H ∗ )3 equivalent to the usual norm f D−H ∗ . As a consequence, the symmetric continuous bilinear form a : (WD−H ∗ )3 × (WD−H ∗ )3 −→ R given by a(u, v) = ρ(u), ρ(v) k,3 +

2 

τm (u, v)m,D−H ∗

m=1

is (WD−H ∗ )3 -elliptic. Besides, the application ϕ : (WD−H ∗ )3 −→ R defined by ϕ(v) = ρ(v) − ((qi )ki=1 ) k,3 is a linear and continuous form, and hence, we conclude by applying Lax-Milgram Lemma.

outside the hole H3D . σ 1 will be the parametric surface used to approximate the set of points P

M.A. Fortes et al. / Mathematics and Computers in Simulation 99 (2014) 71–81

77

The computation of σ 1 follows from Eq. (2): Let N = dim (WD−H ∗ ) and B = {B1 , . . . , BN } be a basis of WD−H ∗ whose elements have local support. Then, dim((WD−H ∗ )3 ) = 3N and {ej Bi } 1 ≤ i ≤ N is a basis of (WD−H ∗ )3 , where {e1 , e2 , e3 } is the canonical 1≤j≤3 R3 .

basis of  3 Let us denote σ 1 = N i=1 αi Bi , where αi ∈ R , to the unique solution of Problem P1 . Then from Eq. (2) we obtain N that α = (αi )i=1 is the unique solution to the linear vectorial system ⎞ ⎛ α1 ⎟ ⎜ . ⎟ C⎜ (3) ⎝ .. ⎠ = B, αN



(Bj ) k + 2m=1 τm (Bi , Bj )m,D−H ∗ , ρ(Bi ), ρ where C = (cij )i,j=1,...,N , cij =  ⎞ ⎛

(B1 ) k (ys )ks=1 , ρ

(B1 ) k (zs )ks=1 , ρ

(B1 ) k (xs )ks=1 , ρ ⎟ ⎜ .. .. .. ⎟, B=⎜ ⎠ ⎝ . . . k k k

(BN ) k (ys )s=1 , ρ

(BN ) k (zs )s=1 , ρ

(BN ) k (xs )s=1 , ρ  

is the operator defined in WD−H ∗ by ρ

(f ) = f (ui , vi )ki=1 . and ρ It can be shown that the coefficient matrix C is symmetric, positive definite and of band type (in [6], Section 4, it is explained how the basis {Bi }N i=1 can be constructed and in Lemma 2 it is proved why C is symmetric, positive definite and of band type). Remark. We could consider optimum values for the parameters τ = (τ 1 , τ 2 ) those that the corresponding filling be as close as possible to the original data set P. To this end, we can proceed as follows: first, we consider the function  3 element N j=1 α0j Bj ∈ (WD−H ∗ ) which minimizes the function  2  q    N    R(α1 , . . . , αN ) =  αj Bj (ui , vi ) − pi  .   i=1 j=1 3 ∗ ∗ = (τ1 , τ2 ) such that the solution of the vectorial system (3) with these parameters , then we will take τ * as optimum parameters. But, in general, there would not exist parameters

If there exist parameters values τ ∗

values be α0 = (α0j )N j=1 such that the solution of the vectorial system (3) be α0 . So, let fτ ∈ (WD−H ∗ )3

 τ τ τ be the filling function obtained for the parameters values τ = (τ 1 , τ 2 ), and denote fτ = N j=1 αj Bj and α = αj . We will * τ consider as optimum parameters τ those minimizing the function φ(τ) = < α0 − α > N . To minimize this function we can use the Powell’s algorithm ([8], Section 7.3) to find the minimum of a vectorial function without using derivatives. Step 2. Let N = {t1 . . . , tn } be the set of all vertices of T laying on the boundary of H* and let us consider the map ϕ(v) = (ϕi (v))3n i=1

for v ∈ (WH ∗ )3 ∪ (WD−H ∗ )3 ,

where ⎧ ϕi (v) = v(ti ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ∂v ϕn+i (v) = (ti ) for i = 1, . . . , n; ∂x ⎪ ⎪ ⎪ ∂v ⎪ ⎪ ⎩ ϕ2n+i (v) = (ti ) ∂y

78

M.A. Fortes et al. / Mathematics and Computers in Simulation 99 (2014) 71–81

the subset Wσ 1 = {v ∈ (WH ∗ )3 : ϕ(v) = ϕ(σ 1 )} and the functional J2 : (WH ∗ )3 −→ R defined by J2 (f) = τ1 |f|1,H ∗ + τ2 |f|2,H ∗ , with τ 0 > 0 . We want to mention that although the parameters involved in the functional J2 could be different from the ones involved in J1 , we have chosen them to coincide just to consider the same weight for the seminorms outside and inside H* . The minimization problem we want to solve is: Problem P2 . We look for an element σ 2 ∈ Wσ 1 such that J2 (σ 2 ) ≤ J2 (v) for all v ∈ Wσ 1 . Proposition 2. Problem P2 has a unique solution which is also the unique solution to the following variational problem: Find σ 2 ∈ Wσ 1 such that 2 

τm (σ 2 , w)m,H ∗ = 0

forall w ∈ W0 = {v ∈ (WH ∗ )3 : ϕ(v) = 0}.

m=1

The proof of Proposition 2 follows the same pattern that the one of Theorem 3 of [7] by taking ((u, v)) =

2 

< ϕin+1 (u), ϕin+1 (v)>3 +

i=0

2 

τm (u, v)m,H ∗ .

m=1

σ 2 will be the parametric surface used to reconstruct the original surface inside the hole H3D . In order to compute σ 2 , let N0 = dim(WH ∗ ). Let us extend the set of basis functions {ωi }3n i=1 introduced in the proof 3 N 0 of Proposition 2 to the usual Hermite basis {ωi }i=1 of WH ∗ . Then, dim((WH ∗ ) ) = 3N0 and {ej ωi } 1 ≤ i ≤ N is a 0

1≤j≤3 basis of (WH ∗ )3 , where {e1 , e2 , e3 } is the canonical basis of R3 . Then {ej ωi } 3n + 1 ≤ i ≤ N is a basis of W0 and hence, 0

1≤j≤3 σ2 =

3n 

ϕi (σ 1 )ωi +

i=1

N0 

cj ω j .

j=3n+1

0 By applying Proposition 2 we get that the matrix (cj )N j=3n+1 is the solution of the linear vectorial system

N0  j=3n+1



2 

τm (ωj , ωt )m,H ∗

m=1

cj = −

3n  i=1

for all t = 3n + 1, . . ., N0 . The global C1 -reconstruction is defined as  σ1 (u, v) if (u, v) ∈ D − H ∗ . σ2 (u, v) if (u, v) ∈ H ∗

 ϕi (σ 1 )

2 

m=1

τm (ωi , ωt )

◦ ∗ m,H

M.A. Fortes et al. / Mathematics and Computers in Simulation 99 (2014) 71–81

79

Fig. 7. Scattered data points from a semi-sphere.

4. Graphical examples In the first example we have considered a set of 900 points on a semi-sphere in which we have done a hole consisting of 90 points, like it is shown in Fig. 7. We have considered a uniform partition of the projection domain D into 128 triangles and the smoothing parameters τ 1 = 10−9 , τ 2 = 10−5 and τ 0 = 10−5 . In Fig. 8 we show the graphics obtained in Step 1 and in Step 2 of the process. In the second example we have considered a set of 2000 points on a whole sphere in which we have done a hole consisting of 100 points. In Fig. 9 we show the set P2D of the projection of the scattered points domain onto the domain plane which lay out of the elliptic projection H of the hole H3D . We have considered a uniform partition of the projection domain D into 162 triangles and the smoothing parameters τ 1 = 10−5 , τ 2 = 10−3 and τ 0 = 10−5 . In Fig. 10 we show the graphics obtained in Step 1 and in Step 2 of the process. Next, we have considered a set of 2000 points in a spheroidal surface, in which we have done a hole consisting of 100 points. We have considered a uniform partition of the projection domain D into 162 triangles and the smoothing parameters τ 1 = 10−5 , τ 2 = 10−3 and τ 0 = 10−5 . In Fig. 11 we show the graphics obtained in Step 1 and in Step 2 of the process. Finally, we have considered a set of 2000 points in the surface given by the graphic of Franke’s function, defined by   (9u − 2)2 + (9v − 2)2 34 34 9u + 1 (9v + 1)2 f (x, y) = − + − + exp 4 exp 10 49  .  12 (9v − 7)2 + (9u − 3)2 15  2 2 − − −(9v − 4) + (9u − 7) exp 4 exp

Fig. 8. Hole reconstruction in a semi-sphere.

80

M.A. Fortes et al. / Mathematics and Computers in Simulation 99 (2014) 71–81

Fig. 9. Hole in the projected points.

Fig. 10. Hole reconstruction in a sphere.

Fig. 11. Hole filling in a spheroidal surface.

M.A. Fortes et al. / Mathematics and Computers in Simulation 99 (2014) 71–81

81

Fig. 12. Scattered 3D points from Franke function and hole reconstruction.

In Fig. 12 we show the graphic of the scattered data in which we have done a hole consisting of 1000 points on the left. We have considered a uniform partition of the projection domain D into 162 triangles and the smoothing parameters τ 1 = 10−6 , τ 2 = 10−3 and τ 0 = 10−5 . In the middle and right of this figure we show the graphics obtained in Step 1 and in Step 2 of the process. Acknowledgements This work has been supported by the Junta de Andalucia (Research Group FQM/191) and by the Ministerio de Educación y Ciencia of Spain (Research Projects MTM2005-01403 and MTM2008-00671). References [1] D. Apprato, R. Arcangèli, R. Manzanilla, Sur la construction de surfaces de classe Ck à partir d’un grand nombre de données de Lagrange (French).[On the construction of Ck -surfaces starting from a large number of Lagrange data], RAIRO. Modélisation Mathématique et Analyse Numérique 21 (4) (1987) 529–555. [2] D. Apprato, C. Gout, Approximation de surfaces à partir de morceaux de surfaces., Comptes Rendus de l’Académie des Sciences, Paris, Série 1 325 (4) (1997) 445–448. [3] D. Apprato, C. Gout, D. Komatitsch, A new method for Ck surface approximation from a set of curves, with application to ship track data in the Marianas trench, Mathematical Geology 34 (7) (2002) 831–843. [4] D. Apprato, C. Gout, P. Sénéchal, Ck -Reconstruction of surfaces from partial data, Mathematical Geology 32 (8) (2000) 969–983. [5] R. Arcangèli, M.C. López de Silanes, J.J. Torrens, Multidimensional Minimizing Splines, Kluwer Academic Publisher, Boston, 2004. [6] D. Barrera, M.A. Fortes, P. González, M. Pasadas, Minimal energy surfaces on Powell-Sabin type triangulations, Applied Numerical Mathematics 58 (5) (2008) 635–645. [7] D. Barrera, M.A. Fortes, P. González, M. Pasadas, Filling polygonal holes with minimal energy surfaces on Powell-Sabin type triangulations, Journal of Computational and Applied Mathematics 234 (4) (2010) 1058–1068. [8] R. Brent, Algorithms for Minimization Without Derivatives, Dover Publications, New York, 2002. [10] C.K. Chui, M-J. Lai, Filling polygonal holes using C1 -cubic triangular patches, Computer Aided Geometric Design 17 (2000) 297–307. [12] J. Duchon, Interpolation des fonctions de deux variables suivant le principle de la flexion des plaques minces, RAIRO 10 (12) (1976) 5–12. [13] C. Gout, Ck surface approximation from surface patches, Computers & Mathematics with Applications 44 (3–4) (2002) 389–406. [14] G. Greiner, Surface construction based on variational principles., in: P.J. Laurent, A. Le Méhauté, L.L. Schumaker (Eds.), Wavelets, Images and Surface Fitting, AK Peters, Wellesley, Boston, 1994, pp. 277–286. [15] J.A. Gregory, J. Zhou, Filling polygonal holes with bicubic patches, Computer Aided Geometric Design 11 (1994) 391–410. [16] A. Kouibia, M. Pasadas, Approximation by discrete variational splines, Journal of Computational and Applied Mathematics 116 (2000) 145–156. [17] M. Laghchim-Lahlou, P. Sablonnière, Cr -finite elements of Powell-Sabin type on the three direction mesh, Advances in Computational Mathematics 6 (1996) 191–206. [18] M.J.D. Powell, The theory of radial basis function approximation in 1990., in: W. Light (Ed.), Advances in Numerical Analysis, Vol II, Oxford Science Publications, Oxford, 1992, pp. 105–210. [19] M.J.D. Powell, M.A. Sabin, Piecewise quadratic approximations on triangles, ACM Transactions on Mathematical Software 3 (4) (1977) 316–325. [20] V. Rayevskaya, L.L. Schumaker, Multi-sided macro-element spaces based on Clough–Tocher triangle splits with applications to hole filling, Computer Aided Geometric Design 22 (1) (2005) 57–79.