Scattered Data Interpolation Using Quartic Triangular Patch for Shape-Preserving Interpolation and Comparison with Mesh-Free Methods

: Scattered data interpolation is important in sciences, engineering, and medical-based problems. Quartic B é zier triangular patches with 15 control points (ordinates) can also be used for scattered data interpolation. However, this method has a weakness; that is, in order to achieve C 1 continuity, the three inner points can only be determined using an optimization method. Thus, we cannot obtain the exact B é zier ordinates, and the quartic scheme is global and not local. Therefore, the quartic B é zier triangular has received less attention. In this work, we use Zhu and Han’s quartic spline with ten control points (ordinates). Since there are only ten control points (as for cubic B é zier triangular cases), all control points can be determined exactly, and the optimization problem can be avoided. This will improve the presentation of the surface, and the process to construct the scattered surface is local. We also apply the proposed scheme for the purpose of positivity-preserving scattered data interpolation. The sufficient conditions for the positivity of the quartic triangular patches are derived on seven ordinates. We obtain nonlinear equations that can be solved using the regula-falsi method. To produce the interpolated surface for scattered data, we employ four stages of an algorithm: (a) triangulate the scattered data using Delaunay triangulation; (b) assign the first derivative at the respective data; (c) form a triangular surface via convex combination from three local schemes with C 1 continuity along all adjacent triangles; and (d) construct the scattered data surface using the proposed quartic spline. Numerical results, including some comparisons with some existing mesh-free schemes, are presented in detail. Overall, the proposed quartic triangular spline scheme gives good results in terms of a higher coefficient of determination (R 2 ) and smaller maximum error (Max Error), requires about 12.5% of the CPU time of the quartic B é zier triangular, and is on par with Shepard triangular-based schemes. Therefore, the proposed scheme is significant for use in visualizing large and irregular scattered data sets. Finally, we tested the proposed positivity-preserving interpolation scheme to visualize coronavirus disease 2019 (COVID-19) cases in Malaysia.


Introduction
Scattered data interpolation and approximation are still active research topics in computer-aided design (CAD) and geometric modeling [1][2][3][4][5][6][7][8][9]. This is because engineers and scientists often face the Symmetry 2020, 12 problem of how to produce smooth curves and surfaces for the raw data obtained from experiments or observations. This is where scattered data interpolation can be used to assist them. To construct smooth curves and surfaces, some mathematical formulations are required. This can be achieved using functions which are well-established, such as the Bézier, B-spline, and radial basis functions (RBFs). All these methods are guaranteed to produce smooth curves and surfaces. The formulation problem in scattered data interpolation can be described as follows: Given functional data x i , y i , z i , i = 1, 2, . . . , N construct a smooth C 1 surface z = F(x, y) such that z i = F x i , y i , i = 1, 2, . . . , N To solve the above problem, there are many methods that can be used, such as meshless methods (e.g., radial basis functions (RBFs) and many types of Shepard's families). However, some meshless schemes are global. Fasshauer [10] gave details on many meshless methods to solve the problems arising in scattered data interpolation and approximation, as well as partial differential equations. Beyond that, another approach that can be used to solve the problem is the triangulation of the given data points. Then, the Bézier or spline triangular can be used to construct a piecewise smooth surface with some degree of smoothness, such as C 1 or C 2 . The Shepard triangular can also be used to produce a continuous surface from irregular scattered data. For instance, Cavoretto et al. [6], Dell'Accio and Di Tommaso [11], and Dell'Accio et al. [12,13] have discussed the application of the Shepard triangular for surface reconstruction. However, their schemes require more computation time in order to produce the interpolated surfaces.
Crivellaro et al. [14] applied RBFs to reconstruct 3D scattered data via new algorithms, which involves an adaptive multi-level interpolation approach based on implicit surface representation. The least squares approximation is used to remove the noise that appears in the scattered data. Chen and Cao [15] employed neural network operators of a logistic function through translations and dilation. Meanwhile, Bracco et al. [2] considered scattered data fitting using hierarchical splines where the local solutions are represented in variable degrees of the polynomial spline. Zhou and Li [16] studied scattered noise data by extending the weighted least squares method via triangulating the data points. Zhou and Li [17] discussed the scattered data interpolation for noisy data by using bivariate splines defined on triangulation. Qian et al. [18] also considered scattered data interpolation by using a new recursive algorithm based on the non-tensor product of bivariate splines. Liu [19] constructed local multilevel scattered data interpolation by proposing a new idea (i.e., nested scattered data sets), and scaled the compactly supported RBFs. Borne and Wende [3] also considered the meshless scheme based on definite RBFs for scattered data interpolation. In their study, they applied the domain decomposition methods to produce a symmetric-saddle point system. Joldes et al. [20] modified the moving least squares (MLS) methods by integrating the polynomial bases to solve the scattered data interpolation problem. Brodlie et al. [5] discussed the constrained surface interpolation by using the Shepard interpolant. The solution to the problem is obtained by solving some optimization. Lai and Meile [21] discussed scattered data interpolation by using nonnegative bivariate triangular splines to preserve the shape of the scattered data. Schumaker and Speleers [22] considered the nonnegativity preservation of scattered data by using macro-element spline spaces including Clough-Tocher macro-elements. Furthermore, they also give general results for range-restricted interpolation. Karim et al. [23] discussed the spatial interpolation for rainfall data by employing cubic Bézier triangular patches to interpolate the scattered data. Karim et al. [24] have constructed a new type of cubic Bézier-like triangular patches for scattered data interpolation. Karim and Saaban [25] constructed the terrain data using cubic Ball triangular patches [23]. In this study, they show that the scattered data interpolation scheme by Said and Rahmat [26] is not C 1 everywhere. Thus, a new condition for C 1 continuity is derived. The final surface is C 1 and provides a smooth surface. Feng and Zhang [27] proposed piecewise bivariate Hermite interpolations based on triangulation. They applied the scheme for large scattered data sets to produce high-accuracy surface reconstruction. Sun et al. [28] constructed bivariate rational interpolation defined on a triangular domain for scattered data lying on a parallel line. They only considered a few data sets, and it was not tested for large data sets. By using a rational spline, the computation time increases. Bozzini et al. [4] proposed a polyharmonic spline to approximate the noisy scattered data.
The main motivation of the present study is described in the following paragraphs. In triangulationbased approaches to scattered data interpolation, cubic Bézier triangular or quintic Bézier triangular patches are the common methods. The quartic Bézier triangular has received less attention due to the need to solve optimization problems in order to calculate the Bézier ordinates. This approach increases the computation time. There are four steps in constructing a surface using a triangulation method: (a) triangulate the domain by using Delaunay triangulation; (b) specify the first partial derivative at the data points (sites); (c) assign the control points or ordinates for each triangular patch; and finally (d) the surface is constructed via a convex combination scheme. Goodman and Said [29] constructed a suitable C 1 triangular interpolant for scattered data interpolation using a convex combination scheme between three local schemes. Their work is different from that of Foley and Opitz [30]. However, both studies developed a C 1 cubic triangular convex combination scheme. Said and Rahmat [26] constructed a scattered data surface using cubic Ball triangular patches [31,32] with the same approach as in Goodman and Said [29]. Based on the numerical results, their scheme gave the same results as cubic Bézier triangular patches. The main advantages of cubic Ball triangular patches are that the required computation is 7% less when compared with the work of Goodman and Said [29]. This is what has been claimed by References [26,29]. However, in the work of Karim and Saaban [25], it was proved that Said and Rahmat [26] is not C 1 continuous everywhere, and Karim and Saaban [25] found that the [26] scheme produced the same surface for scattered data interpolation when the inner coefficient was calculated by using Reference [29]. Hussain and Hussain [33] proposed the rational cubic Bézier triangular for positivity-preserving scattered data interpolation. They claimed that their proposed scheme is C 1 positive everywhere. However, from their results, it is possible that their scheme may not be positive everywhere. Chan and Ong [7] considered range-restricted interpolation using a cubic Bézier triangular comprising three local schemes. All the schemes were implemented by estimating the partial derivatives at the respective knots using the method proposed by Goodman et al. [34].
Other than the use of cubic Ball and cubic Bézier triangular patches for scattered data interpolation, there are some studies that have utilized quartic Bézier triangular and rational quartic Bézier triangular patches for scattered data interpolation. For instance, Saaban et al. [35] constructed C 1 (or G 1 ) scattered data interpolation based on the quartic Bézier triangular. Piah et al. [36] considered C 1 range-restricted positivity-preserving scattered data interpolation by using the quartic Bézier triangular. They employed an optimization method (i.e., the minimized sum of squares) to calculate the inner Bézier points proposed in Saaban et al. [35]. Hussain et al. [37] extended this idea to construct convexity-preserving scattered data interpolation. Hussain et al. [38] constructed a new scattered data interpolation scheme by using the rational quartic Bézier triangular. They applied it to positivity-preserving interpolation. However, to achieve C 1 continuity, we still need to solve some optimization problems. This is the main weakness of quartic Bézier triangular patches when applied to scattered data interpolation. Some good surveys on scattered data interpolation can be found in [39][40][41][42][43].
The present study aims to answer the following research questions: a. Can we construct a scattered data interpolation scheme by using quartic triangular patches but without an optimization method?
b. How can we produce a C 1 surface (everywhere)? c. Is the proposed scheme better than some existing schemes in terms of CPU time, coefficient of determination (R 2 ), and maximum error?
To answer these research questions, we will use the quartic triangular basis initiated by Zhu and Han [44]. The main advantage of using this quartic basis is that it only requires ten control points to construct one triangular patch. This is the same as the number of control points in the cubic Bézier Symmetry 2020, 12, 1071 4 of 26 triangular patch. Thus, in order to construct C 1 scattered data interpolation using the quartic spline triangular, we can employ the Foley and Opitz [30] cubic precision scheme to calculate the inner ordinates. With this, the optimization problem required in a quartic triangular basis will be avoided. Hence, this will show that the proposed scheme is local. Furthermore, the proposed scheme is different from the works of Lai and Meile [21] and Schumaker and Spellers [22], even though all schemes required triangulation of the given data in the first step.
Some contributions from the present study are described below: 1. The proposed scattered data interpolation scheme produces a C 1 surface without any optimization method like Piah et al. [36], Saaban et al. [35] and Hussain et al. [37,38].
3. Based on the CPU time needed to construct the surface, the proposed scheme is faster than quartic Bézier triangular patches. Thus, the reconstruction of scattered surfaces from large data sets can be performed in less time.
4. Furthermore, the proposed positivity-preserving scattered data interpolation is capable of producing a better interpolated surface than quartic Bézier triangular patches. This lies in contrast to scattered data schemes by Ali et al. [1], Draman et al. [9] and Karim et al. [24], which are not positivity-preserving interpolations.
This paper is organized as follows: In Section 2 we give a review of the triangular basis initiated by Zhu and Han [44], and the derivation of the quartic triangular basis with ten control points. Some graphical results are presented, as well as the construction of a local scheme comprising convex combination between three local schemes. The numerical results and the discussion are given in Section 3 with various numerical and graphical results, including a comparison with some existing schemes. Error analysis is also investigated in this section. The construction of the positive scattered data interpolant is discussed in Section 4. Meanwhile, numerical results for positivity-preserving scattered data interpolation are shown in Section 5. Conclusions and future recommendations are given in the final section.

Review of the Cubic Triangular Bases of Zhu And Han
Zhu and Han [44] proposed a new cubic triangular basis with three exponential parameters α, β, γ. Since we are dealing with triangulation, the barycentric coordinate (u, v, w) on the triangle T 1 with vertices V 1 , V 2 and V 3 is defined by u + v + w = 1, where u, v, w ≥ 0. Set the point inside the triangle as V(x, y) ∈ R 2 (as shown in Figure 1), which can be expressed as: Symmetry 2020, 12, x FOR PEER REVIEW 4 of 28 construct one triangular patch. This is the same as the number of control points in the cubic Bézier triangular patch. Thus, in order to construct scattered data interpolation using the quartic spline triangular, we can employ the Foley and Opitz [30] cubic precision scheme to calculate the inner ordinates. With this, the optimization problem required in a quartic triangular basis will be avoided. Hence, this will show that the proposed scheme is local. Furthermore, the proposed scheme is different from the works of Lai and Meile [21] and Schumaker and Spellers [22], even though all schemes required triangulation of the given data in the first step.
Some contributions from the present study are described below: 1. The proposed scattered data interpolation scheme produces a surface without any optimization method like Piah et al. [36], Saaban et al. [35] and Hussain et al. [37,38].
3. Based on the CPU time needed to construct the surface, the proposed scheme is faster than quartic Bézier triangular patches. Thus, the reconstruction of scattered surfaces from large data sets can be performed in less time.
4. Furthermore, the proposed positivity-preserving scattered data interpolation is capable of producing a better interpolated surface than quartic Bézier triangular patches. This lies in contrast to scattered data schemes by Ali et al. [1], Draman et al. [9] and Karim et al. [24], which are not positivitypreserving interpolations.
This paper is organized as follows: In Section 2 we give a review of the triangular basis initiated by Zhu and Han [44], and the derivation of the quartic triangular basis with ten control points. Some graphical results are presented, as well as the construction of a local scheme comprising convex combination between three local schemes. The numerical results and the discussion are given in Section 3 with various numerical and graphical results, including a comparison with some existing schemes. Error analysis is also investigated in this section. The construction of the positive scattered data interpolant is discussed in Section 4. Meanwhile, numerical results for positivity-preserving scattered data interpolation are shown in Section 5. Conclusions and future recommendations are given in the final section.

Review of the Cubic Triangular Bases of Zhu And Han [44]
Zhu and Han [44] proposed a new cubic triangular basis with three exponential parameters , , . Since we are dealing with triangulation, the barycentric coordinate , , on the triangle with vertices , and is defined by 1, where , , 0. Set the point inside the triangle as , ∈ (as shown in Figure 1), which can be expressed as: (1) Zhu and Han's triangular basis functions satisfy the following properties: Non-negativity: . For more details on the other properties, please refer to Zhu and Han [44]. Zhu and Han's triangular patches with three parameters α, β, and γ, and control points b ijk , i + j + k = 3 are defined as Figure 2 shows the Zhu and Han ordinates, and Figure 3 shows one patch of the Zhu and Han triangular with α = β = γ = 3 (i.e., cubic Bézier triangular).

Quartic Zhu and Han Triangular Patches
From Equation (2), let 4, then we obtain the following ten quartic basis functions defined on the triangular domain:

Quartic Zhu and Han Triangular Patches
From Equation (2), let α = β = γ = 4, then we obtain the following ten quartic basis functions defined on the triangular domain: (4) Figure 4 shows the quartic triangular basis on the triangular domain. (4) Figure 4 shows the quartic triangular basis on the triangular domain.  Thus, the quartic Zhu and Han triangular patch can be defined by The main advantage of Zhu and Han's quartic is that it only requires ten control points to construct one triangular patch; meanwhile, the quartic Bézier triangular will require 15 control points to produce one patch. Furthermore, when the quartic Bézier triangular is used for scattered data interpolation, an optimization method is required to produce the interpolated surface, as discussed in Saaban et al. [35], Piah et al. [36] and Hussain et al. [37,38]. However, if we apply the proposed quartic triangular patches for scattered data interpolation, the optimization is not required since we can employ the cubic precision scheme of Foley and Opitz [30] to construct a C 1 interpolated surface everywhere. So far, this is the first study to apply a a quartic triangular basis but with ten control points for scattered data interpolation. Figure 5a shows examples of quartic Zhu and Han, and Figure 5b shows the quartic Bézier triangular patch.

Scattered Data Interpolation Using Quartic Zhu and Han Triangular Patches
To apply the quartic triangular patch defined in Section 2.2 for scattered data, we use the local scheme comprising a convex combination between three local schemes , , and [1,9,24] such that: , , , The local scheme , 1,2,3 is obtained by replacing in (5) with to ensure the condition is satisfied. Given the vertex of the triangle (i.e., , , and ), the derivative along (see Figure 6)-that is, the edge connecting two points and -is defined as [1,9,24,29]: which can be simplified as 1 4 (7)

Scattered Data Interpolation Using Quartic Zhu and Han Triangular Patches
To apply the quartic triangular patch defined in Section 2.2 for scattered data, we use the local scheme comprising a convex combination between three local schemes K 1 , K 2 , and K 3 [1,9,24] such that: The local scheme K i , i = 1, 2, 3 is obtained by replacing b 111 in (5) with b i 111 to ensure the C 1 condition is satisfied. Given the vertex of the triangle (i.e., F(V 1 ) = b 300 , F(V 2 ) = b 030 , and F(V 3 ) = b 003 ), the derivative along e jk (see Figure 6)-that is, the edge connecting two points x j − y j and (x k − y k )is defined as [1,9,24,29]: which can be simplified as Similarly, the other five ordinates are calculated as follows: Symmetry 2020, 12, x FOR PEER REVIEW 9 of 28 The remaining , 1,2,3 is obtained by using the cubic precision of Foley and Opitz [30] as shown in Figure 7. For complete derivation, the reader can refer to [30]. In order to achieve continuity along all edges, the following equations must be satisfied: The remaining b i 111 , i = 1, 2, 3 is obtained by using the cubic precision of Foley and Opitz [30] as shown in Figure 7. For complete derivation, the reader can refer to [30]. The remaining , 1,2,3 is obtained by using the cubic precision of Foley and Opitz [30] as shown in Figure 7. For complete derivation, the reader can refer to [30]. In order to achieve continuity along all edges, the following equations must be satisfied: In order to achieve C 1 continuity along all edges, the following equations must be satisfied: To find c 1 111 in (13) and (14), we need to add these equations together. Thus, we obtain Similarly, with Equations (15) and (16), we obtain Now we establish the theorem for the main result. Theorem 1. The local scheme defined by (6) is a rational function with degree 7, that is, degree five in numerator and degree two in denominator with C 1 continuity everywhere. It has the following form: (18) and the barycentric coordinate satisfies u + v + w = 1.
The following Algorithm 1 can be used to implement the proposed scheme.

Algorithm 1 (Scattered Data Interpolation)
Step 1: Input scattered data points; Step 2: Estimate the partial derivative at the data points by using [25]; Step 3: Triangulate the domain of the data points; Step 4: Calculate the boundary control points using Equations (7)-(12); Step 5: Calculate inner control points for the local scheme, b i 111 , i = 1, 2, 3 by using the cubic precision method as in Foley and Opitz [30]; Step 6: Construct the interpolated surface using the convex combination method of three local schemes defined by (6); Step 7: Calculate CPU time (in seconds), R 2 , and maximum error. Repeat steps 1 through 6 for the other scattered data sets.
Below we give the theorem for scattered data interpolation by using quartic Bézier triangular patches.
Let the quartic Bézier triangular patch (n = 4) with barycentric coordinates u, v, w be expressed as and c ijk are the Bézier ordinates of S. Let the total number of triangles in the whole triangular mesh be n t and the total number of interior edges be n e . S(x,y) which will minimize the functional I(S(x,y)) leads to the optimization problem x T Qx + ex + h, subject to the C 1 continuity constraint: where Q is a 6n t × 6n t sparse matrix, e is a 1 × 6n t row vector, x is a 6n t × 1 column vector consisting of unknown ordinates b m h is a real constant, A is a 3n e × 6n t (3n e ≤ 6n t ) coefficient matrix, x is a 6n t × 1 unknown column vector consisting of the remaining ordinates b 220 , b 202 , b 022 , b 211 , b 121 and b 112 to be determined for the entire triangular mesh, and b is a 3n e × 1 constant column vector. The optimization problem stated in (18) was solved using the optimization toolbox in MATLAB 2017 on Intel ® Core ™ i5-8250U 1.60 GHz.
Note that the optimization problem in Theorem 2 is obtained by using a minimized principal curvature norm with respect to the C 1 continuity constraint, which results in a global method for scattered data interpolation. Meanwhile, by using the proposed scheme in this study, the resulting surface is local.

Results and Discussion for Scattered Data Interpolation
We tested the proposed scheme using two well-known test functions F 1 (x, y) and F 2 (x, y): F 2 (x, y) = (1.25 + cos cos (5.4y)) 6 + 6(3x − 1) 2 We implemented the proposed scheme using MATLAB 2017 version on Intel ® Core ™ i5-8250U 1.60 GHz. MATLAB coding was developed based on Algorithm 1. About 25 MATLAB functions were used to obtain all the results.
We chose 36 data point samples in the domain [0, 1] × [0, 1], as shown in Table 1. Figure 8 shows the Delaunay triangulation for the data. Figure 6 shows examples of surface interpolation for both functions. Comparing Figures 9 and 10, the surfaces produced by the proposed scheme visually look smoother than the surfaces obtained from the quartic Bézier triangular of Saaban et al. [35] and Piah et al. [36]. Figure 10 shows an example of scattered data interpolation using quartic Bézier triangular patches.
To validate the proposed scheme, we calculated the maximum error (Max Error) and coefficient of determination (COD; i.e., R 2 ) for both functions and compared them with those obtained for quartic Bézier triangular for three different numbers of points i.e., 100, 65, and 36 for both functions F 1 (x, y) and F 2 (x, y). Functions 1 and 2 represent F 1 (x, y) and F 2 (x, y), respectively. Table 2 shows the error analysis for both tested functions by using (a) quartic Zhu and Han and (b) quartic Bézier triangular. Meanwhile, Table 3 shows CPU time in seconds. From Table 2, we can see that the proposed quartic triangular patches for scattered data interpolation gave smaller Max Error values than the quartic Bézier triangular. Additionally, the proposed scheme gave higher R 2 values for all numbers of data points (100, 65, and 36). From Table 3, the proposed scheme required less CPU time than the quartic Bézier. For instance, for 100 data points, the proposed scheme only required 0.71 s for data from function F 1 (x, y) and 0.42 s for data from function F 2 (x, y), compared with the quartic Bézier which requires 5.6 s and 3.57 s for 100 data points from functions F 1 (x, y) and F 2 (x, y), respectively. Thus, the proposed scheme in this study gave very good results, and was better at treating scattered data than using the quartic Bézier triangular proposed by Piah et al. [36], Saaban et al. [35], and Hussain et al. [37,38]. We conclude that the proposed scheme required less CPU time than the quartic Bézier triangular. This reduction of CPU time consumption is an advantage when the goal is to construct a surface with thousands of data points or big data.    Table 1.  To validate the proposed scheme, we calculated the maximum error (Max Error) and coefficient of determination (COD; i.e., R 2 ) for both functions and compared them with those obtained for quartic Bézier triangular for three different numbers of points i.e., 100, 65, and 36 for both functions , and , . Functions 1 and 2 represent , and , , respectively. Table 2 shows the error analysis for both tested functions by using (a) quartic Zhu and Han and (b) quartic Bézier triangular. Meanwhile, Table 3 shows CPU time in seconds. From Table 2, we can see that the proposed quartic triangular patches for scattered data interpolation gave smaller Max Error values than the quartic Bézier triangular. Additionally, the proposed scheme gave higher R 2 values for all numbers of data points (100, 65, and 36). From Table 3, the proposed scheme required less CPU time than the quartic Bézier. For instance, for 100 data points, the proposed scheme only required 0.71 s for data from function , and 0.42 s for data from function , , compared with the quartic Bézier which requires 5.6 s and 3.57 s for 100 data points from functions , and , , respectively. Thus, the proposed scheme in this study gave very good results, and was better at treating scattered data than using the quartic Bézier triangular proposed by Piah et al. [36], Saaban et al. [35], and Hussain et al. [37,38]. We conclude that the proposed scheme required less CPU time than the quartic Bézier triangular. This reduction of CPU time consumption is an advantage when the goal is to construct a surface with thousands of data points or big data.   This can be seen clearly from Table 3. With the largest number of data points, the CPU time for the proposed scheme was approximately 12.5% that of the CPU time required for the quartic Bézier triangular based scheme. This is very significant, especially when the user wants to render and reconstruct surfaces obtained from very dense data sets. Many studies in scattered data interpolation usually involve the use of Shepard-type interpolants such as Shepard triangular schemes for scattered data interpolation [6,[11][12][13]. We also implemented the Shepard triangular to the same data sets as listed in Table 1. Tables 4 and 5 show the error analysis for the schemes of Cavoretto et al. [6], Dell'Accio et al. [12,13], and Dell'Accio and Di Tommaso [11]. Based on CPU time, for all tested data sets, the proposed scheme was faster than the schemes in [6,[11][12][13], except for the case with 100 data points for the data from function F 1 (x, y). Considering the Max Error, the proposed scheme was better than all four schemes except for the case with 100 data points from function F 2 (x, y). Therefore, we can conclude that the proposed scheme is better than quartic triangular patch and the Shepard triangular based schemes [6,[11][12][13].  Table 5. CPU time (in seconds).

Positivity-Preserving Scattered Data Interpolation
In this section, we apply the proposed scheme discussed in the previous section to preserve the positivity of scattered data sets. To do this, first we derive the sufficient condition for the positivity of the quartic triangular spline defined in (5). Finally, the rational corrected scheme defined by (19) will be used to construct a positive surface with C 1 continuity.
To derive the sufficient condition for the positivity of the quartic spline triangular patch, we adopted a similar approach to Saaban et al. [35] and Piah et al. [36]. Assume that the quartic ordinates at the vertices are strictly positive such b 300 , b 030 , b 003 > 0. Let A = b 300 , B = b 030 , C = b 003 , therefore A, B, C > 0 (see Figure 11). Meanwhile, let the other ordinates have the same value, that is −t < 0 where t > 0. Thus, Equation (5) becomes Symmetry 2020, 12, x FOR PEER REVIEW 15 of 28 In this section, we apply the proposed scheme discussed in the previous section to preserve the positivity of scattered data sets. To do this, first we derive the sufficient condition for the positivity of the quartic triangular spline defined in (5). Finally, the rational corrected scheme defined by (19) will be used to construct a positive surface with C 1 continuity.
To derive the sufficient condition for the positivity of the quartic spline triangular patch, we adopted a similar approach to Saaban et al. [35] and Piah et al. [36]. Assume that the quartic ordinates at the vertices are strictly positive such , , 0. Let , , , therefore , , 0 (see Figure 11). Meanwhile, let the other ordinates have the same value, that is 0 where 0. Thus, Equation (5) becomes decreases. We want to find the value of when the minimum value of , , 0 . By taking first partial derivatives, we will obtain the following: , , (23) Figure 11. Quartic triangular ordinates arrangement for positivity preservation.
From (23) we can observe that when t = 0 then P(u, v, w) > 0. Meanwhile as t increases, P(u, v, w) decreases. We want to find the value of t = t 0 when the minimum value of P(u, v, w) = 0. By taking first partial derivatives, we will obtain the following: The minimum value of P(u, v, w) occurs when ∂P ∂u − ∂P ∂v = 0 and ∂P ∂u − ∂P ∂w = 0 or equivalently.
From Equation (25) we have Hence: Since u + v + w = 1, we obtain the following relations: Substituting this value into (23) we obtain the minimum value of P(u, v, w) i.e., P(u, v, w) Let s = 1/t, then Then Equation (27) can be written as G(s) = 1, s ≥ 0. This equation can be solved by using regula-falsi method with suitable choice of initial guess.
Since A, B, C > 0 and s ≥ 0 then G (s) < 0 and G (s) > 0. (see Figure 12). Thus, the curve is convex on that region. Let X = max(A, B, C) and Y = min(A, B, C), then the following holds   Now we establish the main theorem for positivity preservation using the proposed scheme. Some observations from Theorem 3 can be made as follows: Such that G 26 X ≥ 1 and G 26 Y ≤ 1. Figure 12 shows the example of the relative locations of 26 X and 26 Y and s 0 . Now we establish the main theorem for positivity preservation using the proposed scheme. Theorem 3. Consider the quartic triangular patch P(u, v, w) with vertex A = b 300 , B = b 030 , C = b 003 , such that A, B, C > 0. If the remaining quartic triangular ordinates are equal to −t 0 where t 0 = 1 s 0 is a unique solution to (28), then P(u, v, w) ≥ 0 for all u, v, w ≥ 0 and u + v + w = 1.
Some observations from Theorem 3 can be made as follows: Therefore, as A → ∞, then G(s) → 1 (s+1) 1/3 . Hence, s 0 → 0 and therefore t 0 → ∞. Thus, the ordinate values are unbounded compared with the work of Chan and Ong [7] in which the Bézier ordinates are bounded by a lower bound −1/3.

Remark 1.
The sufficient condition for the positivity of the quartic triangular patch developed in this study is the same as the sufficient condition for the quartic Bézier triangular patch developed in Saaban et al. [35]. The main difference is that the proposed quartic polynomial only requires ten control points (or ordinates) as compared to the quartic Bézier triangular which requires 15 control points and involves some optimization problems as shown in Saaban et al. [35] and Hussain et al. [37,38]. Therefore, the proposed positivity preservation using a quartic triangular patch requires less computation time than some established schemes for scattered data interpolation.
The final construction of the positive scattered surface is described below: 1. Input positive scattered data points; 2. Triangulate the scattered data using Delaunay triangulation; 3. Assign the first partial derivative at the respective data sites and adjust if necessary, to provide the positivity preservation; 4. The C 1 triangular surface is constructed via convex combination between three local schemes; 5. Repeat Steps 1 through 4 for other positive scattered data sets.

Numerical Results and Discussion for Positivity-Preserving Scattered Data Interpolation
After we derive the sufficient condition for the positivity of quartic triangular patch, the final C 1 scattered data scheme for positivity preservation can be written as follows: with We test the proposed scheme by using four well-known test functions given below: The positive test functions F 1 , F 2 , F 3 , and F 4 were evaluated on 36, 33, 26, and 100 node points respectively (Tables 6-9) where all function values were greater than or equal to zero. The nodes of 36 and 33 points were defined on a rectangular domain (Figure 13a,b), while the 26-and 100-point nodes were defined on a sparse non-rectangular domain (Figure 13c,d). Tables 8 and 9 show examples of irregular scattered data sets. Table 6. Value of F 1 on 36 node points.  For test function F 1 as the data in Table 6, the interpolated surface did not preserve the positivity of the original surface for the C 1 Zhu and Han quartic (from Theorem 1), as shown in Figure 14a with calculated min F 1 (x, y) (x,y)∈D = −0.039975. Observe that these surfaces cross the xy-plane at a number of places. After applying positivity-preserving methods from Theorem 3, the result is shown in Figure 14b, where the interpolated surfaces lie above or on the xy-plane min F 1 (x, y) For test function as the data in Table 6, the interpolated surface did not preserve the positivity of the original surface for the C 1 Zhu and Han quartic (from Theorem 1), as shown in Figure  14a with calculated , , ∈ -0.039975 . Observe that these surfaces cross the xy-plane at a number of places. After applying positivity-preserving methods from Theorem 3, the result is shown in Figure 14b, where the interpolated surfaces lie above or on the xy-plane , , ∈ 0. For test function as the data in Table 7, the interpolated surface did not preserve the positivity of the original surface for the C 1 Zhu and Han quartic as shown in Figure 15a, with calculated , , ∈ -0.039975. These surfaces cross the xy-plane at a number of places. Using the proposed positivity-preserving methods, the interpolated surface lies above or on the xy-plane, as shown in Figure 15b, with calculated , 0.0072657. For test function F 2 as the data in Table 7, the interpolated surface did not preserve the positivity of the original surface for the C 1 Zhu and Han quartic as shown in Figure 15a, with calculated min F 2 (x, y) (x,y)∈D = −0.039975. These surfaces cross the xy-plane at a number of places. Using the proposed positivity-preserving methods, the interpolated surface lies above or on the xy-plane, as shown in Figure 15b, with calculated min F 2 (x, y) For test function as the data in Table 7, the interpolated surface did not preserve the positivity of the original surface for the C 1 Zhu and Han quartic as shown in Figure 15a, with calculated , , ∈ -0.039975. These surfaces cross the xy-plane at a number of places. Using the proposed positivity-preserving methods, the interpolated surface lies above or on the xy-plane, as shown in Figure 15b, with calculated , , ∈ 0.0072657.
For the third test function defined on a sparse non-rectangular domain (data in Table 8), the interpolated surface did not preserve the positivity, as shown in Figure 16a where the surface crosses below the xy-plane with , , ∈ -0.0053288 and the positivity-preserving interpolated surface using the proposed scheme is shown in Figure 16b where the surface lies above or on the xyplane, with calculated .
For the third test function defined on a sparse non-rectangular domain (data in Table 8), the interpolated surface did not preserve the positivity, as shown in Figure 16a where the surface crosses below the xy-plane with min F 3 (x, y) (x,y)∈D = −0.0053288 and the positivity-preserving interpolated surface using the proposed scheme is shown in Figure 16b where the surface lies above or on the xy-plane, with calculated.  Table 8): (a) without positivity preserved; (b) with positivity preserved from Theorem 3.
The interpolated surface of the Zhu and Han C 1 quartic without positivity preservation is given in Figure 17a, with calculated , , ∈ 0.67634, while the positivity-preserved surface lying above the xy-plane is illustrated in Figure 17b with calculated , , ∈ 0.000074928. The interpolated surface of the Zhu and Han C 1 quartic without positivity preservation is given in Figure 17a, with calculated min F 4 (x, y) (x,y)∈D = −0.67634, while the positivity-preserved surface lying above the xy-plane is illustrated in Figure 17b with calculated min F 4 (x, y) (x,y)∈D = 0.000074928.
The interpolated surface of the Zhu and Han C 1 quartic without positivity preservation is given in Figure 17a, with calculated , , ∈ 0.67634, while the positivity-preserved surface lying above the xy-plane is illustrated in Figure 17b with calculated , , ∈ 0.000074928.
We also calculated the CPU time (in seconds), maximum error, and coefficient of determination (R 2 ) for the positivity-preserving scattered data interpolation as shown in Tables 10 and 11. Once again, the proposed scheme was superior to the quartic Bézier triangular patch. For positivity preservation in scattered data interpolation with dense data sets (i.e., 100 data points with 1697 points of evaluation), the proposed scheme only required 0.5168 s, compared with the quartic Bézier which required 18.5996 s. This is about 36 times faster than the times obtained by the schemes of Saaban et  Table 9): (a) without positivity preserved; (b) with positivity preserved from Theorem 3.
We also calculated the CPU time (in seconds), maximum error, and coefficient of determination (R 2 ) for the positivity-preserving scattered data interpolation as shown in Tables 10 and 11. Once again, the proposed scheme was superior to the quartic Bézier triangular patch. For positivity preservation in scattered data interpolation with dense data sets (i.e., 100 data points with 1697 points of evaluation), the proposed scheme only required 0.5168 s, compared with the quartic Bézier which required 18.5996 s. This is about 36 times faster than the times obtained by the schemes of Saaban et al. [35] and Piah et al. [36]. Roughly, the proposed scheme only required about 2.78% of the CPU times of schemes [35,36]. This is very significant when we want to visualize thousands of scattered data points.   Figures 18 and 19 show the example of surface interpolation for COVID-19 scattered data listed in Table 12. Figure 18 shows the interpolated surface without positivity preservation. Figure 19 shows the interpolated surface after we applied the positivity-preserving scheme. Clearly, Figure 19 is suitable for the relevant agency to visualize the number of COVID-19 cases. Then, they can prepare any contingency plan for the spread of COVID-19. They could also try to minimize the spread of COVID-19. This is very crucial, since at the time of writing there are no available vaccines to cure the patients.

Conclusions
Zhu and Han [44] proposed new cubic Bernstein-Bézier basis functions defined on a triangular domain. We implemented quartic triangular bases (with ten control points) for scattered data interpolation. This new quartic basis makes it possible to avoid the optimization problem that appears when the quartic Bézier triangular is used for scattered data interpolation. From the results, we can see that the proposed scheme in this study outperformed the quartic Bézier triangular, having the smallest maximum error, higher R 2 , and requiring only 12.5% of the CPU time needed by the quartic Bézier triangular scheme. This is very significant, especially when the goal is to reconstruct surfaces from large scattered data sets. Furthermore, based on a comparison against the Shepard triangular for scattered data, the proposed scheme was also superior to the schemes of Cavoretto et al. [6], Dell'Accio et al. [12,13] and Dell'Accio and Di Tommaso [11]. Finally, we constructed a positive interpolant based on the proposed quartic triangular spline to preserve the positivity of scattered data. Numerical results suggest that the proposed scheme is better than existing schemes, especially in terms of CPU time-our proposed scheme requires less computation time than positivity schemes proposed by Piah et al. [36] and Saaban et al. [35]. Finally, we implemented our proposed positivity-preserving interpolation to visualize COVID-19 cases in Selangor State and Klang Valley, Malaysia. The resulting surfaces were smooth and positive everywhere. Future works will focus on the construction of a quintic Zhu and Han spline for scattered data interpolation with quintic precision as well as shape-preserving interpolation (e.g., positivity-preserving and range-restricted interpolation). This can be achieved by extending the main idea from Karim et al. [46]. Another potential study could be a comparison between the use of a CPU and a graphical processing unit (GPU) for large scattered data sets. Finally, the proposed scheme can also be applied to visualize large sets of scattered data, such as from geophysical data, medical imaging, and total COVID-19 cases around the world.