*2.3. 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*<sup>3</sup> [1,9,24] such that:

$$P(\mu, v, w) = \frac{vw\mathbf{K\_1} + \mu v\mathbf{K\_2} + \mu v\mathbf{K\_3}}{vw + \mu w + \mu v}, \text{ } \mu + v + w = \mathbf{1} \tag{6}$$

The local scheme *Ki*, *i* = 1, 2, 3 is obtained by replacing *b*<sup>111</sup> in (5) with *b<sup>i</sup>* <sup>111</sup> to ensure the *<sup>C</sup>*<sup>1</sup> 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 *ejk* (see Figure 6)—that is, the edge connecting two points - *xj* − *yj* and (*xk* − *yk*) is defined as [1,9,24,29]:

$$\frac{\partial P}{\partial \boldsymbol{\varepsilon}\_{\vec{j}\vec{k}}} = \left(\mathbf{x}\_{k} - \mathbf{x}\_{\vec{j}}\right) \frac{\partial F}{\partial \boldsymbol{\alpha}} + \left(y\_{k} - y\_{\vec{j}}\right) \frac{\partial F}{\partial y}$$

Thus

$$b\_{210} = F(V\_1) + \frac{1}{4} \frac{\partial P}{\partial \varepsilon\_3}(V\_1)$$

which can be simplified as

$$b\_{210} = b\_{300} + \frac{1}{4} \left( (\mathbf{x}\_2 - \mathbf{x}\_1) F\_\mathbf{x}(V\_1) + (y\_2 - y\_1) F\_y(V\_1) \right) \tag{7}$$

Similarly, the other five ordinates are calculated as follows:

$$b\_{201} = b\_{300} - \frac{1}{4} \left( (\mathbf{x}\_1 - \mathbf{x}\_3) F\_\mathbf{x} \left( V\_1 \right) + (y\_1 - y\_3) F\_\mathbf{y} \left( V\_1 \right) \right) \tag{8}$$

$$b\_{021} = b\_{030} + \frac{1}{4} \left( (\mathbf{x}\_3 - \mathbf{x}\_2) F\_x(V\_2) + (y\_3 - y\_2) F\_y(V\_2) \right) \tag{9}$$

$$b\_{120} = b\_{030} - \frac{1}{4} \left( (\mathbf{x}\_2 - \mathbf{x}\_1) F\_\mathbf{x} (V\_2) + (y\_2 - y\_1) F\_\mathbf{y} (V\_2) \right) \tag{10}$$

$$b\_{102} = b\_{003} + \frac{1}{4} \left( (\mathbf{x}\_1 - \mathbf{x}\_3) F\_x(V\_3) + (y\_1 - y\_3) F\_y(V\_3) \right) \tag{11}$$

$$b\_{012} = b\_{003} - \frac{1}{4} \left( (\mathbf{x}\_3 - \mathbf{x}\_2) F\_x(V\_3) + (y\_3 - y\_2) F\_y(V\_3) \right) \tag{12}$$

**Figure 6.** Side-vertex blending.

The remaining *bi* <sup>111</sup>, *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].

**Figure 7.** Two adjacent quartic triangular patches.

In order to achieve *C*<sup>1</sup> continuity along all edges, the following equations must be satisfied:

$$c\_{201} = r^2 b\_{210} + 2st b\_{021} + 2rs b\_{120} + s^2 b\_{030} + 2rt b\_{111}^2 + t^2 b\_{012} \tag{13}$$

$$c\_{210} = r^2 b\_{201} + 2st b\_{012} + 2rt b\_{102} + s^2 b\_{021} + 2rs b\_{111}^1 + t^2 b\_{003} \tag{14}$$

$$b\_{210} = \mathfrak{u}^2 \mathfrak{c}\_{201} + 2\nu \mathfrak{wc}\_{012} + 2\mu \mathfrak{wc}\_{102} + \mathfrak{v}^2 \mathfrak{c}\_{021} + 2\mu \mathfrak{wc}\_{111}^1 + \mathfrak{w}^2 \mathfrak{c}\_{003} \tag{15}$$

$$b\_{201} = \mathfrak{u}^2 \mathfrak{c}\_{210} + 2\mathfrak{v}\mathfrak{wc}\_{021} + 2\mathfrak{u}\mathfrak{wc}\_{120} + \mathfrak{v}^2 \mathfrak{c}\_{030} + 2\mathfrak{u}\mathfrak{wc}\_{111}^1 + \mathfrak{w}^2 \mathfrak{c}\_{012} \tag{16}$$

To find *c*<sup>1</sup> <sup>111</sup> in (13) and (14), we need to add these equations together. Thus, we obtain

$$\begin{array}{l} \mathbf{c}\_{111}^{1} = \frac{1}{2\underline{\mathbf{u}}(\mathbf{v} + \mathbf{w})} \Big( b\_{201} + b\_{210} - \mu^2 (\mathbf{c}\_{210} + \mathbf{c}\_{201}) - \mathbf{v}^2 (\mathbf{c}\_{030} + \mathbf{c}\_{021}) \Big) \\ + \mathbf{w}^2 (\mathbf{c}\_{012} + \mathbf{c}\_{003}) - 2\underline{\mathbf{v}} \mathbf{w} (\mathbf{c}\_{021} + \mathbf{c}\_{012}) - \underline{\mathbf{u}} \mathbf{v} \mathbf{c}\_{120} - 2\underline{\mathbf{u}} \mathbf{u} \mathbf{c}\_{102}. \end{array}$$

Similarly, with Equations (15) and (16), we obtain

$$\begin{array}{l} b\_{111}^1 = \frac{1}{2r(s+t)} \Big( c\_{201} + c\_{210} - r^2 (b\_{210} + b\_{201}) - s^2 (b\_{030} + b\_{021}) \Big) \\ + t^2 (b\_{012} + b\_{003}) - 2st (b\_{021} + b\_{012}) - rsb\_{120} - 2rtb\_{102}. \end{array}$$

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*<sup>1</sup> *continuity everywhere. It has the following form:*

$$P(u,v,w) = \sum\_{\substack{i+j+k=3 \ i,j,k\neq 1}} b\_{ijk} B\_{i,j,k}^3(u,v,w) + 6uvw \left(a\_1 b\_{111}^1 + a\_2 b\_{111}^2 + a\_3 b\_{111}^3\right) \tag{17}$$

*with*

$$a\_1 = \frac{vw}{vw + uw + uv'}, \\ a\_2 = \frac{uw}{vw + uw + uv'}, \\ a\_3 = \frac{uv}{vw + uw + uv} \tag{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, *bi* <sup>111</sup>, *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), R2, 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.

**Theorem 2.** *C1 quartic Bézier triangular patches using minimized sum of squares of principal curvatures* [35].

Let the quartic Bézier triangular patch (*n* = 4) with barycentric coordinates *u*, *v*, *w* be expressed as

$$S(\boldsymbol{u}, \boldsymbol{v}, \boldsymbol{w}) = \sum\_{i+j+k=4} c\_{i\bar{j}k} B\_{i\bar{j}k}^4(\boldsymbol{u}, \boldsymbol{v}, \boldsymbol{w}) \tag{19}$$

where *B*<sup>4</sup> *ijk*(*u*, *<sup>v</sup>*, *<sup>w</sup>*) <sup>=</sup> 4! *i*!*j*!*k*! *ui vj v<sup>k</sup>* and *cijk* are the Bézier ordinates of S. Let the total number of triangles in the whole triangular mesh be *n*<sup>t</sup> and the total number of interior edges be *ne*. *S*(*x*,*y*) which will minimize the functional *I*(*S*(*x*,*y*)) leads to the optimization problem **x***<sup>T</sup> Q***x** + **ex** + *h*, subject to the *C*<sup>1</sup> continuity constraint:

$$\min \, \text{x}^T \mathbf{Q} \mathbf{x} + \epsilon \mathbf{x} + h \, \text{subject to } \mathbf{A} \mathbf{x} = \mathbf{b} \tag{20}$$

where *Q* is a 6*nt* × 6*nt* sparse matrix, **e** is a 1 × 6*nt* row vector, **x** is a 6*nt* × 1 column vector consisting of unknown ordinates - *bm* <sup>220</sup> *bm* <sup>211</sup> *bm* <sup>121</sup> *bm* <sup>202</sup> *bm* <sup>112</sup> *bm* 022 , *m* = 1, ... , *nt*, *h* is a real constant, *A* is a 3*n*<sup>e</sup> × 6*nt* (3*ne* ≤ 6*nt*) coefficient matrix, **x** is a 6*nt* × 1 unknown column vector consisting of the remaining ordinates *b*220, *b*202, *b*022, *b*211, *b*<sup>121</sup> and b112 to be determined for the entire triangular mesh, and **b** is a *3ne* × 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*<sup>1</sup> 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.
