Computer-aided geometric design (CAGD) provides the mathematical basis when dealing with geometric datasets. The term “CAGD” was coined by Barnhill and Riesenfeld in 1974 [1–5]. One of the main tools in CAGD is a spline function. The spline function was introduced by Schoenberg in 1967 for application in statistics. Since then, the spline has been used extensively in shape designing especially in car industries, in the design of airplane fuselages and wings at Boeing Company, as well as in the development of cartoon movies. General Motors developed its first computer-aided design/computer-aided manufacturing (CAD/CAM) system called DAC-I by utilizing the spline curve and surface techniques initiated by de Boor and Gordon [5]. Furthermore, in some applications, such as data visualization, there is a need to reconstruct the curve or surface for the data collected from the measurement or experimental laboratory. This task can be achieved by using interpolation or approximation schemes.

Shape-preserving interpolation is an important topic in CAGD that can be used to preserve the geometric shape of the dataset such as positivity, convexity, and monotonicity. For instance, given positive data, then the final interpolating curve or surface must be positive everywhere. Any negativity is meaningless and unacceptable. Some examples of positive data are the amount of gas discharge during the experiment (Brodlie and Butt [6], Butt and Brodlie [7], Brodlie et al. [8], Sarfraz et al. [9]). Meanwhile, monotonicity-preserving concerns about preserving the monotone datasets, i.e., either monotonically increasing or monotonically decreasing. For instance, the approximation of couples and quasi-couples in statistics and the empirical option model in finance are always involving monotone datasets [10]. The blood uric acid in patients suffering from gout is also an example of monotone data. Meanwhile, convexity always arises in finance and engineering-based problems such as the optimal control, non-linear programming, and parameter estimation problems.

In general, the standard cubic Hermite polynomial and cubic spline interpolation functions are unable to produce the interpolating curve or surface with shape-preserving properties. Therefore, to cater this weakness, many researchers have proposed several methods for shape-preserving approximation and interpolation. One of the simplest ways is to use rational splines in many different forms such as cubic/cubic, cubic/quadratic, cubic/linear, and quartic/quadratic. Those schemes provided more design parameters which can be used to modify the interpolating curves or surfaces. In the next paragraphs, we provide some literature review on positivity-preserving interpolation schemes.

The modified quadratic Shepard (MQS) method has been used widely for visualization of scattered data. For example, Asim et al. [11, 12] improved the original MQS method for positivity-preserving on scattered datasets. Their strategy is to reduce the deviation from the original shape while at the same time preserving positivity in the sense of least square. Brodlie et al. [8] constructed the MQS method to interpolate scattered data of any dimensionality. Their scheme preserves the positivity of the data for curves or surfaces by forcing the quadratic basis functions to be positive. They also extended the method to handle modeling other types of constraints such as the lower bound of 0 and upper bound of 1 and generalized the constraint to any arbitrary functions as lower and upper bounds. But from the numerical results, they indicated that their schemes have 10% error bound associated with the data values. Wu et al. [13] discussed about positivity-preserving for curve and surface approximation and interpolation by using compactly supported radial basis functions (CSRBFs). To preserve the positivity of data, the optimization problem needs to be solved. In [14, 15], the multi-quadric (MQ) function has been used as the quasi-interpolation operator and the CSRBFs are used to construct the interpolation function that preserves the monotonicity and convexity of planar datasets.

Most shape-preserving surface interpolation methods are constructed by extending the work of Casciola and Romani [16]. Nowadays, their schemes are known as partially blended rational bi-cubic functions. Brodlie and Butt [6] constructed the piecewise cubic Hermite interpolant with C 1   continuity to preserve positivity, monotonicity, and convexity, respectively. Their schemes required one or two extra knots to be inserted in the interval where shape violation exists without needs for derivative modification. However, this technique will increase the total number of data points, and indeed, the computation times will be increased.

Brodlie et al. [8] extended the Butt and Brodlie [7] idea to construct the positivity-preserving scheme of positive surface data over a rectangular grid. The sufficient conditions on first partial derivatives and twist values are derived where the obtained values will be projected to the valid interval through the efficient knot insertion algorithm. However, their scheme has increased unnecessary computation time.

Schmidt and Hess [17, 18] derived the necessary and sufficient conditions for the positivity of rational quadratic splines and for the positivity of cubic spline interpolation on intervals. Schmidt [19] described the positive, monotone, and S -convex surface interpolation for data arranged on rectangular grids. The C 1 rational quadratic spline has been used, and the necessary conditions for positivity, monotonicity, and S -convex surface are derived. However, no numerical results were presented. Hussain et al. [20] constructed C 0 shape-preserving surface schemes for 3D positive and convex surface data. The rational bi-quadratic function with eight parameters was used to construct the positive and monotone rational interpolant. But their scheme was unable to produce smooth surfaces as well as no free parameters to alter the final resulting interpolating surfaces. Chan and Ong [21] described a local scheme for range-restricted scattered data interpolation by using cubic triangular Bézier patches. The piecewise interpolating surface is obtained through a convex combination of three cubic Bézier triangular patches. The sufficient conditions for the non-negativity of a cubic Bézier triangle on the Bézier ordinates were derived with respect to a lower bound. The gradients (first-order partial derivatives) at the data sites are modified if necessary, to ensure that the non-negativity conditions are fulfilled. Luo and Peng [22] described the C 1 rational spline as a piecewise rational convex combination of three cubic Bézier triangular patches sharing the same boundary Bézier ordinates. The sufficient conditions for non-negativity were derived on the boundary Bézier ordinates of the adjacent triangle and the normal derivatives at the data points. If in any triangular patch non-negativity was lost, then the gradients at the data points and normal derivatives at the edge knots will be modified to ensure positivity is preserved. Indeed, their main scheme also requires the modification of the first-order partial derivative and the normal derivatives if shape violation is found.

Piah et al. [23] improved the lower bound of Bézier ordinates obtained in the work of Chan and Ong [21] by proposing an alternative scheme which is simpler and has more relaxed conditions of positivity. The new lower bound for Bézier ordinates can become infinitely smaller compared to the lower bound of Bézier ordinates in [21]. Some numerical results were presented to show the capability of the alternative scheme. The main drawback of their scheme is that the actual first-order partial derivative needs to be adjusted if positivity on each triangular patch is violated. Hussain and Hussain [24] constructed positive scattered data interpolation by using a rational cubic spline defined by the side-vertex scheme.

Saaban et al. [25] constructed the positivity-preserving scheme by using quintic triangular Bézier patches. The surface is constructed by using a convex combination comprising of three local quintic triangular Bézier patches. They implemented their proposed method by using rainfall data collected at various meteorology stations in West (Peninsular) Malaysia.

Hussain et al. [26, 27] extended the application of rational quartic spline given by Wang and Tan [28] to preserve positivity, constrained data, and convexity. The data-dependent sufficient conditions for the rational interpolant to satisfy the shape-preserving properties are derived. It was proved by Harim et al. [29] that their sufficient conditions may not successfully produce positive and convex interpolating curves on entire given intervals.

Liu et al. [30] described a new bivariate rational quartic spline (quartic/quadratic) for positivity- and monotonicity-preserving interpolation. The rational bi-quartic spline has four parameters, and the sufficient conditions for positivity and monotonicity are derived on all parameters. Thus, there is no free parameter for shape modification. They have compared the performance of their schemes against that of Hussain and Sarfraz’s [31] for positivity-preserving interpolation. No error measurement is given to show the superiority of the proposed scheme. Han [32, 33] discussed the rational quartic spline with a quadratic denominator. Harim et al. [29, 34, 35] constructed a new rational quartic spline (quartic/quadratic) with three parameters for shape modification. Karim et al. [31, 36, 37] discussed the positivity-preserving interpolation and constrained surface modeling by using rational bi-cubic spline interpolation with 12 parameters. Qin et al. [38] proposed a new method to derive the sufficient condition for the positivity of rational bi-cubic spline interpolation. However, their method requires extra computation time.

In this section, we have given a comprehensive review on the existing schemes of positivity-preserving interpolation. Table 1 summarizes the comparison of various shape-preserving interpolation schemes related to our studies.

From Table 1, we have come out with our main objectives of the study. There are three main objectives given as follows:

1) We extend the univariate rational quartic spline with three parameters from Harim et al. [34] to the bivariate rational quartic spline interpolation. We employ the partially blended scheme initiated by Casciola and Romani [16]. There are three free parameters on each boundary of the rectangular meshes. This will give 12 parameters that can be used to control the shape of the interpolating surface.

Furthermore, we identified several advantages of the proposed partially blended rational bi-quartic spline interpolation for positivity-preserving:

1) There are 12 parameters in the description of the rational bi-quartic spline, and eight are free parameters for shape modification. Meanwhile, Hussain and Ali [37], Hussain and Hussain [40], Tian et al. [41], and Sarfraz et al. [9] schemes have no free parameters for shape modification.

The remainder of this paper is organized as follows. Materials and Methods is devoted to materials and methods used in this study. This includes the method to estimate first derivatives using the arithmetic mean method (AMM) and geometric mean method (GMM), a review on the rational quartic spline, the derivation of the rational bi-quartic spline, and local control and quadratic reproducing properties. Meanwhile, Positivity-Preserving Interpolation discusses the construction of positive rational bi-quartic spline interpolation on each boundary curve network. Furthermore, an efficient algorithm is presented in detail. Results and Discussion is devoted to the results obtained and discussion. Conclusion gives the concluding remarks.

The rational bi-quartic spline defined in Eq. 17 is not positivity-preserving for positive datasets. As seen in Figure 2A, the partially blended rational bi-quartic spline is not positivity-preserving interpolation. We can obtain a positive surface by manipulating all free parameters. However, this approach is not practical and really time-consuming. Therefore, we want to find an automatic method to calculate the parameter values that will produce a positive rational interpolant. To do this, we derive the sufficient condition for the positivity of the rational quartic splines on all four boundary curve networks (see Figure 1). On each boundary, there will be two free parameters. This results in eight free parameters for shape modification. The derivation of the sufficient condition is described as follows.

Let ( x i , y i , F i , j ) be strictly positive data (to avoid dividing by zero at later stages) defined over a rectangular grid [ x i , x i + 1 ] × [ y j , y j + 1 ] , i = 0,1 , ... , n − 1 ; j = 0,1 , ... , m − 1   , such that F i , j > 0 ,       ∀ i , j . The partially blended rational bi-quartic spline in Eq. 17 is positive if all boundary curves defined in Eqs. 18–21 are positive, i.e., S ( x , y j ) > 0 ,       S ( x , y j + 1 ) > 0 ,       S ( x i , y ) > 0 , and S ( x i + 1 , y ) > 0 ,   respectively. The following theorem summarizes the construction of positivity-preserving interpolation using the proposed rational bi-quartic spline.

Theorem 4. Given strictly positive surface data F i , j > 0 ,       ∀ i , j ,   the piecewise partially blended rational bi-quartic function S ( x , y ) in Eq. 17 defined over the rectangular mesh [ x i , x i + 1 ] × [ y j , y j + 1 ] , i = 0,1 , ... , n − 1 ; j = 0,1 , ... , m − 1 , preserves the positivity of the data if the parameters satisfy the following sufficient conditions: α i , j > 0 , α i , j + 1 > 0 , β i , j > 0 , β i , j + 1 > 0 , α ^ i , j > 0 , α ^ i + 1 , j > 0 , β ^ i , j > 0 , β ^ i + 1 , j > 0 γ i , j > Max { 0 , − α i , j ( h i F i , j x + 2 F i , j ) / F i , j , β i , j ( h i F i + 1 , j x − 2 F i + 1 , j ) / F i + 1 , j } γ i , j + 1 > Max { 0 , − α i , j + 1 ( h i F i , j + 1 x + 2 F i , j + 1 ) / F i , j + 1 , β i , j + 1 ( h i F i + 1 , j + 1 x − 2 F i + 1 , j + 1 ) / F i + 1 , j + 1 } γ ⌢ i , j > Max { 0 , − α ^ i , j ( h ^ i F i , j y + 2 F i , j ) / F i , j , β ^ i , j ( h ^ j F i , j + 1 x − 2 F i , j + 1 ) / F i , j + 1 } γ ⌢ i + 1 , j > Max { 0 , − α ^ i + 1 , j ( h ^ i F i + 1 , j y + 2 F i + 1 , j ) / F i + 1 , j , β ^ i + 1 , j ( h ^ j F i + 1 , j + 1 y − 2 F i + 1 , j + 1 ) / F i + 1 , j + 1 } } . (27)

Proof.

To prove the theorem, we adopted the idea from Karim et al. [22]. As conjectured by Casciola and Romani [16], the surface is shape-preserving if all four curve networks are shape-preserving. This can be summarized as that the rational bi-quartic surface S ( x , y )   is positive if and only if all four curves are positive.

From Equations (18)–(21), S ( x , y j ) > 0 , if A i > 0 ,     i = 0,1,2,3. α i , j > 0 , β i , j > 0 , γ i , j > 0 , A 0 > 0 ,   A 2 > 0 ,   and   A 4 > 0.   Therefore, S ( x , y j ) > 0 , if and only if A 1 > 0 and A 3 > 0. Both conditions can be simplified as follows: γ i , j > − α i , j ( h i F i , j x + 2 F i , j ) / F i , j , (28) γ i , j > β i , j ( h i F i + 1 , j x − 2 F i + 1 , j ) / F i + 1 , j . (29) Combining conditions in Eqs. 28, 29 results in the following sufficient conditions for the positivity of S ( x , y j ) in inequality form: γ i , j > M a x { 0 , − α i , j ( h i F i , j x + 2 F i , j ) / F i , j , β i , j ( h i F i + 1 , j x − 2 F i + 1 , j ) / F i + 1 , j } . (30) Likewise, the remaining three boundary curves S ( x , y j + 1 ) ,   S ( x i , y ) , and S ( x i + 1 , y ) are positive if the following sufficient conditions are satisfied.

S ( x , y j + 1 ) > 0 , if and only if B 1 > 0 and B > 0 : γ i , j + 1 > − α i , j + 1 ( h i F i , j + 1 x + 2 F i , j + 1 ) / F i , j + 1 , (31) γ i , j + 1 > β i , j + 1 ( h i F i + 1 , j + 1 x − 2 F i + 1 , j + 1 ) / F i + 1 , j + 1 . (32)

  S ( x i , y ) , if and only if C 1 > 0 and C 3 > 0 : γ ^ i , j > − α ^ i , j ( h ^ j F i , j y + 2 F i , j ) / F i , j , (33) γ ^ i , j > β ^ i , j ( h ^ j F i , j + 1 y − 2 F i , j + 1 ) / F i , j + 1 . (34) S ( x i + 1 , y ) , if and only if D 1 > 0 and D 3 > 0 : γ ^ i + 1 , j > − α ^ i + 1 , j ( h ^ j F i + 1 , j y + 2 F i + 1 , j ) / F i + 1 , j , (35) γ ^ i + 1 , j > β ^ i + 1 , j ( h ^ j F i + 1 , j + 1 y − 2 F i + 1 , j + 1 ) / F i + 1 , j + 1 . (36)