3.3.1. Process the Single-Sided Region

For the simple single-sided region, the polygon triangulation method based on constrained Delaunay triangulation is used to directly model it, as shown in Figure 6. For the complex single-sided region, we consider transforming it into a four-sided region by analyzing its shape characteristics and then performing the Coons surface interpolation. *Minerals* **2022**, *12*, x FOR PEER REVIEW 10 of 20

**Figure 6.** Schematic diagram of polygon triangulation: (**a**) the original polygon; (**b**,**c**) the process and the result of polygon triangulation. **Figure 6.** Schematic diagram of polygon triangulation: (**a**) the original polygon; (**b**,**c**) the process and the result of polygon triangulation.

#### 3.3.2. Process the Two-Sided and the Three-Sided Regions 3.3.2. Process the Two-Sided and the Three-Sided Regions

For the two-sided and the three-sided regions, we consider generating the virtual sides to transform them into four-sided regions. For the two-sided and the three-sided regions, we consider generating the virtual sides to transform them into four-sided regions.

As shown in Figure 7, for the two-sided region, two virtual edges are generated at two common intersections and on both sides. Then, transition lines are generated between and , which means virtual points are inserted on the two virtual sides. When performing interpolation calculation, the coordinate values of points on the same virtual side are the same. The coordinates of , = 1,2, … , are the same as . The coordinates of , = 1,2, … , are the same as . As shown in Figure 7, for the two-sided region, two virtual edges are generated at two common intersections *A* and *B* on both sides. Then, *n* transition lines are generated between *A* and *B*, which means *n* virtual points are inserted on the two virtual sides. When performing interpolation calculation, the coordinate values of *n* points on the same virtual side are the same. The coordinates of *A<sup>i</sup>* , *i* = 1, 2, . . . , *n* are the same as *A*. The coordinates of *B<sup>i</sup>* , *i* = 1, 2, . . . , *n* are the same as *B*.

> B1 B2 B3 Bi Bn−<sup>2</sup> Bn−<sup>1</sup> Bn

A1 A2 A3 Ai An−<sup>2</sup> An−<sup>1</sup> An

(b) (c) The intersection of the two-sided region The virtual point The vertex of the side The side of the two-sided region The virtual side The transition line

**Figure 7.** Schematic diagram of the two-sided region processing: (**a**) the original two-sided region; (**b**,**c**) the process and result of adding virtual sides and virtual points. The black dotted lines are

For the three-sided region, a virtual side is constructed at the intersection of any two sides. The transition lines are generated, whose number is the same as the number of vertices on the third side. Then, the same number of virtual points are inserted on the virtual side, and the coordinates of the virtual points are the same as that of the

For the five-sided region and regions that have more than five sides, first, search for the shortest side among all sides, and then search for the shorter side adjacent to this side for merging and connection. Repeat the above operation until this region is merged into

After the above preprocessing, all the regions are transformed into four-sided regions. Since the synchronous forward method is used for interpolation, the number of vertices on the opposite side is required to be the same. However, the vertices cannot be directly added to the sides because of the requirement of forming manifold graphics. In this paper, we consider translating the sides of the four-sided region to the middle according to the specified distance, which will generate four new sides. As shown in

A B

omitted lines.

(a)

intersection where the virtual side is inserted.

3.3.4. Process All the Four-Sided Regions

3.3.3. Process Other Regions

a four-sided region.

omitted lines.

the result of polygon triangulation.

3.3.2. Process the Two-Sided and the Three-Sided Regions

same as . The coordinates of , = 1,2, … , are the same as .

sides to transform them into four-sided regions.

(a) (b) (c) **Figure 6.** Schematic diagram of polygon triangulation: (**a**) the original polygon; (**b**,**c**) the process and

For the two-sided and the three-sided regions, we consider generating the virtual

As shown in Figure 7, for the two-sided region, two virtual edges are generated at two common intersections and on both sides. Then, transition lines are generated between and , which means virtual points are inserted on the two

**Figure 7.** Schematic diagram of the two-sided region processing: (**a**) the original two-sided region; (**b**,**c**) the process and result of adding virtual sides and virtual points. The black dotted lines are **Figure 7.** Schematic diagram of the two-sided region processing: (**a**) the original two-sided region; (**b**,**c**) the process and result of adding virtual sides and virtual points. The black dotted lines are omitted lines.

For the three-sided region, a virtual side is constructed at the intersection of any two sides. The transition lines are generated, whose number is the same as the number of vertices on the third side. Then, the same number of virtual points are inserted on the virtual side, and the coordinates of the virtual points are the same as that of the intersection where the virtual side is inserted. For the three-sided region, a virtual side is constructed at the intersection of any two sides. The transition lines are generated, whose number is the same as the number of vertices on the third side. Then, the same number of virtual points are inserted on the virtual side, and the coordinates of the virtual points are the same as that of the intersection where the virtual side is inserted.

#### 3.3.3. Process Other Regions

*3.4. Sided Region Modeling* 

*v0(0)*

*v0(1)*

*v0(2)*

*v0(j)*

*v0(n−3)*

*v0(n−2)*

*v0(n−1)*

3.3.3. Process Other Regions For the five-sided region and regions that have more than five sides, first, search for the shortest side among all sides, and then search for the shorter side adjacent to this side for merging and connection. Repeat the above operation until this region is merged into For the five-sided region and regions that have more than five sides, first, search for the shortest side among all sides, and then search for the shorter side adjacent to this side for merging and connection. Repeat the above operation until this region is merged into a four-sided region.

#### a four-sided region. 3.3.4. Process All the Four-Sided Regions

3.3.4. Process All the Four-Sided Regions After the above preprocessing, all the regions are transformed into four-sided regions. Since the synchronous forward method is used for interpolation, the number of vertices on the opposite side is required to be the same. However, the vertices cannot be directly added to the sides because of the requirement of forming manifold graphics. In this paper, we consider translating the sides of the four-sided region to the middle according to the specified distance, which will generate four new sides. As shown in After the above preprocessing, all the regions are transformed into four-sided regions. Since the synchronous forward method is used for interpolation, the number of vertices on the opposite side is required to be the same. However, the vertices cannot be directly added to the sides because of the requirement of forming manifold graphics. In this paper, we consider translating the sides of the four-sided region to the middle according to the specified distance, which will generate four new sides. As shown in Figure 8b, the same numbers of vertices are inserted on the opposite sides, which makes the same *t* of the vertices at the same position on the opposite sides. The *t* is the ratio of the distance from the vertice to the starting point to the side length. At last, the polygon triangulation method is used to model the newly formed sides and old sides, as shown in Figure 8c. *Minerals* **2022**, *12*, x FOR PEER REVIEW 11 of 20 Figure 8b, the same numbers of vertices are inserted on the opposite sides, which makes the same of the vertices at the same position on the opposite sides. The is the ratio of the distance from the vertice to the starting point to the side length. At last, the polygon triangulation method is used to model the newly formed sides and old sides, as shown in Figure 8c.

**Figure 8.** Schematic diagram of the side processing of the four-sided region: (**a**) the original foursided region; (**b**) the formed sides and the points inserted on them; (**c**) the result of the polygon triangulation of the new and old sides. **Figure 8.** Schematic diagram of the side processing of the four-sided region: (**a**) the original foursided region; (**b**) the formed sides and the points inserted on them; (**c**) the result of the polygon triangulation of the new and old sides.

After preprocessing, the four-sided regions can be modeled using bicubic Coons surface interpolation. Since the synchronous forward method is used to perform the

*u0(m−2)*

*u0(m−1)*

in the grid will be ×. As shown in Figure 9, set the first group of opposite sides respectively as and ଵ, the second group of opposite sides respectively as and ଵ.

*u0(0) u0(1) u0(2) u0(i) u0(m−3)*

*u0*

*v0 v1*

*Aij*

*u1*

**Figure 9.** Schematic diagram of the Coons surface interpolation.

#### *3.4. Sided Region Modeling 3.4. Sided Region Modeling*

triangulation of the new and old sides.

The original side

Figure 8c.

*Minerals* **2022**, *12*, x FOR PEER REVIEW 11 of 20

After preprocessing, the four-sided regions can be modeled using bicubic Coons surface interpolation. Since the synchronous forward method is used to perform the Coons surface interpolation, the numbers of the vertices on the opposite sides are the same, which are *m* and *n*. This means that the number of the formed interpolation points in the grid will be *m* × *n*. As shown in Figure 9, set the first group of opposite sides respectively as *u*<sup>0</sup> and *u*1, the second group of opposite sides respectively as *v*<sup>0</sup> and *v*1. After preprocessing, the four-sided regions can be modeled using bicubic Coons surface interpolation. Since the synchronous forward method is used to perform the Coons surface interpolation, the numbers of the vertices on the opposite sides are the same, which are and . This means that the number of the formed interpolation points in the grid will be ×. As shown in Figure 9, set the first group of opposite sides respectively as and ଵ, the second group of opposite sides respectively as and ଵ.

(a) (b) (c)

**Figure 8.** Schematic diagram of the side processing of the four-sided region: (**a**) the original foursided region; (**b**) the formed sides and the points inserted on them; (**c**) the result of the polygon

The inserted point

The new formed side

Figure 8b, the same numbers of vertices are inserted on the opposite sides, which makes the same of the vertices at the same position on the opposite sides. The is the ratio of the distance from the vertice to the starting point to the side length. At last, the polygon triangulation method is used to model the newly formed sides and old sides, as shown in

**Figure 9.** Schematic diagram of the Coons surface interpolation. **Figure 9.** Schematic diagram of the Coons surface interpolation.

Based on Equation (9), the *x*, *y*, and *z* of each vertex can be calculated respectively by harmonic function, according to the coordinates of four corners and four relevant points, as well as the *t* of the four relevant points. Because the synchronous forward method is adopted, the *t* values of the opposite sides are the same.

Set the *t* of the first group sides *u*<sup>0</sup> and *u*<sup>1</sup> as {*T<sup>i</sup>* |0 ≤ *i* < *m*, 0 ≤ *T<sup>i</sup>* ≤ 1}, and the *t* of the second group sides *v*<sup>0</sup> and *v*<sup>1</sup> as {*t<sup>i</sup>* |0 ≤ *i* < *n*, 0 ≤ *t<sup>i</sup>* ≤ 1}. Then, set the point interpolated by the *i* point on the *u*<sup>0</sup> and the *j* point on the *v*<sup>0</sup> as *Aij*, whose calculation formulas of coordinates are:

*x* = *u*0(*i*)*<sup>x</sup>* × *F*<sup>0</sup> *tj* +*v*0(*j*)*<sup>x</sup>* × *F*0(*Ti*) + *u*1(*i*)*<sup>x</sup>* × *F*<sup>1</sup> *tj* + *v*1(*j*)*<sup>x</sup>* × *F*1(*Ti*) −(*u*0(0)*<sup>x</sup>* × *F*0(*Ti*) + *u*0(*m* − 1)*<sup>x</sup>* × *F*1(*Ti*)) × *F*<sup>0</sup> *tj* +(*u*1(0)*<sup>x</sup>* × *F*0(*Ti*) + *u*1(*m* − 1)*<sup>x</sup>* × *F*1(*Ti*)) × *F*<sup>1</sup> *tj* (11) *y* = *u*0(*i*)*<sup>y</sup>* × *F*<sup>0</sup> *tj* +*v*0(*j*)*<sup>y</sup>* × *F*0(*Ti*) + *u*1(*i*)*<sup>y</sup>* × *F*<sup>1</sup> *tj* + *v*1(*j*)*<sup>y</sup>* × *F*1(*Ti*) −(*u*0(0)*<sup>y</sup>* × *F*0(*Ti*) + *u*0(*m* − 1)*<sup>y</sup>* × *F*1(*Ti*)) × *F*<sup>0</sup> *tj* +(*u*1(0)*<sup>y</sup>* × *F*0(*Ti*) + *u*1(*m* − 1)*<sup>y</sup>* × *F*1(*Ti*)) × *F*<sup>1</sup> *tj* (12) *z* = *u*0(*i*)*<sup>z</sup>* × *F*<sup>0</sup> *tj* +*v*0(*j*)*<sup>z</sup>* × *F*0(*Ti*) + *u*1(*i*)*<sup>z</sup>* × *F*<sup>1</sup> *tj* + *v*1(*j*)*<sup>z</sup>* × *F*1(*Ti*) −(*u*0(0)*<sup>z</sup>* × *F*0(*Ti*) + *u*0(*m* − 1)*<sup>z</sup>* × *F*1(*Ti*)) × *F*<sup>0</sup> *tj* (13)

+(*u*1(0)*<sup>z</sup>* × *F*0(*Ti*) + *u*1(*m* − 1)*<sup>z</sup>* × *F*1(*Ti*)) × *F*<sup>1</sup>

 *tj*  、

(0)௬

、

(0)௫

、、

、

; ()௫

respectively,

where the harmonic function is: (0)௭ refers to respectively on the ; (−1)௫(−1)௬(−1)௭ refer to, respectively,

൜

、

*Minerals* **2022**, *12*, x FOR PEER REVIEW 12 of 20

、

method is adopted, the values of the opposite sides are the same.

by harmonic function, according to the coordinates of four corners and four relevant points, as well as the of the four relevant points. Because the synchronous forward

 of the second group sides and ଵ as ሼ|0 ≤ < , 0 ≤ ≤ 1ሽ. Then, set the point interpolated by the point on the and the point on the as , whose

> − ((0)௫ × () + (−1)௫ × ଵ())×൫൯ + (ଵ(0)௫ × () + ଵ(−1)௫ × ଵ())×ଵ൫൯

> − ((0)௬ × () + (−1)௬ × ଵ())×൫൯ + (ଵ(0)௬ × () + ଵ(−1)௬ × ଵ())×ଵ൫൯

> − ((0)௭ × () + (−1)௭ × ଵ())×൫൯ + (ଵ(0)௭ × () + ଵ(−1)௭ × ଵ())×ଵ൫൯

> > 、、

=()௫ × ൫൯+()௫ × () + ଵ()௫ × ଵ൫൯+ଵ()௫ × ଵ()

=()௬ × ൫൯+()௬ × () + ଵ()௬ × ଵ൫൯+ଵ()௬ × ଵ()

=()௭ × ൫൯+()௭ × () + ଵ()௭ × ଵ൫൯+ଵ()௭ × ଵ()

() = 2<sup>ଷ</sup> − 3<sup>ଶ</sup> + 1

、

Set the of the first group sides and ଵ as ሼ|0 ≤ < , 0 ≤ ≤ 1ሽ, and the

Based on Equation (9), the

calculation formulas of coordinates are:

where the harmonic function is:

、

()௬

、、

、

$$\begin{cases} F\_0(t) = 2t^3 - 3t^2 + 1 \\ F\_1(t) = -2t^3 + 3t^2 \end{cases} \tag{14}$$

coordinates of the first point

、

(11)

(12)

(13)

、、

, and of each vertex can be calculated respectively

*u*0(0)*x*, *u*0(0)*y*, *u*0(0)*<sup>z</sup>* refers to respectively *x*, *y*, *z* coordinates of the first point on the *u*0; *u*0(*m* − 1)*<sup>x</sup>* , *u*0(*m* − 1)*<sup>y</sup>* , *u*0(*m* − 1)*<sup>z</sup>* refer to, respectively, *x*, *y*, *z* coordinates of the last point on the *u*0; *u*1(0)*x*, *u*1(0)*y*, *u*1(0)*<sup>z</sup>* refers to respectively *x*, *y*, *z* coordinates of the first point on the *u*1; *u*1(*m* − 1)*<sup>x</sup>* , *u*1(*m* − 1)*<sup>y</sup>* , *u*1(*m* − 1)*<sup>z</sup>* refer to, respectively, *x*, *y*, *z* coordinates of the last point on the *u*1; *u*0(*i*)*x*, *u*0(*i*)*y*, *u*0(*i*)*<sup>z</sup>* refer to, respectively, *x*, *y*, *z* coordinates of the point on the *u*0, whose *t* is *T<sup>i</sup>* ; *v*0(*j*)*x*, *v*0(*j*)*y*, *v*0(*j*)*<sup>z</sup>* refer to, respectively, *x*, *y*, *z* coordinates of the point on the *v*0, whose *t* is *t<sup>j</sup>* ;*u*1(*i*)*x*, *u*1(*i*)*y*, *u*1(*i*)*<sup>z</sup>* refer to, respectively, *x*, *y*, *z* coordinates of the point on the *u*1, whose *t* is *T<sup>i</sup>* ; *v*1(*j*)*x*, *v*1(*j*)*y*, *v*1(*j*)*<sup>z</sup>* refer to, respectively, *x*, *y*, *z* coordinates of the point on the *v*1, whose *t* is *t<sup>j</sup>* . refer to, respectively, 、、 coordinates of the last point on the ଵ; ()௫、()௬ ()௭ refer to, respectively, 、、 coordinates of the point on the , whose is ()௭ refer to, respectively, 、、 coordinates of the point on the , whose is ; ଵ()௫、ଵ()௬、ଵ()௭ refer to, respectively, 、、 coordinates of the point on the ଵ , whose is ; ଵ()௫、ଵ()௬、ଵ()௭ refer to, coordinates of the point on the ଵ, whose is .

ଵ() = −2<sup>ଷ</sup> + 3ଶ (14)

、

Performing the interpolation calculation according to Equations (10)–(13), the modeling of this four-sided region can be completed, as shown in Figure 10. Performing the interpolation calculation according to Equations (10)–(13), the modeling of this four-sided region can be completed, as shown in Figure 10.

**Figure 10.** Schematic diagram of the four-sided region modeling:(**a**) the preprocessed four-sided region; (**b**,**c**) the process and the result of the Coons surface interpolation and modeling. **Figure 10.** Schematic diagram of the four-sided region modeling:(**a**) the preprocessed four-sided region; (**b**,**c**) the process and the result of the Coons surface interpolation and modeling.

#### *3.5. Combine Sub-Meshes*

*3.5. Combine Sub-Meshes*  Multiple polygon data sets are set into one complete data set using the mature polygon data merging technology, which means all the sub-meshes are merged into a complete triangular mesh model. In the process of merging, all the polygon data will be extracted, but the extraction and expansion of point and cell attributes (scalar, vector, and normal) can only be carried out when multiple data sets contain them. As shown in Figure 11, we can obtain a complete triangular mesh model.

> We use C++ to implement the relevant algorithms of the modeling method based on the Coons surface interpolation and test them on a 64-bit computer with Windows 10 professional edition. This computer has 3.70 Ghz Intel (R) Core (TM) i7-8700u, 32GB ram, and NVIDIA GeForce RTX 2080 GPU.

model.

*4.1. Examples* 

Multiple polygon data sets are set into one complete data set using the mature polygon data merging technology, which means all the sub-meshes are merged into a complete triangular mesh model. In the process of merging, all the polygon data will be extracted, but the extraction and expansion of point and cell attributes (scalar, vector, and

**Figure 11.** Schematic diagram of the sub-mesh combination: (**a**) the sub-meshes; (**b**–**d**) the complete **Figure 11.** Schematic diagram of the sub-mesh combination: (**a**) the sub-meshes; (**b**–**d**) the complete model.
