Construction of new cubic Bézier-like triangular patches with application in scattered data interpolation

This paper discusses the functional scattered data interpolation to interpolate the general scattered data. Compared with the previous works, we construct a new cubic Bézier-like triangular basis function controlled by three shape parameters. This is an advantage compared with the existing schemes since it gives more flexibility for the shape design in geometric modeling. By choosing some suitable value of the parameters, this new triangular basis is reduced to the cubic Ball and cubic Bézier triangular patches, respectively. In order to apply the proposed bases to general scattered data, firstly the data is triangulated using Delaunay triangulation. Then the sufficient condition for C1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$C^{1}$\end{document} continuity using cubic precision method on each adjacent triangle is implemented. Finally, the interpolation scheme is constructed based on a convex combination between three local schemes of the cubic Bézier-like triangular patches. The detail comparison in terms of maximum error and coefficient of determination r2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$r^{2}$\end{document} with some existing meshfree methods i.e. radial basis function (RBF) such as linear, thin plate spline (TPS), Gaussian, and multiquadric are presented. From graphical results, the proposed scheme gives more visually pleasing interpolating surfaces compared with all RBF methods. Based on error analysis, for all four functions, the proposed scheme is better than RBFs except for data from the third function. Overall, the proposed scheme gives r2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$r^{2}$\end{document} value between 0.99920443 and 0.99999994. This is very good for surface fitting for a large scattered data set.


Introduction
Scattered data interpolation arises in many scientific and engineering fields. This method can be used to represent the observed or computed values of some physical quantities such as global temperature, rainfall distribution at some country or station, digital elevation, or the stress measurement in finite element methods. Furthermore, it can be used for spatial data interpolation.
There are two methods that can be used in scattered data interpolation i.e. meshfree and triangulation-based schemes. For instance, one of the simplest meshfree methods is Shepard's global surface scheme [2]. Many authors have improved the work of Shepard by proposing many new ideas on the extension of the original method. For instance, RBF, thin plate spline, and compactly supported positive definite function [2,10,11]. Crivellaro et al. [5] applied RBFs to reconstruct 3D scattered data. To achieve this, they proposed new algorithms involving an adaptive multi-level interpolating approach based on implicit surface representation and least square approximation to filter the noisy data. Liu [23] proposed local multilevel scattered data interpolation by employing nested scattered data sets and scaled the RBFs compactly. The method guarantees convergence. Majdisova and Skala [25,26] discussed the applications of RBFs for big geo data as well as finding the best basis for function approximations. A good survey on scattered data interpolation using meshfree methods and other methods can be found in Lodha and Franke [24], Franke and Nielson [10,11], and Franke [9]. MATLAB implementation for meshfree methods can be found in Fasshauer [7].
In triangulation-based approach for scattered data interpolation, cubic Bézier triangular or quantic Bézier triangular patches are the common methods that have been used. Quartic Bézier triangular has received less attention due to the need to apply an optimization problem in solving the scattered data interpolation problem. This increases the computation time. There are four steps in constructing surface by using a triangulation method: (a) start with triangulation of the domain by using Delaunay triangulation; (b) specify the first partial derivative at the data points (sites) [16]; (c) assign the control points or ordinates for each triangular patch; and (d) finally, construct the surface via convex combination scheme. Lawson in 1977 gave the idea on construction of the surface via triangulation-based approach [13]. Goodman and Said [13] constructed the suitable C 1 triangular interpolant for scattered data interpolation using a convex combination scheme from three local schemes. Their work is different from that of Foley and Opitz [8], even though both studies developed a C 1 cubic triangular convex combination scheme. Said and Rahmat [29] constructed a scattered data surface by using cubic Ball triangular patches of Goodman and Said [13]. Considering numerical results, the results are almost the same as by using cubic Bézier triangular patches except that the computation is less by 7% by using cubic Ball triangular [14,15]. Karim and Saaban [21] have proved that the scheme of Said and Rahmat [29] is not producing C 1 surface everywhere in the given triangular domain. Karim et al. [20] discussed spatial interpolation for rainfall scattered data by extending the results from Chan and Ong [3]. The cubic Bézier triangular patches with three local schemes are blended to produce C 1 surface everywhere. Sometimes in certain applications, the partial derivatives are not given, hence, they must be estimated. Thus, Goodman et al. [16] proposed a method to estimate the partial derivatives of the data for the scattered data interpolation. Throughout our study, we implement [16] scheme to estimate the partial derivatives at data sites.
Zhou and Li [31] also considered scattered noisy data by using bivariate splines on a triangulation domain. Chen and Cao [4] discussed scattered data approximation by employing neural network operators via translations and dilation of logistic function. Bracco et al. [1] discussed the scattered data fitting by utilizing hierarchical splines. They have extended the main idea from local least squares approximation where the local solutions are described in variable degree polynomial spline. Lai and Meile [22] discussed the scattered data interpolation by using bivariate splines with nonnegative property. The spline is constructed on a triangular domain i.e. the given data are triangulated first. Qian et al. [27] discussed scattered data interpolation based on the bivariate recursive polynomials defined on triangulation. Zhu et al. [32] constructed another version of triangular patch with three exponent parameters. Hussain and Hussain [17] constructed C 1 scattered data interpolation based triangulation by using side-vertex method where the rational cubic function is used to construct rational interpolant for each side of the triangle. Then the final scheme is a blend between three local schemes. Sarfraz et al. [30] extended the idea in [17] but by using different rational cubic functions. Both papers are quite similar to each other, and the results are also not much different. Hussain et al. [18] also discussed positivity-preserving scattered data interpolation by using cubic trigonometric spline functions. They have tested their scheme to two irregular scattered positive data sets. Unfortunately, in the paper, they showed different surface interpolation for the given data sets.
Even though the study on interpolation based on Bézier and Ball representation is already thirty years old, many researchers are still focusing on how to improve both representations by adding more flexibility to the control point to control the shape of curves or surfaces. For instance, [26] constructed an explicit parametric curve to be taken as the limitation curve of progressive iteration approximation (PIA) which can interpolate some scattered data points by using normalized totally positive (NTP) basis by specially choosing two kinds of NTP bases, Said-Bézier type generalized Ball (SBGB) basis and DP basis. Their results avoid the tedious calculation of the inverse matrix and hence will gain extensive application in reverse engineering. In 2013, [12] solved the parameterization problem for polynomial Bézier surfaces by applying the firefly algorithm, a powerful nature-inspired metaheuristic algorithm introduced recently to address difficult optimization problems. The method has been successfully applied to some illustrative examples of open and closed surfaces, including shapes with singularities. Their method performs very well, being able to yield the best approximating surface with a high degree of accuracy. But in order to obtain the final solution, we need to train the nodes; besidesit is time consuming for certain data.
The outcome of our current study is motivated by the works of Said [28], Goodman and Said [13], and Said and Rahmat [29]. We propose a new cubic Bézier-like triangular basis function with three parameters by extending the univariate cubic Bézier-like of Said [28]. Several properties of the new cubic Bézier-like triangular basis are derived. This new cubic triangular basis reduces to the cubic Ball and Bézier triangular bases with suitable choices of the three parameters. This new triangular basis is extended to the scattered data interpolation. The sufficient condition for C 1 continuity along the adjacent triangles with cubic precision method is constructed. Each triangular patch of the interpolating surface is constructed by using convex combination of three local schemes of cubic Bézier-like triangular patches. Several numerical results are presented including comparison with some existing schemes including meshfree methods such as radial basis functions (RBFs) i.e. linear, thin plate spline, Gaussian, and multiquadric.

Construction of new cubic Bézier-like triangular patches with shape parameters
In this section, a new Bézier-like triangular basis is constructed. Let the barycentric coordinate u, v, w on triangle T with vertices V 1 , V 2 , and V 3 be such that u + v + w = 1 and u, v, w ≥ 0. Any point V (x, y) ∈ R 2 inside the triangle (including at the vertices) ( Fig. 1) can Now we can establish the construction of new Bézier-like triangular patches with three parameters α, β, and γ as follows.
The following ten functions are a new cubic Bézier-like triangle basis on a triangular domain (2)

Definition 2 The basis function
can be simplified as Thus a new cubic Bézier-like triangular patch can be written as follows.
Since, for α, β, γ ∈ [0, 1], |i+j+k|=3 B 3 i,j,k (u, v, w) = 1 and B 3 i,j,k (u, v, w) ≥ 0, then the resulting cubic triangular patches will satisfy the convex hull as well as affine invariance. (e) Boundary property: When one of the barycentric coordinates is zero, for instance w = 0, and v = 1u, then a new cubic Bézier-like triangular basis will degenerate to a standard univariate Bézier-like basis function such that . Figure 4 shows the examples of few basis functions where α, β, and γ vary.

Scattered data interpolation using new Bézier-like triangular bases 4.1 Local scheme and C 1 continuity
The scheme comprising convex combination of three local schemes P 1 , P 2 , and P 3 is defined as follows: There are two methods that can be used to construct three local schemes P i , i = 1, 2, 3, i.e. those by Goodman and Said [13] and Foley and Opitz [8], where in [13] the crossderivative is used, meanwhile Foley and Opitz [8] have implemented the cubic precision method. The local scheme P i , i = 1, 2, 3, is obtained by replacing b 111 in (4) with b i 111 so that C 1 condition is satisfied only on the boundary e i of the triangle.
Assume that the barycentric coordinates at the vertices are given as V 1 = (1, 00), V 2 = (0, 1, 0), and V 3 = (0, 0, 1) where u + v + w = 1. The direction vectors e i , i = 1, 2, 3, on the side opposite to the vertex V i are given as Fig. 5 [21]). The directional derivatives for P(u, v, w) on the direction z = (z 1 , z 2 , z 3 ), where z 1 + z 2 + z 3 = 0, are given as Now, we describe the method to determine the edge coordinates on the triangle. From Fig. 5, the directional derivatives at V 1 i.e. along the edges e 2 and e 3 are calculated as follows: Vertices and edges of the triangle from [21] where F(V 1 ) = b 300 , F x (V 1 ) and F y (V 1 ) are given at the vertex V 1 . Similarly, we can obtain that F(V 2 ) = b 030 , F x (V 2 ) and F y (V 2 ) and F(V 3 ) = b 003 , F x (V 3 ) and F y (V 3 ). The partial derivatives are estimated by using Goodman et al. [16] method. Then, by simple derivation, the following are obtained: By substituting (7) to (8), we obtain Similarly, by considering the other directional derivatives at V 2 i.e. along the edges e 1 and e 3 as well as at V 3 i.e. along the edges e 1 and e 2 , the cubic Bézier-like ordinates b 021 , b 120 , b 102 , and b 012 are given as follows: Now, we just need to determine the inner Bézier-like ordinate b i 111 for each local scheme P i , i = 1, 2, 3. To achieve this, we have adopted the main idea from the scheme proposed by Foley and Opitz [8]. Consider the two adjacent triangles (as shown in Fig. 6), let the first triangle T 1 be represented in terms of the cubic Bézier-like triangular defined in (3) and T 2 be represented by There are three cases that can be considered as follows.
Note that the coefficients (s + t)r and u(v + w) in Eqs. (10)- (15) are not equal to zero because r, u < 0 and s + t, v + w > 0 [8]. Thus, if e i , i = 1, 2, 3, is a common edge to two tri-  In this case, the hybrid patch will be a standard cubic patch. By using all calculated Bézier ordinates of the three local schemes, the final interpolating surface P on each triangle can now be constructed. It is well described in the following paragraphs.

Graphical examples
The construction of scattered data interpolation using new Bézier-like triangular patches is described as follows: (a) Triangulate the domain by using Delaunay triangulation [13]. (b) Specify the derivatives at the data points, then assign Bézier-like ordinate (control point) values for each triangular patch by using [16]. (c) Generate the triangular patches for each of the triangle domains to form composite C 1 surface by using the local scheme defined by (16) with three parameters α, β, and γ . (d) Calculate the error measurement such as maximum error (Max. Error) and coefficient of determination (COD) i.e. r 2 . The final C 1 scheme can be written as with In Goodman and Said [13], the other version of convex combination is used such that The main difference between (18) and (17) is that the final degree of the rational patches obtained from (17) is rational quintic with quadratic denominator (5/4), while form (18) will give a rational heptic with quartic denominator (7/4). But both schemes require only 10 control points.
Remark 1 Although the local scheme given in (16) has singularities at the triangle vertices, the scheme in (16) still satisfies C 1 with removable singularities at the vertices. This is consequent from the work of Goodman and Said [13]. As an example of the surface interpolation, by using the proposed scheme, we use the following two well-known test functions:  The triangular domain defined on 36 data points in the domain of [0, 1] × [0, 1] is illustrated in Fig. 7, while the surface interpolation with different values of three parameters α, β, and γ defined on triangulation mesh shown in Fig. 7 is given in Fig. 8(a) for the test function F 1 (x, y) and in Fig. 9(a) for the test function F 2 (x, y) respectively. We compare the The interpolating surfaces are given in Fig. 8(b)-8(d) for the test function F 1 (x, y) and in Fig. 9(b)-9(d) for the test function F 2 (x, y). Based on the observation from the interpolating surfaces, the proposed scheme gives a very visually pleasing surface compared to all RBF techniques. This is significant in scattered data interpolation.

Error analysis
To test the capability of the new cubic Bézier-like triangular scheme for scattered data interpolation, we calculate the value of maximum error (Max. Error) and coefficient of determination (COD) i.e. r 2 based on the test functions F 1 and F 2 , together with the fol-   Fig. 7, while for the 65 and 100 data points, sets are given in Fig. 10(a) and (b), respectively. Tables 1, 2, 3 show the error analysis for all four functions with respective nodes i.e. 36, 65, and 100 points. We compare the performance between the proposed scheme and RBFs methods such as linear, thin plate spline, Gaussian, and multiquadric. The statistical goodness-fit used are (a) maximum error and (b) coefficient of determination r 2 . For 36 data points, the proposed scheme is better than RBFs functions except for the data from function F 3 . Similarly, for 65 and 100 nodes, the proposed scheme is better except for F 3 . This is understandable, since function F 3 looks like a Gaussian-type function. Therefore, maybe this is the main reason why RBFs are better for the data from that function. But if we refer to Table 1, the maximum error for function F 3 is 0.00662 (Gaussian) and 0.009586 (the proposed scheme), and the value of r 2 is not too much different. For more details on numerical results including graphical images for the proposed method including the comparison between the implementation using Goodman and Said [13] and Foley and Opitz [8] schemes, the reader can refer to Karim and Saaban [19]. Table 4 summarizes the main results for error analysis shown in Tables 1, 2, 3.
Our final example considers the irregular scattered data from Ibraheem et al. [18]. The total number of the data is 72 as shown in Table 5. Figure 11 shows the interpolating surfaces for irregular data obtained by using the proposed scheme with various values of α, β, and γ . Figure 11(a) shows the Delaunay triangulation for data in Table 5. Figure 1(b) shows the linear interpolant for the irregular data. Figure 11(c) to 11(e) show the surface interpolation with different values of the parameters.   Fig. 11(f ) shows the solid version for Fig. 11(c). Clearly, the produced surface is smooth and visually pleasing. Interestingly, the given data is positive, and the proposed scheme has the capability to preserve the positivity of the data without the need to apply any positivity-preserving as discussed in Ibraheem et al. [18], Hussain and Hussain [17], and Sarfraz et al. [30]. It should be noted that we cannot compare the results with the work of Ibraheem et al. [18], since in [18] their results are shown for different data sets, not for the data given in their paper.

Conclusion and recommendation
In this study, a new cubic triangular patch with three parameters is constructed. This new basis function is called cubic Bézier-like triangular patches. Some properties of the proposed basis are studied in detail. An application in scattered data interpolation is investigated in detail. We have adopted Goodman and Said [13] method in order to find all Bézier-like ordinates that will ensure that the sufficient condition for C 1 is satisfied. The ordinates b i 111 , i = 1, 2, 3, respectively for local scheme P i are obtained by using Foley and Opitz [8] method as discussed in Sect. 4. Throughout this study, we test the proposed scheme with regular and irregular scattered data. For regular data sets, we compare the performance with meshfree based methods i.e. RBFs in terms of maximum error, coefficient of determination r 2 , and visually pleasing for graphical displays. Based on the validation, the proposed scheme is better than RBFs methods since, for all data sets, the proposed scheme has smaller maximum error and higher r 2 value i.e. between 0.999204443 and 0.99999994. Furthermore, based on graphical images, the proposed scheme is more visually pleasing compared with all RBFs. For irregular scattered data, the proposed scheme can reconstruct the surface that is visually pleasing. We would like to stress that the main reason we have adopted the Foley and Opitz [8] method to calculate the inner ordinates b i 111 , i = 1, 2, 3, is that, it will give different surface for different values of the shape parameters (of the same scattered data). Meanwhile if we apply Goodman and Said [13] scheme to calculate the inner ordinates, we will obtain same interpolating surface for different values of the parameters. This results as well as comprehensive numerical comparison can be obtained in Karim and Saaban [19]. The proposed scheme can be used for surface reconstruction for cloud data i.e. thousands of data or more. Furthermore, by using the proposed scheme, we can preserve the positivity of the scattered data by manipulating the three parametric values. Work on shape-preserving interpolation is underway. Furthermore, we also can extend the main results in this study to construct more general triangulation surface based on Clough-Tocher splitting method. This will avoid the use of rational corrected as shown in formulas (16) and (17). Such results will be explored in the future.

Figure 11
Interpolating scattered data surface for irregular data