This study deals with constructing a convexity-preserving bivariate C1 interpolants to scattered data whenever the original data are convex. Sufficient conditions on lower bound of Bezier points are derived in order to ensure that surfaces comprising cubic Bezier triangular patches are always convex and satisfy C1 continuity conditions. Initial gradients at the data sites are estimated and then modified if necessary to ensure that these conditions are satisfied. The construction is local and easy to be implemented. Graphical examples are presented using several test functions

Abd Majid, Ahmad

Mt Piah, Abd Rahni

Saaban, Azizan

UUM Repository

MATEMATIKA, 2008, Volume 24, Number 1, 31–42c©Department of Mathematics, UTM.Convexity-Preserving Scattered Data Interpolation1Abd. Rahni Mt. Piah, 2Azizan Saaban & 3Ahmad Abd. MajidSchool of Mathematical Sciences, Universiti Sains Malaysia 11800 USM Pulau Pinang, Malaysiae-mail: 1arahni@cs.usm.my, 2azizan.s@uum.edu.my, 3majid@cs.usm.myAbstract This study deals with constructing a convexity-preserving bivariate C1interpolants to scattered data whenever the original data are convex. Sufficient con-ditions on lower bound of Bézier points are derived in order to ensure that surfacescomprising cubic Bézier triangular patches are always convex and satisfy C1 continuityconditions. Initial gradients at the data sites are estimated and then modified if nec-essary to ensure that these conditions are satisfied. The construction is local and easyto be implemented. Graphical examples are presented using several test functions.Keywords Scattered data, interpolation, convexity, continuity.1 IntroductionIn Computer Aided Geometric Design, we often need to consider designing interpolatingsurfaces with convexity or concavity properties. Several authors ([4], [5], [2], [10]) havederived sufficient conditions to ensure the convexity of either triangular or tensor-productBernstein-Bézier surfaces. Examples of the methods used to construct convexity-preservinginterpolating surfaces are given in ([1], [6], [10], [12], [13], [15]). [1] derived the convexity-preserving conditions for piecewise quadratic C1 interpolation based on Powell-Sabin splitsbut are not complete as shown by a counter example in [13]. The complete set of condi-tions on convexity-preserving interpolation by piecewise quadratic and cubic C1 functionswas then derived by [13] by combining local convexity conditions with an optimizationscheme which takes into account the global nature of the problem. [10] used bivariateC1 cubic splines to deal with convexity-preserving scattered data interpolation. A neces-sary and sufficient condition on Bézier-Bernstein polynomials were derived. The problemof convexity-preserving is set into quadratically constraints in a quadratic programmingproblem. Then the quadratic constraints were replaced by three linear constraints and theproblem was formulated into linearly constraint quadratic programming. [12] described aC1 convexity-preserving scattered data interpolation. The method consists of constructinga triangulation of the nodes for which the triangular-based piecewise linear interpolant isconvex, computes a set of nodal gradients for which a convex Hermite interpolant exists,and constructs a smooth convex surface that interpolates the nodal values and gradients.It involves two data-dependent triangulations with a straight-line dual of each. [6] de-fined new triangular macro-elements based on Clough-Tocher cubic elements which allowthe construction of a C1 surface connecting patches originally of non-equal degree. The32 Abd. Rahni Mt. Piah, Azizan Saaban & Ahmad Abd. Majidmain tool in this construction relies on the degree elevation process associated with care-ful modification of the Bézier nets. In addition, [6] proposed a shape-preserving criterionfor automatic selection of the tension parameters, which permits preserving monotonicityand/or convexity of the data in domains containing, and along directions parallel to theedges of the triangulation.Motivated by sufficient condition for convexity of Bernstein-Bézier surfaces over trianglesgiven in [2] and [10], this paper will propose sufficient conditions by imposing a lower boundon Bézier points in order to ensure that surfaces comprising cubic Bézier triangular patchesare always convex if given data are convex and satisfy C1 continuity conditions. Eachtriangular patch of the interpolating surface is formed as a convex combination of three cubicBézier triangular patches. This method is local and the convexity-preserving interpolatingsurfaces are controlled by certain free parameters. In Section 2, we present some technicalprerequisites concerning convex data, the Bernstein-Bézier representation of polynomials,the convexity conditions of triangular Bézier surfaces and the sufficient conditions for C1convexity-preserving interpolant. An outline of the surface construction process is given inSection 3, while Section 4 presents some examples. Finally, the conclusion will be given inSection 5.2 Background2.1 Convex DataScattered data, (xi, yi, zi), i = 1, 2, ..., N , are called strictly convex if there exists a strictlyconvex function F such that F (xi, yi) = zi andF (µp + (1− µ)q) = µF (p) + (1− µ)F (q), ∀µ ∈ [0, 1] (1)for every pair of points p, q ∈ D where D is the convex hull of the nodes. Equivalently,the data set is convex if there is a triangulation of the nodes for which the triangular-basedpiecewise linear interpolant is convex (Such a triangulation is unique and the data are saidto be strictly convex if no four data points are coplanar [12]). An algorithm to constructsuitable triangulation T for convex data is given in [16].2.2 Triangular Bézier PatchA triangular Bernstein-Bézier surface of degree n over triangle T =< V1, V2, V3 > is definedasP (x, y) =∑i+j+k=ni≥0, j≥0, k≥0bijkBnijk(x, y) (2)where bijk are the Bézier ordinates or control points of P andBnijk(x, y) =n!i!j!k!risjtk (3)with (x, y) = rV1 + sV2 + tV3 and r, s, t , are barycentric coordinates (i.e. r + s + t = 1).Convexity-Preserving Scattered Data Interpolation 332.3 Convexity ConditionsThe convexity of a polynomial can be conveniently characterized through its Bézier rep-resentation. There are many sufficient conditions for convexity-preserving triangle in theliterature, such as in [1], [2], [4], [5], [10]. For our study, we use convexity-preserving condi-tions in [10]. For any direction given as vector d = (d1, d2, −(d1 + d2)), the second orderdirectional derivative is given byD2dP (x, y) = n(n− 1)∑i+j+k=n−2Cmi,j,k(d)Bn−2ijk (x, y) (4)where Bnijk(x, y) ) is defined in (3) andCmi,j,k(d) = d1(Cm−1i+1,j,k(d)− Cm−1i,j,k+1(d) ) + d2(Cm−1i, j+1,k(d) − Cm−1i,j,k+1(d) ) (5)for m ≥ 1, with C0i,j,k(d) = bijk, i + j + k = n.Using (4) and (5) on a triangular patch represented by (2), we obtainD2dP (x, y) = n(n− 1)∑i+j+k=n−2(d21(bi+2, j, k + bi,j,k+2 − 2bi+1, j, k+1)+2d1d2(bi, j, k+2 + bi+1, j+1, k − bi+1, j, k+1 − bi, j+1, k+1)+ d22(bi, j+2, k + bi,j,k+2 − 2bi,j+1,k+1))Bn−2ijk (x, y)= n(n− 1) ∑i+j+k=n−2(d1 d2)Q(i, j, k)(d1d2)Bn−2ijk (x, y)with Q(i, j, k), a 2 x 2 matrix defined byQ(i, j, k) =(∆1bijk ∆2bijk∆2bijk ∆3bijk)(6)where difference operators, ∆lbijk, l = 1, 2, 3 is defined as∆1bijk = bi+2,j,k + bi,j,k+2 − 2bi+1, j, k+1 ,∆2bijk = bi,j,k+2 + bi+1,j+1,k − bi+1, j, k+1 − bi, j+1, k+1and ∆3bijk = bi,j+2,k + bi,j,k+2 − 2 bi, j+1, k+1 for i + j + k = n − 2, i, j, k ≥ 0. Notethat P (x, y) is convex on T if and only if D2dP (x, y) ≥ 0 for all (x, y) ∈ T for all directiond. P (x, y) is convex when Q(i, j, k) in (6) is positive semi-definite for all i, j, k ≥ 0 withi + j + k = n− 2 which is equivalent to the conditions∆1bijk ≥ 0, ∆3bijk ≥ 0and(∆1bijk )(∆3bijk) ≥ (∆2bijk)2. (7)34 Abd. Rahni Mt. Piah, Azizan Saaban & Ahmad Abd. Majid2.4 C1 Cubic Convexity-Preserving InterpolantsA cubic Bézier triangular patch P on T is defined as,P (x, y) =∑i+j+k=3i≥0, j≥0, k≥0bijkB3ijk(x, y)= b300r3 + b030s3 + b003t3 + 3b210r2s + 3b201r2t + 3b120rs2 + 3b102rt2+3b021s2t + 3b012st2 + 6b111rst (8)where bijk are the Bézier ordinates of P as shown in Figure 1.Given the Bézier ordinates at vertices i.e. b300, b030, b003, the sufficient conditions on theremaining Bézier ordinates will be derived to ensure the convexity of the entire patch. LetA = b300, B = b030, and C = b003, A, B, C ∈ < . Note that in this paper we willonly focus on positive and convex data i.e. for the case of A > 0, B > 0 and C > 0.Our approach is to find lower bounds on the remaining Bézier ordinates, so that seconddirectional derivative, D2dP (x, y) ≥ 0 for all (x, y) ∈ T for any direction d. Let Bézierordinates b210, b201, b120, b021, b012, b102 have the same value−p and b111 = −q, with p, q > 0.Using (7), we obtain the following set of conditions for convexity-preserving on each trian-gular patch,−p + q ≥ 0, (9)p ≥ −A, p ≥ −B , p ≥ −C, (10)(−p + q)(2A − q + 3p) ≥ 0, (11)(−p + q)(2B − q + 3p) ≥ 0, (12)(−p + q)(2C − q + 3p) ≥ 0. (13)Convexity-Preserving Scattered Data Interpolation 35Since A > 0, B > 0, C > 0, equations (9) – (13) clearly hold when p = q = 0.From (10), we obtainp ≥ max(−A,−B,−C). (14)Thus, the lower bound of p is given byp0 = α max(−A, −B, −C) (15)where α is a free parameter.We need α < 0 in order for p0 > 0. If (9) holds, then conditions(11), (12) and (13) will just reduce to q ≤ min (3p + 2A, 3p + 2B, 3p + 2C).It is then easy to show that min (3p + 2A, 3p + 2B, 3p + 2C) ≥ p . Thusp ≤ q ≤ min (3p + 2A, 3p + 2B, 3p + 2C) . (16)From (15), when p0 is a lower bound of p, then (16) becomesp0 ≤ q0 ≤ min (3p0 + 2A, 3p0 + 2B, 3p0 + 2C) . (17)The lower bound of q can then be written asq0 = βp0 + (1− β) min (3p0 + 2A, 3p0 + 2B, 3p0 + 2C) (18)where 0 ≤ β ≤ 1 is a free parameter.PropositionConsider the cubic Bézier triangular patch P (x, y) with b300 = A, b030 = B, b003 = C. Ifb021, b012, b201, b102, b210, b120 ≥ −p0, and b111 ≥ −q0 ,where p0, q0 are given by (15) and(18) respectively then P (x, y) is convex on a given convex domain D ∈ <2.3 Construction of Convexity-Preserving Interpolating SurfaceWe want to construct a C1 convexity-preserving functional surface F (x, y) which interpo-lates a given positive and convex scattered data, (xi, yi, zi), i = 1, 2, ..., N, zi > 0. Thesurface comprises cubic Bézier triangular patches, each of which is guaranteed to remainconvex. We use Delaunay triangulation to triangulate the convex hull of the data points.Estimations of the first order partial derivatives of F with respect to x and y are obtainedusing the methods of quadratic approximation and least squares i.e.min∑i(ax2i + by2i + cxi + dyi + e− zi)2where a, b, c, d and e are coefficients to be determined.Then first order partial derivatives at vertex Vi can be estimated as Fx = 2ax + c andFy = 2by + d. Note that, in this paper we use a quadratic function without cross-productterm because it is suitable for convex data. Let Vi = (xi, yi), i = 1, 2, 3 be vertices of atriangle, such that F (Vi) = zi, and the first partial derivatives, Fx(Vi) and Fy(Vi). Foreach triangular patch P as in (8), the derivative along edge ejk joining (xj , yj) to (xk, yk)is given by∂P∂ejk= (xk − xj)∂F∂x+ (yk − yj)∂F∂y.36 Abd. Rahni Mt. Piah, Azizan Saaban & Ahmad Abd. MajidThe given data and estimated derivative values at (xi, yi), enable us to determine all theBézier ordinates brst except for b111. For example, we will have:b300 = F (V1), b210 = F (V1) +13∂F∂e12(V1) and b201 = F (V1) − 13∂F∂e31(V1).However, the initial estimate for each edge ordinate may not satisfy the convexity conditionsfor P . In view of the Proposition, we need these Bézier ordinates to be greater or equal to−p0. If they are not, then the magnitudes of Fx, Fy at the vertices need to be reduced sothat the convexity conditions are satisfied. The modification of these partial derivatives ata vertex Vi, is achieved by multiplying each derivative at that vertex, by a scaling factor0 < γ < 1. The smallest value of γ is obtained by considering all triangles that meet atvertex Vi, which satisfy the convexity conditions of all the triangles. For example,(b210)j = F (V1) + γ∂F/∂e123≥ (y3)jand(b201)j = F (V1) − γ∂F/∂e313≥ (y2)j ,where subscript j represents a quantity which corresponds to triangle j. Having adjustedthese derivatives, if necessary, the Bézier ordinates are recalculated using the formulaeabove. The above process is repeated at all the nodes, Vi (see [3] for further details).For each triangle, the inner Bézier ordinate b111, remains to be calculated, in such a way toguarantee the preservation of convexity and to ensure C1 continuity across patch boundaries.Let T1 and T2 be two adjacent cubic Bézier triangular patches with a common boundarycurve (see Figure 2). C1 continuity conditions between adjacent patches are given byc102 = µ1b120 + µ2b030 + µ3b021 (19)c111 = µ1b111 + µ2b021 + µ3b012 (20)c120 = µ1b102 + µ2b012 + µ3b003 (21)where (µ1, µ2, µ3) is the barycentric coordinates of W3 with respect to triangle T1. (19) and(21) are automatically satisfied since the Bézier ordinates b201, b210, b012, b102, b120, and b021have already been obtained.We shall use similar method as in [9] to determine initial value of inner Bézier ordinate,bi111 which is C1 across boundary ei, i = 1, 2, 3 where ei is the edge of the triangle oppositevertex Vi. Local scheme Pi is defined by replacing b111 in (8) with bi111. The values of innerBézier ordinates b1111, b2111 and b3111 for local schemes P1, P2, and P3 are given byb1111 =12(b120 + b102 +12(b012 + b021 − b003 − b030)),b2111 =12(b210 + b012 +12(b102 + b201 − b003 − b300))andb3111 =12(b021 + b201 +12(b120 + b210 − b300 − b030))Convexity-Preserving Scattered Data Interpolation 37respectively. Initial estimate of the inner Bézier ordinates as above in each triangle may notsatisfy the convexity conditions of P (x, y) as stated in the Proposition i.e. b111 ≥ −q0. Forexample, when e1 is a common edge to two triangles as in Figure 2, and the initial valueb1111 or c1111 is less than −q0, then these ordinates are both modified to ensure that b1111,c1111 ≥ −q0, and C1 continuity condition along the common boundary curve of adjacentpatches is satisfied. When e1 is on the boundary of the domain, and b1111 is less than −q0then b1111 is reset to be equal to q0. The ordinates b2111 and b3111 of P2 and P3 are adjustedaccordingly using a similar approach (see [3] for further details). The interpolating surfaceP on triangle T is then defined as a convex combination of all the local schemes such thatsufficient conditions on all sides of the triangles are satisfied, i.e.P (u, v, w) = c1P1(u, v, w) + c2P2(u, v, w) + c3P3(u, v, w)orP (u, v, w) =∑i+j+k=3i 6=1,j 6=1,k 6=1bijkB3ijk (u, v, w) + 6uvw(c1b1111 + c2b2111 + c3b3111)with c1 =vwvw + vu + uw, c2 =uwvw + vu + uwand c3 =vuvw + vu + uw, andu, v, w are the barycentric coordinates.4 ExamplesWe will illustrate our interpolating scheme using the following convex test functions:f(x, y) = x4 + y4 and g(x, y) = x3 + 5(y − 0.6)2 + 1.The first data set consists of eight points of convex triangulation defined on rectangular grid[−1, 1]× [−1, 1] are taken from [15] with data values taken from function f . Figure 3 showsa triangular domain of the data set. The corresponding linear interpolant is in Figure 4and the convexity-preserving interpolating surface of our method is given in Figure 5. Thesecond data set of nine data points are defined on rectangular grid [0, 1]× [0, 1] ([13]) with38 Abd. Rahni Mt. Piah, Azizan Saaban & Ahmad Abd. Majiddata values taken from function g. Figure 6 shows the triangular domain, Figure 7 is thecorresponding linear interpolant while Figure 8 shows the interpolating surface using ourconvexity-preserving scheme.The final data set consists of 26 scattered points from [11] with data values taken fromfunction f . Figures 9 to 11 show a triangular domain, linear interpolant and convexity-preserving interpolating surface respectively. Note that for the first and second data set, wechoose α = −0.028 and β = 0.2, while for the third data set we take α = −0.3 and β = 0.1as parameters to ensure convexity-preserved of the given functions.Convexity-Preserving Scattered Data Interpolation 3940 Abd. Rahni Mt. Piah, Azizan Saaban & Ahmad Abd. MajidConvexity-Preserving Scattered Data Interpolation 415 ConclusionWe have considered generating non-parametric surfaces which interpolate positive and con-vex scattered data using convex combination of cubic triangular patches. We have imposedlower bounds on Bézier points in order to ensure that surfaces comprising cubic Béziertriangular patches satisfy C1 continuity conditions and are always convex when the givendata are convex.Acknowledgements The authors are very grateful to the Ministry of Higher Educationand Universiti Sains Malaysia for providing us with the fundamental research grant scheme(FRGS) account number 203/PMATHS/671040 to enable us to pursue this research.References[1] J.M. Carnicer & W. Dahmen, Convexity preserving interpolation and Powell-Sabinelements, Computer Aided Geometric Design 9 (1992), 279-289.[2] J.M. Carnicer, M.S. Floater & J.M. Pena, Linear convexity conditions for rectangu-lar and triangular Bernstein-Bézier surfaces, Computer Aided Geometric Design 15(1997), 27-38.[3] E.S. Chan & B.H. Ong, Range restricted scattered data interpolation using convexcombination of cubic Bézier triangles, J. Comp. Appl. Math. 136 (2001), 135-147.[4] G. Chang & P.J. Davis, The convexity of Bernstein polynomials over triangles, J.Approx. Theory 40 (1984), 11-28.[5] G. Chang & Y.Y. Feng, An improved condition for the convexity of Bernstein-Béziersurfaces over triangles, Computer Aided Geometric Design 1(1984), 279-283.[6] P. Costantini & C. Manni , A local shape-preserving interpolation scheme for scattereddata, Computer Aided Geometric Design 16 (1999), 385-405[7] T.P. Fang & L.A. Piegl, Algorithm for Delaunay triangulation and convex-hull compu-tation using a sparse matrix. Computer Aided Design 24 (1992), 425-436.[8] G. Farin, Curves and Surfaces for Computer Aided Geometric Design-A PracticalGuide 4th Edition, Academic Press, San Diego (1996).[9] T.N.T. Goodman & H.B. Said, C1 triangular interpolant suitable for scattered datainterpolation, Communications in Applied Numerical Methods 7 (1991), 479-485.[10] M.J. Lai, Convex preserving scattered data interpolation using bivariate C1 cubicsplines, J. Comput. Applied Math. 119 (2000), 249-258.[11] C.L. Lawson, Software for C1 surface interpolation, in Mathematical Software III, J.R. Rice (ed)., Academic Press, New York (1977), 161-194.[12] N.K. Leung & R.J. Renka, C1 convexity-preserving interpolation of scattered data,SIAM J. Sci. Comput. 20(5) (1999), 1732-1752.42 Abd. Rahni Mt. Piah, Azizan Saaban & Ahmad Abd. Majid[13] A. Li, Convexity preserving interpolation, Computer Aided Geometric Design 16(1999), 127-147.[14] A.R.M. Piah, T.N.T. Goodman & K. Unsworth, Positivity-preserving scattered datainterpolation. In Mathematics of Surfaces 2005, LNCS 3604, R. Martin et. al. (eds),Springer-Verlag Berlin Heidelberg (2005) , 336-349.[15] R.J. Renka, Algorithm 833: CSRFPACK-Interpolation of Scattered Data with a C1Convexity-Preserving Surface, ACM Transactions on Mathematical Software 30(2)(2004), 200-211.[16] D.S. Scott, The complexity of interpolating given data in three space with a convexfunction of two variables, J. Approx. Theory 42 (1984), 52-63.

Convexity-preserving scattered data interpolation

Abstract

Similar works

Full text

Available Versions

UUM Repository