1.1.3. Orebody Modeling

The 3D modeling of the orebody is to interpret the topological information and determine the geometric shape of the orebody using mine geological exploration data and production exploration data [41]. According to the different principles of model construction, it can be divided into explicit modeling and implicit modeling [42]. In explicit modeling, first, the contour polylines of the orebody should be delineated manually, and then the 3D reconstruction of two-dimensional contour polylines can be realized by the contour stitching method [37]. Because explicit modeling requires a great deal of manual interaction, the limitation of low efficiency of this method is gradually highlighted for orebodies with a large number of sections and complex shapes. In addition, when the modeling data change locally, the modelers need to reinterpret the data, delineate the orebody, and model through contour stitching. This process is complex, which is not conducive to the dynamic updating of the model [43]. In contrast with explicit modeling, implicit modeling is based on the implicit function. The relationships of the spatial sampling data are used to construct an implicit function, and the 3D model can be constructed by the spatial interpolation algorithm [18]. This method not only greatly reduces the manual interaction but also realizes the dynamic editing of the model, which also improves the automation and efficiency of modeling [44].

#### *1.2. Solution Strategy*

In this paper, we consider transforming the problem of orebody modeling into a modeling problem based on the Coons surface interpolation of n-sided regions. The proposed method requires the contour polylines to be in an approximate plane within the tolerance range. Firstly, we creatively propose a closed-loop division method for contour polylines interpreted from data obtained in mine production and exploration. In this method, the approximate plane of each contour polyline is used to cut all polylines, and the polylines in the final subspace are connected and grouped within the tolerance range to form closed loops. Secondly, we preprocess the formed n-sided region. For the simple single-sided region, the method based on constrained Delaunay triangulation is used for modeling. In terms of complex single-sided regions and non-quadrilateral regions, they are transformed into four-sided regions by different methods. Finally, the Coons surface is used for interpolation and modeling based on the processed four-sided regions to realize the 3D visualization of the orebody. In summary, we creatively propose a closed-loop division and modeling method of orebody contour polylines.

### **2. Overview of the Method**

The basic idea of this modeling method using Coons surface interpolation is as follows. Through the steps of n-sided region processing and function solving, the Coons surface interpolation is used to construct sub-meshes through the formed closed loops, which are obtained by interpreted cross-contour polylines. Finally, the sub-meshes will be combined to realize the modeling of the complete orebody.

The proposed method mainly consists of the following three steps. Figure 1 shows the overall process of the method.


**Figure 1.** Overall flow chart. **Figure 1.** Overall flow chart.

#### 1. Use approximate planes of each contour polyline to cut all polylines for the closed-**3. Method**

derived:

[10]:

(, ) = −[−1 () ଵ() () ଵ()] ×

#### loop division. *3.1. Coons Surface*

2. Preprocess the formed loops and separately model the processed four-sided regions through Coons surface interpolation. 3. Combine all the sub-meshes to construct a complete orebody model. **3. Method**  *3.1. Coons Surface*  With the continuous development of computer technology and modern industry, the application scope of Coons surfaces has gradually expanded from the shape design of cars, ships, and aircraft to various fields such as architectural design, medical research, and geological modeling. Coons surfaces can be divided into three categories according to boundary control conditions. The boundary control conditions of the first kind of Coons With the continuous development of computer technology and modern industry, the application scope of Coons surfaces has gradually expanded from the shape design of cars, ships, and aircraft to various fields such as architectural design, medical research, and geological modeling. Coons surfaces can be divided into three categories according to boundary control conditions. The boundary control conditions of the first kind of Coons surface only contain four boundary curves or part of them. The boundary control conditions of the second kind of Coons surface include both boundary curves and boundary tangent vectors. For the third kind of Coons surface, the boundary control conditions include boundary curves and boundary tangent vectors, as well as boundary second-order derivative vectors. The bicubic Coons surface used in interpolation in this paper belongs to the second kind of Coons surface, which can ensure the interpolated patches are continuous both in position and slope.

surface only contain four boundary curves or part of them. The boundary control conditions of the second kind of Coons surface include both boundary curves and boundary tangent vectors. For the third kind of Coons surface, the boundary control conditions include boundary curves and boundary tangent vectors, as well as boundary To satisfy the slope continuity and the position continuity at the four boundaries, the four boundaries *P*(*u*, *j*) and *P*(*i*, *v*) and their cross-border tangent vectors *Pv*(*u*, *j*) and *Pu*(*i*, *v*) must be given at the same time, as shown in Figure 2. *Minerals* **2022**, *12*, x FOR PEER REVIEW 6 of 20

**Figure 2.** Schematic diagram of the interpolation of the second kind Coons surface: (**a**) boundaries and their cross-border tangent vectors of the surface patch; (**b**) the interpolated surface. **Figure 2.** Schematic diagram of the interpolation of the second kind Coons surface: (**a**) boundaries and their cross-border tangent vectors of the surface patch; (**b**) the interpolated surface.

Under these constraints, the second kind of Coons surface equation [10] can be

This fifth-order matrix is the boundary information matrix of the second kind of Coons patches (Coons patches with given boundaries and cross-border tangent vectors). , ଵ, and ଵ are mixed functions, which meet the following conditions

> () = ൜1, = 0, ≠

The above methods and formulas are extremely rigorous in theory, but too much boundary information is needed, which is not convenient for calculation and application. To meet its application requirements and simplify the calculation process as much as possible, Professor Coons proposed using corner information and mixed functions to

If the point vectors, tangent vectors of direction, tangent vectors of direction, and mixed partial derivative vectors of the four corners are known, as shown in Figure 3, by applying the mixed function , ଵ, and ଵ, we can express the four boundaries

() = () = 0, , = 0,1

⎧0() = 2<sup>3</sup> − 3<sup>2</sup> + 1 1() = −2<sup>3</sup> + 3<sup>2</sup> 0() = <sup>3</sup> − 2<sup>2</sup> + 1() = <sup>3</sup> − <sup>2</sup>

௩(, 0) ௩(, 1) ௩(0,0) ௩(0,1) ௩(1,0) ௩(1,1)

௨௩(0,0) ௨௩(0,1) ௨௩(1,0) ௨௩(1,1)⎦

⎥ ⎥ ⎥ ⎤ ⎣ ⎢ ⎢ ⎢ <sup>⎡</sup> −1 () ଵ() () ଵ()⎦ ⎥ ⎥ ⎥ ⎤

, , ∈ [0,1] (1)

(2)

(3)

(1, ) (1,0) (1,1)

௨(0, ) ௨(0,0) ௨(0,1) ௨(1, ) ௨(1,0) ௨(1,1)

> () = ᇱ

 ᇱ

Based on Equation (2), we can solve that [10]:

⎩ ⎪ ⎨ ⎪

define the boundary curves and their cross-border tangent vectors.

of the surface patch and their cross-border tangent vectors as follows [10].

⎣ ⎢ ⎢ ⎢ Under these constraints, the second kind of Coons surface equation [10] can be derived:

$$Q(u,v) = -\begin{bmatrix} -1 & F\_0(u)F\_1(u)G\_0(u)G\_1(u) \end{bmatrix} \times \begin{bmatrix} 0 & P(u,0) & P(u,1) & P\_0(u,0) & P\_0(u,1) \\ P(0,v) & P(0,0) & P(0,1) & P\_0(0,0) & P\_0(0,1) \\ P(1,v) & P(1,0) & P(1,1) & P\_0(1,0) & P\_0(1,1) \\ P\_\text{in}(0,v) & P\_\text{in}(0,0) & P\_\text{in}(0,1) & P\_\text{in}(0,0) & P\_\text{in}(0,1) \\ P\_\text{in}(1,v) & P\_\text{in}(1,0) & P\_\text{in}(1,1) & P\_\text{in}(1,0) & P\_\text{in}(1,1) \end{bmatrix} \begin{bmatrix} -1 \\ F\_0(v) \\ F\_1(v) \\ G\_0(v) \\ G\_1(v) \end{bmatrix}, \mu, v \in [0,1] \tag{1}$$

This fifth-order matrix is the boundary information matrix of the second kind of Coons patches (Coons patches with given boundaries and cross-border tangent vectors). *F*0, *F*1, *G*<sup>0</sup> and *G*<sup>1</sup> are mixed functions, which meet the following conditions [10]:

$$\begin{aligned} F\_i(j) &= G\_i'(j) = \begin{cases} 1, i = j \\ 0, i \neq j \end{cases} \\ F\_i'(j) &= G\_i(j) = 0, i, j = 0, 1 \end{aligned} \tag{2}$$

Based on Equation (2), we can solve that [10]:

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

The above methods and formulas are extremely rigorous in theory, but too much boundary information is needed, which is not convenient for calculation and application. To meet its application requirements and simplify the calculation process as much as possible, Professor Coons proposed using corner information and mixed functions to define the boundary curves and their cross-border tangent vectors.

If the point vectors, tangent vectors of *u* direction, tangent vectors of *v* direction, and mixed partial derivative vectors of the four corners are known, as shown in Figure 3, by applying the mixed function *F*0, *F*1, *G*<sup>0</sup> and *G*1, we can express the four boundaries of the surface patch and their cross-border tangent vectors as follows [10].

$$\begin{cases} P(0,v) = F\_0(v)P(0,0) + F\_1(v)P(0,1) + G\_0(v)P\_v(0,0) + G\_1(v)P\_v(0,1) \\ P(1,v) = F\_0(v)P(1,0) + F\_1(v)P(1,1) + G\_0(v)P\_v(1,0) + G\_1(v)P\_v(1,1) \\ P(u,0) = F\_0(u)P(0,0) + F\_1(u)P(0,1) + G\_0(u)P\_u(0,0) + G\_1(u)P\_u(0,1) \\ P(u,1) = F\_0(u)P(1,0) + F\_1(u)P(1,1) + G\_0(u)P\_u(1,0) + G\_1(u)P\_u(1,1) \\ P\_u(u,0) = F\_0(u)P\_u(0,0) + F\_1(u)P\_u(0,1) + G\_0(u)P\_{uv}(0,0) + G\_1(u)P\_{uv}(0,1) \\ P\_u(u,1) = F\_0(u)P\_u(1,0) + F\_1(u)P\_u(1,1) + G\_0(u)P\_{uv}(1,0) + G\_1(u)P\_{uv}(1,1) \\ P\_v(u,0) = F\_0(u)P\_v(0,0) + F\_1(u)P\_v(0,1) + G\_0(u)P\_{uv}(0,0) + G\_1(u)P\_{uv}(0,1) \\ P\_v(u,1) = F\_0(u)P\_v(1,0) + F\_1(u)P\_v(1,1) + G\_0(u)P\_{uv}(1,0) + G\_1(u)P\_{uv}(1,1) \end{cases} \tag{4}$$

⎩ ⎨ ⎧

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

**Figure 3.** Schematic diagram of the interpolation of the bicubic Coons surface: (**a**) corner information, boundaries, and their cross-border tangent vectors of the surface patch; (**b**) the interpolated surface. **Figure 3.** Schematic diagram of the interpolation of the bicubic Coons surface: (**a**) corner information, boundaries, and their cross-border tangent vectors of the surface patch; (**b**) the interpolated surface.

⎪ ⎪ <sup>⎧</sup> (0, ) = 0()(0,0) + 1()(0,1) + 0()(0,0) + 1()(0,1) (1, ) = 0()(1,0) + 1()(1,1) + 0()(1,0) + 1()(1,1) (, 0) = 0()(0,0) + 1()(0,1) + 0()(0,0) + 1()(0,1) After bringing Equation (2) into Equation (1) and simplifying it, we can obtain the simplified surface equation [10]:

ଵ()

$$\begin{cases} Q(u,v) = \left[F\_{\mathbb{D}}(u)F\_{\mathbb{1}}(u)G\_{\mathbb{D}}(u)G\_{\mathbb{1}}(u)\right] \times \begin{bmatrix} P(0,0) & P(0,1) & P\_{\mathbb{D}}(0,0) & P\_{\mathbb{D}}(0,1) \\ P(1,0) & P(1,1) & P\_{\mathbb{D}}(1,0) & P\_{\mathbb{D}}(1,1) \\ P\_{\mathbb{u}}(0,0) & P\_{\mathbb{u}}(0,1) & P\_{\mathbb{u}\mathbb{D}}(0,0) & P\_{\mathbb{u}\mathbb{U}}(0,1) \\ P\_{\mathbb{u}}(1,0) & P\_{\mathbb{u}}(1,1) & P\_{\mathbb{u}\mathbb{D}}(1,0) & P\_{\mathbb{U}\mathbb{U}}(1,1) \end{bmatrix} \begin{bmatrix} F\_{0}(v) \\ F\_{1}(v) \\ G\_{0}(v) \\ G\_{1}(v) \end{bmatrix}, u,v \in [0,1] \tag{5}$$

simplified surface equation [10]: where *F*0, *F*1, *G*<sup>0</sup> and *G*<sup>1</sup> are the same as Equation (3).

> ௨௩(0,0) ௨௩(1,0)

(1,0) (1,1) ௩(1,0) ௩(1,1)

−3

௨(0,1) ௨(1,1)

⎡ (0,0) (0,1) ௩(0,0) ௩(0,1) ⎤ ⎡ () ⎤ The cubic mixed function and Equation (3) can be written as [10]:

$$Q(u,v) = \mathcal{U}MBM^TV^T\tag{6}$$

where , ଵ, and ଵ are the same as Equation (3). where:

where:

⎣ ⎢ ⎢

௨(0,0) ௨(1,0)

=

⎢ ⎢ ⎡

(, ) = [() ଵ() () ଵ()] ×

⎪

$$\mathcal{U} = \begin{bmatrix} u^3 & u^2 & u & 1 \end{bmatrix}, \; V = \begin{bmatrix} v^3 & v^2 & v & 1 \end{bmatrix}, \mu, v \in [0, 1] \tag{7}$$

$$M = \begin{bmatrix} 2 & -2 & 1 & 1 \\ -3 & 3 & -2 & -1 \\ 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 \end{bmatrix} \tag{8}$$

$$B = \begin{bmatrix} P(0,0) & P(0,1) & P\_v(0,0) & P\_v(0,1) \\ P(1,0) & P(1,1) & P\_v(1,0) & P\_v(1,1) \\ P\_u(0,0) & P\_u(0,1) & P\_{uv}(0,0) & P\_{uv}(0,1) \\ P\_u(1,0) & P\_u(1,1) & P\_{uv}(1,0) & P\_{uv}(1,1) \end{bmatrix} \tag{9}$$

⎣ ௨(1,0) ௨(1,1) ௨௩(1,0) ௨௩(1,1)⎦ The parametric equation of the Coons surface is [10]: The parametric equation of the Coons surface is [10]:

$$\begin{cases} x(\mu, v) = \mathcal{U}MB\_{\mathcal{X}}M^TV^T\\ y(\mu, v) = \mathcal{U}MB\_{\mathcal{Y}}M^TV^T\\ z(\mu, v) = \mathcal{U}MB\_{\mathcal{Z}}M^TV^T \end{cases} \tag{10}$$

where *Bx*, *By*, and *B<sup>z</sup>* are the coordinate components of boundary information matrix *B*.

#### *3.2. Closed-Loop Division*

Cut all contour polylines using the approximate planes where each contour polyline is located. The contour polyline on the positive side (the positive and negative sides can be specified arbitrarily without affecting the final dividing result) is coded as 1 in this plane bit. The contour polyline on the negative side is coded as 2 in this plane bit. The contour polyline on the approximate plane is coded as 0 because it belongs to both positive and negative sides, and it will be expanded in the subsequent steps. Each time an approximate plane is added for cutting, the length of each polyline mark string is increased by 1, whose value is 0, 1, or 2.

、、

①

②

As shown in Figure 4, there are six contour polylines in the figure, of which the bottom two are coplanar. Cut the contour polylines with the five approximate planes *α*, *β*, *γ*, *δ*,*ε* where the contour polylines are located and divide them into four subspaces <sup>1</sup> , <sup>2</sup> , <sup>3</sup> , <sup>4</sup> . The direction indicated by the arrow is the positive direction of the plane, and the tag array length of each contour polyline is 5. As shown in Figure 4, there are six contour polylines in the figure, of which the bottom two are coplanar. Cut the contour polylines with the five approximate planes 、、 where the contour polylines are located and divide them into four subspaces ①、②、③、④. The direction indicated by the arrow is the positive direction of the plane, and the tag array length of each contour polyline is 5. As shown in Figure 4, there are six contour polylines in the figure, of which the bottom two are coplanar. Cut the contour polylines with the five approximate planes 、、、、 where the contour polylines are located and divide them into four subspaces ①、②、③、④. The direction indicated by the arrow is the positive direction of the plane, and the tag array length of each contour polyline is 5.

(10)

(10)

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

ቐ

*3.2. Closed-Loop Division* 

increased by 1, whose value is 0, 1, or 2.

increased by 1, whose value is 0, 1, or 2.

*3.2. Closed-Loop Division* 

(, ) = (, ) = (, ) =

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

ቐ

where ௫, ௬, and ௭ are the coordinate components of boundary information matrix .

Cut all contour polylines using the approximate planes where each contour polyline is located. The contour polyline on the positive side (the positive and negative sides can be specified arbitrarily without affecting the final dividing result) is coded as 1 in this plane bit. The contour polyline on the negative side is coded as 2 in this plane bit. The contour polyline on the approximate plane is coded as 0 because it belongs to both positive and negative sides, and it will be expanded in the subsequent steps. Each time an approximate plane is added for cutting, the length of each polyline mark string is

(, ) = (, ) = (, ) =

where ௫, ௬, and ௭ are the coordinate components of boundary information matrix .

Cut all contour polylines using the approximate planes where each contour polyline is located. The contour polyline on the positive side (the positive and negative sides can be specified arbitrarily without affecting the final dividing result) is coded as 1 in this plane bit. The contour polyline on the negative side is coded as 2 in this plane bit. The contour polyline on the approximate plane is coded as 0 because it belongs to both positive and negative sides, and it will be expanded in the subsequent steps. Each time an approximate plane is added for cutting, the length of each polyline mark string is

**Figure 4.** Schematic diagram of the closed-loop division: (**a**) the example model and the crosscontour polylines of it; (**b**) the formed subspace after using approximate planes to cut polylines. **Figure 4.** Schematic diagram of the closed-loop division: (**a**) the example model and the cross-contour polylines of it; (**b**) the formed subspace after using approximate planes to cut polylines. contour polylines of it; (**b**) the formed subspace after using approximate planes to cut polylines. As shown in Figures 4 and 5, the contour polylines , , , are cut by the

As shown in Figures 4 and 5, the contour polylines , , , are cut by the approximate planes and coded according to the above rules. The initial coding result is shown in Table 1. As shown in Figures 4 and 5, the contour polylines *a*, *b*, *c*, *d* are cut by the approximate planes and coded according to the above rules. The initial coding result is shown in Table 1. approximate planes and coded according to the above rules. The initial coding result is shown in Table 1.

**Figure 5.** The top view of formed subspaces: (**a**) the top view of the complete subspaces; (**b**) the top view of the divided contour polylines. **Figure 5.** The top view of formed subspaces: (**a**) the top view of the complete subspaces; (**b**) the top view of the divided contour polylines. **Figure 5.** The top view of formed subspaces: (**a**) the top view of the complete subspaces; (**b**) the top view of the divided contour polylines.


**Table 1.** Initial coding table after cutting polylines with the approximate planes.

Copy the contour polyline containing 0 in the code array. On the plane bit where 0 is located, two contour polylines are coded as 1 and 2, and the other codes remain unchanged. For example, expand *a*: 10,211 to *a*1: 11,211 and *a*2: 12,211. After processing all the contour polylines, the final coding result is obtained as shown in Table 2.

The contour polylines with the same code array are the contour polylines in the same subspace. As shown in Table 2, *a*1, *b*1, *c*2, *d*<sup>1</sup> marked by red have the same code arrays, which means that they are in the same subspace. Then, in the same subspace, the contour polylines within the specified tolerance are connected to form the closed loops.


**Table 2.** Final coding table after extending the side whose code array contains 0.

It is worth noting that if the topological adjacency relationship between some contour polylines is too complex, the above automatic closed-loop division algorithm may not be able to obtain the required results. At this time, it is necessary to manually process and group the contour polylines to complete the closed-loop division.
