Construction of C1 Rational Bi-Quartic Spline With Positivity-Preserving Interpolation: Numerical Results and Analysis

From the observed datasets, we should be able to produce curve surfaces that have the same characteristics as the original datasets. For instance, if the given data are positive, then the resulting curve or surface must be positive on entire given intervals, i.e., everywhere. In this study, a new partial blended rational bi-quartic spline with C1 continuity is constructed through the partially blended scheme. This rational spline is defined on four corners of the rectangular meshes. The sufficient condition for the positivity of rational bi-quartic spline is derived on four boundary curve networks. There are eight free parameters that can be used for shape modification. The first-order partial derivatives are estimated by using numerical techniques. We also show that the proposed scheme is local quadratic reproducing such that it can exactly reproduce the quadratic surface. We test the proposed scheme to interpolate various types of positive surface data. Based on statistical indicators such as the root mean square error (RMSE) and coefficient of determination (R 2), we found that the proposed scheme is on par with some established schemes. In fact, it requires less CPU times (in seconds) to generate the interpolating surface on rectangular meshes. Furthermore, by combining the statistical indicators’ result and graphically visualizing the test functions, the proposed method has the capability to reconstruct very comparable smoothing interpolating positive surfaces compared to some existing schemes. This finding is significant in producing a better interpolating surface for computer graphics applications since the proposed scheme has a smaller error compared with existing schemes.


INTRODUCTION
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][2][3][4][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, nonlinear 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 positivitypreserving 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 positivitypreserving 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 shapepreserving 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 nonnegativity 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 shapepreserving 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 positivitypreserving 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. 2) To derive the sufficient condition for the positivity of the rational bi-quartic spline on each boundary, we also derive the quadratic polynomial reproducing. 3) Finally, to compare the performance between the proposed scheme for positivity-preserving interpolation and some established schemes such as Abbas et al. [50] and Karim et al. [36], the root mean square error (RMSE) and coefficient of determination (R 2 ) are used as statistical goodness-of-fit measurements [53].
Furthermore, we identified several advantages of the proposed partially blended rational bi-quartic spline interpolation for positivity-preserving:  [40], and Tian et al. [41] No free parameters for shape refinement. Therefore, we cannot control the final interpolating curves and surfaces. 2 Hussain et al. [20,42] Cannot produce positivity-preserving interpolation on entire given intervals as shown by Harim et al. [29]. 3 Asim and Brodlie [43], Asim et al. [11], Brodlie and Butt [6], Butt and Brodlie [7], Fiorot and Tabka [44], and Peng et al. [45] Additional knots are inserted. This will increase the CPU times. 4 Fristch and Butland [46], Fristch and Carlson [47], and Wang and Tan [28,48,49] Shape-preserving with modification of derivatives. Cannot interpolate first derivatives. 5 Hussain et al. [20] The interpolating curve and surface are C 0 continuous. Thus, the surface is not smooth and not suitable for shape-preserving interpolation. 6 Abbas et al. [50,51] and Karim et al. [36] RMSE is a bit higher and smaller R 2 value. They did not prove whether their schemes are local quadratic polynomial reproducing or not. 7 Wu et al. [13,14] Requires solving an optimization problem that will increase unnecessary CPU times (in seconds) to produce the curves. 8 Ibraheem et al. [52] The interpolating surface is not smooth as well as not visually pleasing as seen in their study. 9 Ashraf et al. [43] Use the subdivision scheme to preserve the shape of the data. 1) There are 12 parameters in the description of the rational biquartic 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. 2) Unlike existing schemes such as Karim et al. [36], Abbas et al. [50,51], Hussain and Ali [37], Hussain and Hussain [40], Tian et al. [41], and Sarfraz et al. [9], the proposed scheme is local quadratic reproducing. This is one of the main novelties of the proposed scheme. 3) Based on statistical indicators, i.e., root mean square error (RMSE) and coefficient of determination (R 2 ), the proposed partially blended rational bi-quartic spline is the best compared with the works of Abbas et al. [50] and Karim et al. [36]. In fact, the proposed scheme has a smaller RMSE error and higher R 2 compared with Brodlie et al.'s scheme [7]. 4) Furthermore, the proposed scheme is direct and not involving any iterative scheme via the subdivision scheme as discussed by Ashraf et al. [43].
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.

Derivative Estimation
In most applications, the first derivative needs to be estimated by using some numerical methods. Three common methods are the arithmetic mean method (AMM), geometric mean method (GMM), and harmonic mean method (HMM) [54].

Arithmetic Mean Method
The first derivatives d i at the data point x i , i 0, 1, ..., n, can be estimated as follows.
At both end points x 0 and x n , the first derivatives are [54] Meanwhile, at the interior points, x i , i 1, 2, ..., n − 1, the values of d i are estimated using the following formula: Formulas in Equations (1)-(3) can be extended to 3D data as follows.
The first partial derivatives are given as follows.
The first partial derivative in the x direction is given as The first partial derivative in the y direction is given as for i 0, 1, 2, ..., n, j 1, 2, ..., m − 1. Meanwhile, the twists (mixed derivative) at the interior points can be estimated by

Geometric Mean Method
The first derivatives also can be calculated by using the GMM.
.., n}, at the end points x 0 and x n , the first derivatives are given as with Δ n,n−2 fn−fn−2 x n −x n−2 . Meanwhile, at interior points, x i , i 1, 2, ..., n − 1, the values of d i are estimated by Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 555517 The advantages of using this choice are the values of d i are always positive. This method is suitable for monotonically increasing data.

Rational Quartic Spline Interpolation
. < x n and the first derivative d i , at the respective point. Then, the RQS scheme with three shape control parameters α i > 0, β i > 0 , and c i ≥ 0 on the interval is given by [29,34,35] where hi , and θ x−xi hi . The RQS scheme defined in Eq. 14 gives the C 1 continuity at the knots and satisfies the following conditions: From Eq. 15 and after some derivations, the unknown A i , B i , and C i are given as [34] A

Partially Blended Rational Bi-Quartic Spline Interpolation
We construct a new rational bi-quartic spline based on the partially blended scheme of Casciola and Romani [16]. There are four curve networks on each edge of the rectangular domain. Figure 1 shows the arrangement of the partially blended rational bi-quartic spline interpolation on a rectangular domain. Let the function S(x, y) be arranged over a rectangular domain ,j indicate the data values and the first partial derivatives with respect to x and y at the node (x i , y j ). The partially blended rational bi-quartic function over each with S(x, y j ), S(x, y j+1 ), S(x i , y), and S(x i+1 , y) are rational quartic functions defined in Equations (18) Figure 1). They are defined as S x, y j with with with The partial derivatives F x i,j and F y i,j are estimated by using the AMM given in Equations (4)-(10).

Local Control Shape Property Analysis
The proposed rational bi-quartic spline defined by Eq. 17 has local control properties. Some observations can be made as follows: 1) The proposed scheme is symmetrical, i.e., the surface can be constructed either in x and y directions or in y and x directions. This is one of the main properties in CAGD. 2) When the parameters α i,j α i,j+1 β i,j β i,j+1 1, α i,j α i+1,j β i,j β i+1,j 1, and c i,j c i,j+1 c i,j c i+1,j 2, the partially blended rational bi-quartic spline in Eq. 17 reduced to a standard bi-quartic polynomial spline. Usually, in this case, the resulting interpolating surface is not shape-preserving.
3) When c i,j → ∞ or α i,j , β i,j → 0, the partially blended rational bi-quartic spline defined by Eq. 17 is reduced to a bilinear surface, i.e., flat surface, since all four curve networks are reduced to a straight line (see [34]): • Similarly, we can obtain the following observations: Based on Equations (22)- (25), we conclude that the final interpolating surface is a bilinear surface. Theorem 1. The partially blended rational bi-quartic spline interpolation defined in Eq. 17 is C 1 continuous everywhere provided that all free parameters are positive.
Theorem 2. The proposed rational bi-quartic spline interpolation given in Eq. 17 is shape-preserving provided that all four curve networks defined in Eqs. 18-21 are shapepreserving.
Theorem 3. The partially blended rational bi-quartic spline S(x, y) defined by Eq. 17 is a C 1 −continuous degree-nine piecewise rational surface (degree seven in the numerator and degree two in the denominator).
The proof of Theorems 1-3 can be obtained by generalizing Proposition 3 in the study of Casciola and Romani [16].
The behavior of local control of the proposed partially blended rational bi-quartic spline based on the observations given in Equations (21)-(25) is given in Example 1.
Example 1. Data from the following function are truncated to five decimal places [50]: F 1 x,y e − ( x 2 +y 2 ) 15 sin(x) + cos y + 0.33, 0 ≤x,y ≤6. (26) Figure 2 shows the interpolating surface by using the proposed scheme for data listed in Table 2. Figure 2A shows the interpolating surface with α i,j β i,j 1, α i,j β i,j 1, and c i,j c i,j 2. Meanwhile, Figure 2B shows the surface with α i,j β i,j 1, α i,j β i,j 1, and c i,j c i,j 50. Clearly, the surfaces are bilinear or flat. A similar effect can be seen in Figure 2C, i.e., when α i,j β i,j 0.01, α i,j β i,j 0.01, and c i,j c i,j 1. Meanwhile, Figure 2D shows the behavior of the interpolating surface when α i,j , β i,j , α i,j , and β i,j become a larger number, i.e., equal to 50. Finally, Figure 2E shows the surface when c i,j and c i,j become a smaller number, i.e., equal to 0.005. Furthermore, the parameters cannot be negative because the interpolating surface will become a rational surface, i.e., having poles on the surface. This will result in the discontinuities of the surface. This fact can be seen in Figure 3. Thus, the parameters are restricted to having positive values only. This is in line with the existing works such as Karim et al. [36], Abbas et al. [50,51], Hussain and Ali [37], Hussain and Hussain [40], Tian et al. [41], and Sarfraz et al. [9].
Finally, as the parameters α i,j , β i,j , α i,j , and β i,j are increased and c i,j and c i,j are fixed, we found that the rational interpolant is stable, i.e., the produced surfaces are almost identical to each other. Meanwhile, as the parameters c i,j and c i,j keep increasing while parameters α i,j , β i,j , α i,j , and β i,j are fixed, we also found that the rational interpolant is stable. These results are shown in Figure 4.

Local Quadratic Reproducing
The proposed partially blended rational bi-quartic spline defined in Eq. 17 is local quadratic reproducing. To prove this property, we use numerical simulations by using the following quadratic test function:  Table 3 shows the values of RMSE, R 2 , and maximum absolute error for each figure. Based on Table 3, we found that as c i,j → 0 or α i,j , β i,j → ∞, the proposed rational bi-quartic spline is local quadratic reproducing. This property is important since we can use our scheme to reconstruct any conic sections. Meanwhile, the works by Karim et al. [36] and Abbas et al. [50] are not local quadratic reproducing.

POSITIVITY-PRESERVING INTERPOLATION
The rational bi-quartic spline defined in Eq. 17 is not positivitypreserving 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.
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.
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: Combining conditions in Eqs. 28, 29 results in the following sufficient conditions for the positivity of S(x, y j ) in inequality form: S(x i , y), if and only if C 1 > 0 and C 3 > 0: S(x i+1 , y), if and only if D 1 > 0 and D 3 > 0: Combining the inequality conditions stated in Eqs. 31-36 leads us to the following conditions: Then from Eqs. 30, 37, the sufficient conditions for the positivity of all the four curve networks S(x, y j ), S(x, y j+1 ), S(x i , y), and S(x i+1 , y) are obtained. This completes the proof.The sufficient conditions in Eq. 27 can be rewritten as where the parameters Theorem 5. The positive partially blended rational bi-quartic spline S(x, y) constructed by using shape parameters calculated using Eq. 38 is a C 1 −continuous degree-nine piecewise rational surface with positivity-preserving properties with respect to the positive datasets.
Proof. This is a consequence from Theorem 3.The following algorithm can be used to implement the proposed positivitypreserving interpolation.
Step 1: Calculate the first partial derivative values F x i,j , F y i,j , F xy i,j by using numerical techniques, i.e., AMM.
Step 3: Construct the positive surfaces by substituting all the required parameters into the partially blended rational bi-quartic spline function in Eq. 17. Repeat Steps 2 and 3 to generate a different positive surface through different positive datasets.

RESULTS AND DISCUSSION
We test the proposed scheme by using two positive datasets. Example 2 uses the same data as in Example 1. All computational, simulation, and graphical results are obtained by using MATLAB Version 2019a installed on Windows 10, AMD Ryzen 3 2200G with Radeon Vega Graphics 3.50 GHz. As validation, we calculate RMSE and R 2 values and compare the results with those of Karim et al. [36] and Abbas et al. [50].
(39) Figure 6A shows the default bi-quartic Hermite spline for the positive data given in Table 4. Figures 6B,C show the xz-view and yz-view for Figure 6A, respectively. Figure 6D shows positivity-preserving by using the proposed rational bi-quartic spline with α i,j β i,j 1, α i,j β i,j 1. Meanwhile, Figures 6E,F show the xz-view and yz-view for Figure 6D, respectively. Figure 6G shows the other view for Figure 6D. Clearly, the interpolating surfaces are positive on entire given intervals. Meanwhile, Figures 6H,I show the positive interpolating surfaces using Karim et al. [36] and Abbas et al. [50] schemes. Example 3. Positive data from the following function are truncated to four decimal places [51]: Figure 7A shows the default bi-quartic Hermite spline for the positive data given in Table 5. Figures 7B,C show the xz-view and yz-view for Figure 7A, respectively. Figure 7D shows positivitypreserving by using the proposed rational bi-quartic spline with α i,j β i,j 1, α i,j β i,j 1. The parameters c i,j and c i,j are obtained by applying Theorem 4, i.e., Eq. 27. Meanwhile, Figures  8E,F show the xz-view and yz-view for Figure 8D, respectively. Figure 7G shows the other view for Figure 7A. Meanwhile, Figure 7H shows the other view for the surface in Figure 8D.
We found that the interpolating surfaces are positive everywhere.
Figures 7I,J show the interpolating surfaces using Karim et al. [36] and Abbas et al. [50] schemes. We compare the performance between the proposed partially blended rational bi-quartic spline interpolation and the works of Abbas et al. [50] and Karim et al. [36] by calculating the root mean square error (RMSE) and coefficient of determination (R 2 ). Tables 6, 7 summarize the main results. It can clearly be seen that the proposed rational bi-quartic spline interpolation is on par with the works of Abbas et al. [50] and Karim et al. [36]. Based on R 2 , for both examples, the proposed scheme is the best. Furthermore, according to Renka and Brown [53] criteria, the proposed scheme is excellent since the R 2 value is in between 0.91 and 0.997. Meanwhile, based on the RMSE value, we found that the proposed scheme is the best when the parameters are α i,j β i,j 0.1 and α i,j β i,j 0.1. From all numerical simulations, we found that the proposed scheme requires less computation times (in seconds) compared with Abbas et al. [50] and Karim et al. [36] schemes. The speedup is around 10% faster than the other schemes. Furthermore, based on graphical results, all schemes are looking competent to each other. Most probably, the proposed scheme is the best for both examples. Based on this, we conclude that the rank from the best scheme to the lowest is as follows: The Proposed Scheme > Karim et al. [36] > Abbas et al. [50] Our final example shows that the proposed rational bi-quartic spline for positivity-preserving has the capability to  reduce the CPU times (in seconds) to construct the interpolating surface compared with interpolation without positivitypreserving. Table 8 shows the positive data from the following function: F 3 x, y 0.25 Figure 8 shows the examples of various interpolating surfaces for data given in Table 8. Figure 8A shows the interpolating surface when the free parameters are α i,j β i,j 1, α i,j β i,j 1, and c i,j c i,j 2. Figures 8B,C show the xz-view and yz-view for Figure 8A.
Clearly, the surfaces are not positive everywhere. Figure 8D shows the positive interpolating surface after applying Theorem 4. The free parameters are α i,j β i,j 1, and α i,j β i,j 1. Meanwhile, c i,j and c i,j are obtained through Eq. 27. Figures 8E,F show the xz-view and yz-view for Figure 8D. From both graphs, the surfaces are positive and lie above both x and y axes. Figure 8G shows another interpolating surface with parameters α i,j β i,j 0.1, α i,j β i,j 0.1. However, the surface is not positive as seen in Figure 8H. Finally, Figure 8I shows the interpolating surface using the proposed scheme. From Figure 8J, it is seen that the interpolating surface shown in Figure 8I is positive everywhere. Therefore, Theorem 4 is confirmed to produce a positive interpolating surface on entire given intervals. This is in contrast with what was claimed by Qin et al. [38]. Therefore, the proposed scheme indeed has reduced the number of mathematical derivations needed to produce a positive rational interpolant. This is very significant in numerical analysis aspects. We only require lesser computation time compared with the work of Qin et al. [38]. CPU times (in seconds) are measured by using MATLAB tic and toc commands.
As seen in Table 9, by applying the positivity-preserving interpolation scheme to the given data, the CPU times (in seconds) needed to construct the interpolating surface are reduced. This is true for all set of parameters. We also noticed that the proposed sufficient condition for the positivity of the rational quartic interpolant is guaranteed to produce a positive interpolating surface everywhere, i.e., on entire given intervals. However, we also found that, by applying positivity-preserving interpolation, R 2 values will be reduced, while at the same time, the RMSE value will be increased. However, the values are not much different.
Final Remark: Qin et al. [38] have proposed a new approach to derive the sufficient condition for the positivity of the partially blended rational spline. They claimed that their condition is guaranteed to produce a positive interpolating surface everywhere. However, their method is mathematically more complicated to apply since it will increase the unnecessary computation time to produce the positive surface. If we apply the Qin et al. [38] method, the sufficient condition for the positivity of the rational quartic spline will take the following form:  Parameters' arrangement: A-α i, j β i, j 1, α i, j β i, j 1, and c i, j c i, j 2. B-αi,j β i, j 2.5, αi,j β i, j 2.5, and c i, j c i, j 0.5. C-αi, j β i, j 0.5, αi, j β i, j 0.5, and c i, j c i, j 0.5. D-αi, j β i, j 0.1, αi, j β i, j 0.1, and c i, j c i, j 0.5.
It is very complicated to solve Eq. 43 until Eq. 46. This is because we need to simplify all equations to form a new rational function, i.e., quintic numerator and quadratic denominator. Compared with our derivation for the positivity of the rational quartic spline which is direct and easy as shown in Eq. 28 until Eq. 36. Furthermore, based on all numerical and graphical results presented in this study, we found that the sufficient condition on all four curve networks is sufficient to produce positive interpolating surfaces on the entire given domain and not just on the four boundary curves only. This fact can be seen clearly in the final example. When we applied the proposed scheme, we found that, without positivity-preserving interpolation, the z values are negative. Meanwhile, after we apply the sufficient condition stated in Theorem 4, we obtain the interpolating surface as positive everywhere as seen in Figures 6-8 as well as in Table 9. Therefore, we encouraged the user to use our proposed rational bi-quartic spline for positivity-preserving interpolation with the following merits: 1) It will guarantee to produce a positive interpolating surface on entire given intervals and not just on the four boundary curve networks as claimed by Qin et al. [38]. From all numerical examples, we found that, after the positivity conditions are applied, all points on the surface are positive. 2) It can reproduce an exactly quadratic polynomial unlike Abbas et al. [50], Karim et al. [36], and Qin et al. [38]. 3) Overall, based on R 2 and RMSE values, the proposed scheme is the best compared with some existing schemes.

CONCLUSION
In this study, a new partially blended rational bi-quartic spline interpolation with C 1 continuity is constructed for positivitypreserving application. The sufficient condition for the positivity of the rational quartic spline is derived on four boundary curves defined on rectangular meshes. This condition results in eight free parameters for shape modification. The proposed scheme is tested for many positive surface datasets. Based on numerical analysis and graphical comparison, the proposed scheme is the best, i.e., higher R 2 and smaller RMSE. We also found that the proposed scheme has the capability to reduce the CPU times (in seconds) when we want to construct a positive interpolating surface. Furthermore, the proposed positivity-preserving interpolation scheme is guaranteed to deliver a positive surface on the entire given interval. The derivation for the positivity sufficient condition in this study is easier than that in the study of Qin et al. [38]. By using our proposed scheme, we have reduced the required numbers of mathematical derivations, and our scheme is easy to use. On top of that, the proposed scheme is local quadratic polynomial reproducing. This is significant when the user wants to reconstruct any quadric surfaces. Future studies will be focusing on the construction of general shape-preserving for scattered data arranged on triangular meshes as well as its applications in image processing and image denoising such as in the studies of Karim and Saaban [55], Karim et al. [56], Walther and Schmidt [57], and Zulkifli et al. [58].

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, and further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
SA, AS, and VN curated the data, involved in formal analysis, obtained the resources, and wrote the original draft. SA and VN performed the methodology. SA acquired the funding and investigated the data. VN involved in project administration. SA and AS ran the software, validated and visualized the results, and reviewed and edited the paper.  PP: positivity-preserving. Note: sets of parameters. I-αi,j β i,j 1, αi,j β i,j 1, and c i,j c i,j 2. II-α i,j β i,j 1, α i,j β i,j 1, and c i,j c i,j 3. III-αi,j β i,j 0.5, αi,j β i,j 0.5, and c i,j c i,j 0.5. IV-αi,j β i,j 0.1, αi,j β i,j 0.1, and c i,j c i,j 2.