Range restricted positivity-preserving scattered data interpolation

The construction of a range restricted bivariate C 1 ( or G 1 ) interpolant to scattered data is considered in which the interpolant is positive everywhere if the original data are positive. This study is motivated by earlier work in which sufficient conditions are derived on Bézier points in order to ensure that surfaces comprising cubic Bézier triangular patches are always positive and satisfy C 1 continuity conditions. In the current work, simpler and more relaxed conditions are derived on the Bézier points. The gradients at the data sites are then calculated (and modified if necessary) to ensure that these conditions are satisfied. Each triangular patch of the interpolating surface is formed as a convex combination of three quartic Bézier triangular patches. Its construction is local and easily extended to include as upper and lower constraints to the interpolant surfaces of the form z = P(x,y) where P is a polynomial of degree less or equal to 4. Moreover, C 1 ( or G 1 ) piecewise polynomial surfaces consisting of polynomial pieces of the form z = P(x,y) on the triangulation of the data sites are also admissible constraints. A number of examples are presented.


Introduction
The properties that are most often used to quantify "shape" in shape preserving interpolation are positivity, convexity and monotonicity.The problem of positivity preserving interpolation, i.e. interpolation to positive data by a positive function, is often of interest.This problem could arise if one has data points on one side of a plane, and wishes to have an interpolating surface which is also on the same side of this plane.For instance, when the data from physical experiment are measured in the form of concentration or pressure or meteorology data such as amount of rainfall where negative values are meaningless, it is important for the interpolant to preserve positivity.This paper will propose sufficient conditions on the bivariate quartic function upon triangulation of the data in order to visualize the positive scattered data which may come from certain scientific phenomena.Various methods concerning visualization of positive data or range restricted interpolation using bivariate functions can be found in [1], [2], [8], [11], [12] and [13].
This study is motivated by previous works in [2], [13] and [15].[2] describes the construction of range restricted bivariate C 1 interpolants to scattered data where sufficient non-negativity condition on the Bézier ordinates are derived to ensure the non-negativity of a cubic Bézier triangular patch.
In [2] and [13], a C 1 non-parametric surface is constructed comprising of cubic Bézier triangular patches.Each triangular patch of the interpolating surface is formed as a convex combination of three cubic Bézier triangular patches.An initial value of inner Bezier ordinates in each triangle are computed using the cubic precision method.
In this paper, we will construct a range restricted positivity-preserving bivariate C 1 (or G 1 ) quartic interpolants to scattered data using similar approach adopted in [13].However, each triangular patch of the interpolating surface is formed as a single quartic Bézier triangular patch.Initial values of Bezier ordinates except for the values determined by gradients at vertices of each triangle are computed using the method of minimized sum of squares of principles curvartures ( [7], [9]) with respect to the G 1 continuity conditions [5] on each non-boundary edge over the triangular mesh using the quadratic form of an objective function [15].

Sufficient Positivity Conditions for a Quartic Bézier Triangular Patch
Consider a triangle T, with vertices V 1 , V 2 , V 3 , and barycentric coordinates u,v, w such that any point V on the triangle can be expressed as V = uV 1 + vV 2 + wV 3 , where u + v + w = 1 and u,v,w ≥ 0. A quartic Bézier triangular patch P on T is defined as, P(u,v,w) = u 4 b 400 + v 4 b 040 + w 4 b 004 + 4u 3 vb 310 + 4u 3 wb 301 + 4v 3   (1) where b ijk are the Bézier ordinates of P for i, j, k = 0, 1, 2, 3, 4 and i + j + k = 4 as shown in Figure 1.
We assume that the Bézier ordinates at vertices are strictly positive, i.e. b 400 , b 040 , b 004 > 0. The sufficient conditions on the remaining Bézier ordinates will be derived to ensure the positivity of the entire patch.Let A = b 400 , B = b 040 , and C = b 004 where A, B, C > 0. Our approach is to find the lower bounds of the remaining Bézier ordinates, so that P(u,v,w) = 0. We also assume that, apart from the vertices, the remaining Bézier ordinates have the same value -r (where r > 0).Thus, (1) can be written as, P(u,v,w) = Au 4 + Bv 4 + Cw 4 -r(4u 3 v+ 4u 3 w +4v 3 u +4 v 3 w + 4w 3 u +4w 3  ( From (2), clearly when r = 0, P(u, v, w) > 0. We are interested to find the value of r = r 0 when the minimum value of P(u, v, w) is equal to zero.The derivatives of P in (2) with respect to u, v and w are given by, We know that at the minimum value of P(u, v, w), . Thus, we have Since u + v + w = 1, we obtain . Substituting these u, v and w into (2), we obtain the minimum value of We need to choose a value r = r 0 so that this minimum value of P is zero.From ( 5) and r > 0, we know that, P(u, v, w) = 0 when If A, B and C are strictly positive, then s = s 0 = 1/r 0 is the solution of G(s) = 1.
Now, let us describe the method to determine the value of s 0 for each triangular patch.Since A, B, C > 0, it is easy to show that for s ≥ 0, G' (s) < 0 and G''(s) > 0. Let M = max(A, B, C) and N = min (A, B, C) and it can be verified that

G(s)
To obtain the value of s 0 for given values of A, B and C, we need to find the root of G(s) in ( 7) and G(s) = 1 that will give us a lower bound on the remaining Bézier ordinates, i.e. r 0 = 1/s 0 using simple iterative scheme.In choosing these scheme, we must ensure of one sided convergence of the root.This can be achieved by the method of false-position [3] with an initial estimate for the root will be the value of s for which the line joining 26/N and 26/M has the value 1.The following proposition prescribes the lower bounds for the Bézier ordinates to ensure the positivity of Bézier patch.≥ -r 0 = -1/s 0 , (r, s, t) ≠ (4,0,0), (0,4,0) and (0,0,4), where s 0 is the unique solution of ( 7) and G(s

Proposition
Note that, if any of the values of A, B or C are zero (i.e. the given data are not strictly positive), we will assign the value of r 0 to be zero for that triangle.

Construction of Positivity-Preserving Interpolating Surface
Having established sufficient conditions on the Bézier points as described in section 2, we are now able to construct the interpolating surface.Given a positive scattered data (x i, y i ,z i ), z i ≥ 0, i = 1,2,…,N, we want to construct a C 1 (or G 1 ) positivity-preserving surface z = F(x,y) that interpolate the given data.The surface comprises of quartic Bézier triangular patches each of which is guaranteed to remain positive.The construction process consists of the following steps: 1. Triangulation of the domain data, (x i, y i ), i=1,2,…,N.

2.
Initialise partial derivatives at each data point, then modify them (if necessary) to be consistent with the positivity constraints imposed on the Bezier ordinates.3. Generate the final surface.
We use Delaunay triangulation [4] to triangulate the convex hull of the data points and the estimation of the first order partial derivative of F with respect to x and y is obtained using the method proposed in [6].Let V i , i = 1, 2, 3 be the vertices of a triangle, such that F(V i ) = z i , and the first partial derivatives F x (V i ) and F y (V i ).Then for each triangular patch P as given by (1), the derivative along the triangle edge e jk joining (x j , y j ) to (x k , y k ) is given by De jk P( V j ) = (x k -x j ) F x (V j ) + (y k -y j ) F y (V j ).From the given data, together with estimated derivatives at all (x i , y . By symmetry, we can obtain the remaining 6 control points.
However, the initial estimate of the above edge ordinates may not satisfy the positivity conditions for P. In If it does not, the magnitudes of F x , F y at the vertices need to be reduced so that the conditions are satisfied.The modification of these partial derivatives is achieved by multiplying each derivative at that vertex by a scaling factor 0 < α < 1.The smallest value of α is obtained by considering all triangles that meet at vertex V that will guarantee satisfaction of the positivity condition for all these triangles.For example where subscript j represents quantities corresponding to triangle j.Having adjusted these derivatives, if necessary, the edge Bézier ordinates are recalculated using the above formulae.which need to be done in order to guarantee preservation of positivity and to ensure C 1 (or G 1 ) continuity across patch boundaries.We adopt a similar approach presented in [13] which guarantee C 1 continuity and has minimized sum of squares of principle curvartures (for G 1 ).
Two patches with a common boundary curve satisfy G 1 continuity if both have continuously varying tangent plane along the common curve.Figure 3 shows an example of Bézier control points of two adjacent quartic Bézier triangular patches (denoted by R and S respectively), where h 0 and h 4 are Bézier points of the vertices, g 0 , f 0 , h 1 , g 3 , h 3 and f 3 are obtained from the patch gradients while g 1 , f 1 , h 2 , g 2 and f 2 are points to be determined.We only have to consider {h i , i = 0, 1,…4} as the common boundary curve and {g i , f i , i = 0, 1,…3} which consist of the control points in each patch.Details of derivation with regard to the C 1 (or G 1 ) conditions can be found in [5].

C 1 continuity conditions
The necessary and sufficient conditions for C 1 continuity along the common boundary curve as shown in Figure 3 are given by [ 5] as follows: where W 1 = uV 1 + vV 2 + wV 3 and u + v + w = 1.
Note that, ( 8) and (11) will be automatically satisfied since the Bézier ordinates g 0 , g 3 , f 0 , f 3 , h 1 , h 3 are determined by the gradients at the vertices V 2 and V 3 respectively.

G 1 continuity conditions
The conditions satisfying G 1 continuity between the two adjacent patches as shown in Figure 3 can be written as (see [5] for further details) where α and β are constants.Since the values of g 0 , f 0 , h 0 , h 1 , g 3 , h 3 , h 4 and f 3 are already known, α and β for each triangle can thus be determined from ( 12) and (15).
We can also obtain expressions ( 9) and ( 10), ( 13) and ( 15) on all non-boundary edges over the whole triangular mesh and be represented as where A and C are l x n (l < n) coefficient matrices, x is a n x 1 unknown vector consisting of all remaining Bézier ordinates to be determined for the entire triangular mesh, and b and d are l x 1 constant vectors.
Note that ( 16) and ( 17) are underdetermined linear systems that will always have a solution since the rows of A and C are linearly independent.We can obtain an approximation to the unknown vector x by solving the linear system as in ( 16) and ( 17) using the least square method, but we might have a poor approximation of the coefficients and the resulting surface might have over-flatness or undesired undulation.To overcome this problem, we shall follow a similar approach as in [15] using minimized sum of squares of principle curvatures with respect to the constraint (17).Our aim is to find the function F(u,v) which will minimize the functional , where m is the number of triangles in a mesh with u + v + w = 1 or in the form of matrix-vector representation where M is a real (n x n) symmetric matrix, e is a (1x n) row vector, x is a (n x 1) column vector representing the unknown Bézier points for the entire triangular mesh and c 1 as a real constant.
In order to find F(u, v) which will minimize I(F(u, v)) leads us to an optimisation problem, x T Mx + ex + c subject to the G 1 continuity constraints Cx = d (refer to [15] for further details of this method).
An initial estimate of the edge ordinates b 220 , b 202 , b 022 and inner Bézier ordinates b 211 , b 121 , b 112 in each triangle obtained by the above method may not satisfy the positivity conditions for P(u, v, w) as stated in Proposition.We may need to adjust the above coefficients in order to fulfill Proposition and the C 1 continuity (( 9) and ( 10)) or G 1 continuity (( 13) and ( 14)).Thus, for every common edge of the two patches R and S as in Figure 3, we may need to adjust the coefficients of g 1 , f 1 , g 2 , f 2 and h 2 .We have adopted similar approach previously done for C 1 triangular cubic patches in [2] but with slight modification for the triangular quartic patches because more coefficients are involved.We also note that, this modification is local.For each patch P, when all the Bézier ordinates have been assigned, the final positive interpolating surface can be generated using (1) for the barycentric coordinates u, v and w.

Range Restricted Interpolation
In the previous section, we have described the construction of C 1 (or G 1 ) interpolating surface which is constrained to lie above the plane z = 0. We shall extend our scheme to include a larger set of constraint surfaces that are of the form z = C(x, y) where C(x, y) is a constant, linear, quadratic, cubic or quartic polynomial, i.e.C(x, y) = ax 4 + bx 3 y + cx 2 y 2 + dy 4 + exy 3 + fx 3 + gx 2 y + hy 3 + ixy 2 + jx 2 + kxy + ly 2 + mx + ny + o where a, b, c, d, e, f, g, h, j, k, l, m, n and o are real numbers.These surfaces are considered because they can be expressed as quartic Bézier triangular patches on each triangle of the triangular mesh.Thus, C 1 (or G 1 ) piecewise polynomial surfaces consisting of polynomial pieces of the form z = C(x , y) on the triangulation of data sites can also be treated as constraint surfaces.
We would like to generate a C 1 (or G 1 ) interpolating surface z = F(x, y) through the data points (x i , y i , z i ) , i = 1, 2, . .., N which lies either above or below the constraint surface or lie between both the constraint surfaces.This problem can easily be reduced to the positivity preserving interpolation case which we have considered earlier.Assume the data points lie above the constraint surface.The initial problem of constructing the interpolation surface F(x, y) with respect to the constraint surface C(x, y) is similar to the construction of a function G(x, y) = F(x, y) -C(x, y) such that G is positive and C 1 (or G 1 ) with G(x i , y i ) = * i z where * i z = z i -C(x i , y i ) is a new set of data points.The initial values of gradients of G at the data sites are obtained from G x (x i , y i ) = F x (x i , y i )-C x (x i , y i ) and G y (x i , y i ) = F y (x i , y i ) -C y (x i , y i ).The gradients of G are modified if necessary using the method described in the previous section.Thus, the positivity-preserving interpolating surface F(x, y) is constructed piecewise as a single quartic triangular patch, where G(x, y) is also a single piecewise quartic triangular patch.We can use a similar construction method if the data points lie below the constraint surface by writing G(x, y) as C(x, y) -F(x, y).The above construction method can also be extended to describe the interpolating surface that lie between both the upper and lower constraint surfaces (see [2] for further details).

Examples
In this section, we will illustrate our interpolating scheme using two test functions.The first example from the well-known data set taken from [10] comprises 36 data points of For this example, we shall use z = 1.001 as the upper constraint and z = -0.001as the lower constraint similar to [2].The result of unconstrained C 1 and G 1 interpolating surfaces are shown in Figure 4, which clearly show unwanted oscillations at the upper and bottom shelf regions.Figures 5 and 6 show the result when the positivity conditions with either the upper or lower constraints are imposed and unwanted oscillations still occurs at the upper or bottom shelf region.When we imposed positivity conditions together with both the upper and lower constraints, the interpolating surface does not oscillate unnecessarily and stays in-between the two constrained surfaces as shown in Figure 7.The second example is the function g taken from [2] where In this example, we use 33 data points which g interpolates.A linear function C1(x,y) = 0.1xy + 0.825x -0.215 and a cubic polynomial C2(x,y) = -3x 3 + 5.

Summary and Conclusion
In this study, we have considered the generation of non-parametric surfaces that interpolate positive scattered data.We have imposed relaxed and simpler conditions on Bézier ordinates by modifying the previous work in [13] and [15].We also extend the problem of C 1 (or G 1 ) positivity preserving interpolants to the range restricted scattered data cases similar to those considered for C 1 interpolants described in [2].This scheme (and that of [2] and [13]), however, suffer from the constraint that in order to prevent a surface patch becoming negative all Bézier points within the corresponding triangle (except those at the vertices) may be assigned the same value according to Proposition.This makes the problem tractable and achieves the desired outcomes.More flexibility would be achieved if they could be adjusted independently while still ensuring positivity.Preliminary investigations indicate that this approach leads to an optimisation problem, however the overall solution procedure is significantly more complex than that described in this paper.A study of this approach is ongoing.

Figure 1 :
Figure 1: Relative locations of Bezier ordinates for P(u, v, w)

Figure 2 Figure 2 :
Figure 2: Function G(s) for s ≥ 0 Consider the quartic Bézier triangular patch P(u, v, w) with b 400 = A, b 040 = B, b 004 = C, where A, B, C > 0. If b rst i ) we can now determine all b rst except of b 220, b 202 , b 022 , b 211 , b 121 and b 112 .For example, we have:
Next, we shall calculate the remaining edge ordinates b 220 , b 202 , b 022 and inner Bézier ordinates b 211 , b 121 , b 112

Figure 3 :
Figure 3: Control points of adjacent quartic Bézier triangular patches

Figure 4 :Figure 5 :Figure 6 :Figure 7 :
Figure 4: (a) Unconstrained C 1 interpolating surface (b) Unconstrained G 1 interpolating surface 55x 2 + 0.2xy -2.25x + 0.207 are used as lower constraint surfaces respectively.The unconstrained interpolating surfaces are shown in Figures 8 and 10, while Figures 9 and 11 show the constrained interpolants.Figures 8 and 10 show part of the unconstrained interpolating surfaces cross the two constraint surfaces at two different regions but the surface interpolant stays above the two lower constrained surfaces given in Figures 9 and 11 .