# Nonnegativity Preserving Interpolation byC1 Bivariate Rational Spline Surface

Authors: Xingxuan Peng; Zhihong Li; Qian Sun
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

This paper is concerned with the nonnegativity preserving interpolation of data on rectangular grids. The function is a kind of bivariate rational interpolation spline with parameters, which isC1 in the whole interpolation region. Sufficient conditions are derived on coefficients in the rational spline to ensure that the surfaces are always nonnegative if the original data are nonnegative. The gradients at the data sites are modified if necessary to ensure that the nonnegativity conditions are fulfilled. Some numerical examples are illustrated in the end of this paper.


---

## Body

## 1.  Introduction

Interpolation to the scientific data is of great significance in the area of computer-aided geometric design. Particularly, there is often some property inherent in the data which one wishes to preserve when the interpolant is visualized. One useful shape property is nonnegativity: one may have all data values nonnegative and seek an interpolant that is everywhere nonnegative. In this paper, we are concerned with the nonnegativity preserving bivariate interpolation to data on a rectangular grid.Several kinds of surfaces are concerned with nonnegativity preserving interpolation on rectangular grid. For example,C1 biquadratic splines on a refined rectangular grid have been considered in paper [1]. In paper [2, 3], Brodlie et al. followed the same approach (but for bicubic interpolation) of maintaining nonnegativity by modifying estimated slopes at data points. In [4], the interpolant is piecewise an average of two blending surfaces. In [5], C1 interpolating surface is constructed piecewise as a convex combination of two bicubic Bézier patches with the same set of boundary Bézier ordinates. Sufficient nonnegativity conditions on the Bézier ordinates are derived to ensure the nonnegativity of a bicubic Bézier patch. The Bézier ordinates are modified locally to fulfill the sufficient nonnegativity conditions.The surface we consider here is rational spline. Rational spline with parameters has been considered in recent years. Those kinds of interpolation splines have simple and explicit mathematical representation. In paper [6], a bivariate rational interpolation is constructed using both function values and partial derivatives of the function being interpolated as the interpolation data. The convexity control method and point control method have been studied in [7, 8]. In paper [9], a rational cubic function is extended to a rational bicubic partially blended functions (Coons-patches) and the constraints on parameters are derived to visualize the shape of nonnegative surface data. In this paper, we consider the rational spline from different points of view. The sufficient nonnegativity conditions of the bicubic Bézier patch are introduced to consider this problem. We get the sufficient nonnegativity conditions of the rational spline. Gradients at data sites are modified if necessary to ensure that the nonnegativity conditions are fulfilled. It is designed in such a way that no additional points need to be supplied and the sufficient conditions are simple and explicit.The paper is organized as follows. In Section2, the bivariate rational spline is introduced. Section 3 deals with the smoothness of the interpolating surface. In Section 4, the sufficient nonnegativity conditions of the rational spline are derived and a local scheme for C1 nonnegativity preserving interpolation is described. We conclude in Section 5 by illustrating the method with some graphical examples.


## 2.  Rational Interpolation Spline

In this section, the univariate rational cubic spline is introduced which was developed by Hussain and Sarfraz [10]. We extend it to a bivariate rational interpolation spline function.LetΩ:[a,b;c,d] be a planar region, {(xi,yj,fi,j,∂fi,j/∂x,∂fi,j/∂y),(i=1,2,…,n;j=1,2,…,m)}  a given set of data points, where a=x1<x2<···<xn=b and c=y1<y2<···<ym=d are the knot sequences. And let fi,j,∂fi,j/∂x,∂fi,j/∂y represent f(xi,yj),∂f(xi,yj)/∂x,∂f(xi,yj)/∂y, respectively. ∂f(xi,yj)/∂x and ∂f(xi,yj)/∂y are not always given; they can be estimated by the method in [11].For any point(x,y)∈[xi,xi+1;yj,yj+1] in the (x,y)-plane, we construct the x-direction interpolating curve Si,j*(x) in [xi,xi+1] for each y=yj,j=1,2,…,m:(2.1)Si,j*(x)=pi,j*(x)qi,j*(x),i=1,2,…,n-1,
where(2.2)pi,j*(x)=(1-θ)3fi,j+αi,jθ(1-θ)2Vi,j*+βi,jθ2(1-θ)Wi,j*+θ3fi+1,j,qi,j*(x)=(1-θ)3+αi,jθ(1-θ)2+βi,jθ2(1-θ)+θ3,Vi,j*=fi,j+hiαi,j∂fi,j∂x,Wi,j*=fi+1,j-hiβi,j∂fi+1,j∂x,hi=xi+1-xi,θ=x-xihi,
with free parameters αi,j>0,βi,j>0.Obviously, the rational cubic interpolation is unique for the given data(xr,fr,j,∂fr,j/∂x,r=i,i+1) and the parameters αi,j,βi,j. It has the following properties:(2.3)Si,j*(xi)=fi,j,Si,j*(xi+1)=fi+1,j,Si,j′*(xi)=∂fi,j∂x,Si,j′*(xi+1)=∂fi+1,j∂x.For each pair(i,j),i=1,2,…,n-1 and j=1,2,…,m-1, using the x-direction interpolation function Si,j*(x),i=1,2,…,n-1;j=1,2,…,m, we can define the bivariate rational interpolating function Si,j on [xi,xi+1;yj,yj+1] as follows:(2.4)Si,j(x,y)=pi,j(x,y)qi,j(y),i=1,2,…,n-1,j=1,2,…,m-1,
where(2.5)pi,j(x,y)=(1-η)3Si,j*(x)+μi,jη(1-η)2Vi,j+ωi,jη2(1-η)Wi,j+η3Si,j+1*(x),qi,j(y)=(1-η)3+μi,jη(1-η)2+ωi,jη2(1-η)+η3,(2.6)Vi,j(x,y)=Si,j*(x)+ljμi,jdi,j(x,yj),Wi,j(x,y)=Si,j+1*(x)-ljωi,jdi,j+1(x,yj+1),di,s(x,ys)=[(1-θ)3+θ(1-θ)2](∂fi,s/∂y)+[θ2(1-θ)+θ3](∂fi+1,s/∂y)(1-θ)3+αi,jθ(1-θ)2+βi,jθ2(1-θ)+θ3,θ∈[0,1],s=j,j+1.μi,j>0 and ωi,j>0 are free parameters, and lj=yj+1-yj,η=(y-yj)/lj.It is obvious thatdi,s(x,ys) satisfies di,s(xr,ys)=∂fr,s/∂y,r=i,i+1,s=j,j+1. And the bivariate rational function Si,j(x,y) satisfies the interpolation conditions Si,j(xr,ys)=f(xr,ys),∂Si,j(xr,ys)/∂x=∂fr,s/∂x,∂Si,j(xr,yj)/∂y=∂fr,s/∂y,r=i,i+1 and s=j,j+1.


## 3.  The Smoothing Conditions

In this section, the smoothing conditions of the rational splineSi,j(x,y) defined by (2.4) are derived. The rational interpolating function Si,j*(x) defined by (2.1) is a piecewise Hermite interpolant, and it has continuous first-order derivative when x∈[x1,xn]. So the bivariate interpolating function Si,j defined by (2.4) has continuous first-order partial derivatives ∂Si,j(x,y)/∂x and ∂Si,j(x,y)/∂y in the interpolating region [x1,xn;y1,ym] except ∂Si,j(x,y)/∂x at the points (xi,y),i=1,2,…,n-1 for every y∈[yj,yj+1],j=1,2,…,m-1. So it is sufficient for Si,j(x,y)∈C1 in the whole interpolating region [x1,xn;y1,ym] if ∂Si,j(xi+,y)/∂x=∂Si,j(xi-,y)/∂x holds. This leads to the following theorem.Theorem 3.1.
If the knots are equally spaced for variablex; that is,  hi=(b-a)/n,  a sufficient condition for the interpolating function Si,j(x,y),i=1,2,…,n-1;j=1,2,…,m-1, to be C1 in the whole interpolating region [x1,xm;y1,yn] is that the parameters μi,j= constant, ωi,j= constant, and αi,j+βi-1,j=2  for each j∈[1,2,…,m-1] and all i=1,2,…,n-1.Proof.
From the analysis above, without loss of generality, for any pair(i,j),1≤i≤n-1,1≤j≤m-1, and y∈[yj,yj+1], it is sufficient to prove
(3.1)∂Si,j(xi+,y)∂x=∂Si,j(xi-,y)∂x.
Since
(3.2)∂Si,j(x,y)∂x=1qi,j(y)[(1-η)3Si,j′*(x)+μi,jη(1-η)2Vi,j′(x)+ωi,jη2(1-η)Wi,j′(x)+η3Si,j+1′*(x)],
from (2.5), we get
(3.3)∂Si,j(x,y)∂x=1qi,j(y){∂fi,j∂x[(1-η)3+μi,jη(1-η)2]Si,j′*(x)+[ωi,jη2(1-η)+η3]Si,j+1′*(x)+η(1-η)2ljdi,j′*(x,yj)-ljη2(1-η)di,j+1′*(x,yj+1)∂fi,j∂x}.
And since
(3.4)Si,j′*(xi+)=∂fi,j∂x,Si-1,j′*(xi-)=∂fi,j∂x,di,j′*(xi+,yj)=1-αi,jhi∂fi,j∂y,di,j′*(xi-,yj)=βi,j-1hi∂fi+1,j∂y,
we have
(3.5)∂Si,j(xi,y)∂x|x=xi+=1qi,j(y)[(1-η)3+μi,jη(1-η)2]∂fi,j∂x+[ωi,jη2(1-η)+η3]∂fi,j+1∂x+η(1-η)2lj1-αi,jhi∂fi,j∂y-η2(1-η)lj1-αi,j+1hi∂fi,j+1∂y,
where
(3.6)qi,j(y)=(1-η)3+μi,jη(1-η)2+ωi,jη2(1-η)+η3.
Similarly we get
(3.7)∂Si-1,j(x,y)∂x|x=xi-=1qi-1,j[(1-η)3+μi-1,jη(1-η)2]∂fi,j∂x+[ωi-1,jη2(1-η)+η3]∂fi,j+1∂x+η(1-η)2ljβi-1,j-1hi-1∂fi,j∂y-η2(1-η)ljβi-1,j+1-1hi-1∂fi,j+1∂y,
where
(3.8)qi-1,j(x,y)=(1-η)3+μi-1,jη(1-η)+ωi-1,jη2(1-η)+η3.
If (3.1) holds, it needs (3.5) = (3.7); it can be seen that μi,j=μi-1,j,ωi,j=ωi-1,j,αi,j+βi-1,j=2, and hi=hi-1.
This completes the proof.The interpolating scheme above begins inx-direction first. If the interpolation begins with y-direction first, we would get a restriction on the data in the y-direction.


## 4. Construction of Nonnegativity Preserving Interpolating Surface

In this section we first introduce the following theorem in paper [5], which is the basis of our discussion below.Theorem 4.1.
LetP(u,v)=∑i=03∑j=03bi,jBi3(u)Bj3(v),u,v∈[0,1], where {b0,0,b3,0,b0,3,b3,3}={αι,βι,γι,ι}, with ι>0 and α≥β≥γ≥1. Let λ=γ if ι and γι are values at the diagonal vertices; otherwise λ=β. If b1,0,b2,0,b1,3,b2,3,b0,1,b0,2,b3,1,b3,2,b1,1,b2,1,b1,2,b2,2≥-l/3a, where a is the smallest solution in (1,5] of
(4.1)-27λ2a4+108λ2a3+(288λ-162λ2)a2+(108λ2-320λ+256)a-27λ2+32λ=0,
then p(u,v)≥0, for all u,v∈[0,1].Theorem4.1 gives us the sufficient nonnegativity conditions on the Bézier ordinates, which ensure the nonnegativity of a bicubic Bézier patch. Now we consider the nonnegativity condition for the rational spline Si,j(x,y) defined by (2.4). For the parameters μi,j≥0,ωi,j≥0,qi,j(x,y) defined in (2.5) is positive obviously. Now we consider pi,j(x,y) defined in (2.5). Assume αi,j and βi,j are constant for each i∈1,2,…,n-1 and all j=1,2,…m-1. Function pi,j(x,y) can be expressed as follows:(4.2)pi,j(x,y)=1qi,j*{(1-θ)3[(1-η)3fi,j+η(1-η)2(μi,jfi,j+lj∂fi,j∂y)+η2(1-η)(ωi,jfi,j+1-lj∂fi,j+1∂y)+η3fi,j+1]+θ(1-θ)2[(1-η)3αi,jVi,j*+η(1-η)2(μi,jαi,jVi,j*+lj∂fi,j∂y)+η2(1-η)(ωi,jαi,jVi,j+1*-lj∂fi,j+1∂y)+η3αi,jVi,j+1*]+θ2(1-θ)[(1-η)3βi,jWi,j*+η(1-η)2(μi,jβi,jWi,j*+lj∂fi+1,j∂y)+η2(1-η)(ωi,jβi,jWi,j+1*-lj∂fi+1,j+1∂y)+η3βi,jWi,j+1*]+θ3[(1-η)3fi+1,j+η(1-η)2(μi,jfi+1,j+lj∂fi+1,j∂y)+η2(1-η)(ωi,jfi+1,j+1-lj∂fi+1,j+1∂y)+η3fi+1,j+1]}=ti,jqi,j*,
where(4.3)qi,j*=(1-θ)3+αi,jθ(1-θ)2+βi,jθ2(1-θ)+θ3.It can be seen thatqi,j* is positive for αi,j>0 and βi,j>0. The function ti,j is a bicubic Bézier patch. It is nonnegative if the Bézier ordinates of it satisfy the conditions in Theorem 4.1.The Bézier ordinatesbp,q(p,q=0,1,2,3) of ti,j are(4.4)b0,0=fi,j,b0,1=μi,jfi,j+lj∂fi,j∂y,b0,2=ωi,jfi,j+1-lj∂fi,j+1∂y,b0,3=fi,j+1,b1,0=αi,jfi,j+hi∂fi,j∂x,b1,1=μi,jαi,jfi,j+μi,jhi∂fi,j∂x+lj∂fi,j∂y,b1,2=ωi,jαi,jfi,j+1+ωi,jhi∂fi,j+1∂x-lj∂fi,j+1∂y,b1,3=αi,jfi,j+1+hi∂fi,j+1∂x,b2,0=βi,jfi+1,j-hi∂fi+1,j∂x,b2,1=μi,jβi,jfi+1,j-μi,jhi∂fi+1,j∂x+lj∂fi+1,j∂y,b2,2=ωi,jβi,jfi+1,j+1-ωi,jhi∂fi+1,j+1∂x-lj∂fi+1,j+1∂y,b2,3=βi,jfi+1,j+1-hi∂fi+1,j+1∂x,b3,0=fi+1,j,b3,1=μi,jfi+1,j+lj∂fi+1,j∂y,b3,2=ωi,jfi+1,j+1-lj∂fi+1,j+1∂y,b3,3=fi+1,j+1.Bézier ordinates ofti,j need not to satisfy the nonnegative conditions. To ensure this, we shall impose upon these Bézier ordinates the conditions bp,q≥-ι/3a(p,q=0,1,2,3) according to Theorem 4.1, where ι=min{fi,j,fi,j+1,fi+1,j,fi+1,j+1} and a is determined by (4.1). This can be achieved by modifying if necessary the gradients at vertices Vi,j=(xi,yj). The modification of derivatives ∂fi,j/∂x and ∂fi,j/∂y at a vertex is performed by scaling each of them with a positive factor α< 1. The scaling factor α is obtained by taking into account all the rectangular patches sharing that vertex. For the four rectangles which share vertex Vi,j, denote vertex Vi,j as O and its adjacent vertices as A,B,C,D,E,F,G,H, respectively (see Figure 1).Figure 1
VertexO with its associated rectangles.Consider rectangle 1 and the lower bound-ι1/3a1 where ι1=min{S(O),S(A),S(B),S(C)} and a1 is obtained by solving (4.1) in Theorem 4.1. Scalar αOA1 is defined as follows. If  μi,jfi,j+lj(∂fi,j/∂y)≥-ι1/3a1, then αOA1=1.  Otherwise, αOA1 is defined by the equation μi,jfi,j+αOA1lj(∂fi,j/∂y)=-ι1/3a1. Similarly for the scalar αOC1, if αi,jfi,j+hi(∂fi,j/∂x)≥-ι1/3a1, then αOC1=1; otherwise αOC1 is given by the equation αi,jfi,j+αOC1hi(∂fi,j/∂x)=-ι1/3a1. For the scalar αOAC1, if μi,jαi,jfi,j+μi,jhi(∂fi,j/∂x)+lj(∂fi,j/∂y)≥-ι1/3a1, then αOAC1=1; otherwise αOAC1 is determined by the equation μi,jαi,jfi,j+αOAC1(μi,jhi(∂fi,j/∂x)+lj(∂fi,j/∂y))=-ι1/3a1. Then define α11=min{αOA1,αOAC1},α12=min{αOC1,αOAC1}. By using the same argument above, we get α21 and α22,α31 and α32,α41 and α42 for rectangles 2, 3, 4, respectively. In order for all the Bézier ordinates adjacent to O to fulfill the positivity preserving conditions stated in Theorem 4.1, set αO1=min{α11,α21,α31,α41}, αO2=min{α12,α22,α32,α42}.IfαO1<1 (or αO2<1), the x (or y) partial derivatives at O are redefined as αO1(∂fO/∂x) and αO2(∂fO/∂y). Then the Bézier ordinates adjacent to O in each rectangle are redetermined. The above process is repeated at all the nodes Vi,j.For the boundary nodeO, which belongs to one or two rectangles of the rectangular grid, αO is defined in a similar way. The only difference is that we consider only one or two rectangles instead of four.Now all the Bézier ordinates are determined, soti,j is nonnegative. Thus we can get the C1 rational interpolant S, which is nonnegative.


## 5.  Numerical Examples

We shall illustrate our discussion with the following examples.Example 5.1.
In the first example the data are given as follows:(5.1)f1,1=0.1,f1,2=1.5,f2,1=2,f2,2=2.5,∂f1,1∂x=-3,∂f1,2∂x=0.5,∂f2,1∂x=-0.1,∂f2,2∂x=-0.1,∂f1,1∂y=-0.1,∂f1,2∂y=0.01,∂f2,1∂y=-0.02,∂f2,2∂y=-0.01.
Figure2 shows the rational interpolation surface, which loses the nonnegativity in its display. Figure 3 is a different view of Figure 2 after making a rotation. It confirms that the surface is not preserving nonnegativity feature. Figures 4 and 5 show the surface after modifying by the scheme of this paper, which are indeed nonnegative. In fact, they are changed as follows:
(5.2)∂f1,1∂x=-0.2876,∂f1,2∂x=0.5,∂f2,1∂x=-0.1,∂f2,2∂x=-0.1,∂f1,1∂y=-0.0496,∂f1,2∂y=0.01,∂f2,1∂y=-0.02,∂f2,2∂y=-0.01.Figure 2
The unconstrained interpolating surface to data.Figure 3
A different view of the surface, of Figure2, after rotation.Figure 4
Nonnegativity-preserving interpolating surface to data.Figure 5
A different view of the surface, of Figure4, after rotation.Example 5.2.
In the second example, data points are generated from functiong(x,y) [12]:
(5.3)g(x,y)={2(y-x),0≤y-x≤0.5,1,y-x≥0.5,0.5cos⁡(4π(x-1.5)2+(y-0.5)2)+0.5,(x-1.5)2+(y-0.5)2≤116,0,elsewhere on[0,2]×[0,1].
Figure6 shows the unconstrained interpolating surface, which loses the nonnegativity in its display. Output from the nonnegativity preserving scheme of this paper is shown in Figure 7. It clearly shows that the surface remains nonnegative everywhere and visually pleasant.Figure 6
The unconstrained interpolating surface to data fromg.Figure 7
Nonnegativity-preserving interpolating surface to data fromg.


## 6. Conclusion

In this paper, we construct a nonnegativity preserving interpolant of data on rectangular grids by using the rational spline. The sufficient nonnegativity conditions are derived, which are simple and explicit. And the experimental results illustrate the feasibility and validity of our method.



---

*Source: 624978-2012-04-10.xml*