Orebody Modeling Method Based on the Coons Surface Interpolation

Abstract

This paper presents a surface modeling method for interpolating orebody models based on a set of cross-contour polylines (geological polylines interpreted from the raw geological sampling data) using the bi-Coons surface interpolation method. The method is particularly applicable to geological data with cross-contour polylines acquired during the geological and exploration processes. The innovation of this paper is that the proposed method can automatically divide the closed loops and automatically combine the sub-meshes. The method solves the problem that it is difficult to divide closed loops from the cross-contour polylines with complex shapes, and it greatly improves the efficiency of modeling based on complex cross-contour polylines. It consists of three stages: (1) Divide closed loops using approximate planes of contour polylines; each loop is viewed as a polygon combined with several polylines, that is the n-sided region. (2) After processing the formed n-sided regions, Coons surface interpolation is improved to complete the modeling of every single loop (3) Combine all sub-meshes to form a complete orebody model. The corresponding algorithm was implemented using the C++ programing language on 3D modeling software. Experimental results show that the proposed orebody modeling method is useful for efficiently recovering complex orebody models from a set of cross-contour polylines.

1. Introduction

Three-dimensional geological modeling and visualization is one of the hot topics in the field of geosciences [1]. As a geological body with high economic value, orebody modeling is also one of the current research hot topics and is of great significance for promoting mine digitization and improving mining efficiency [2]. In 2002, Mallet [3] proposed a new method called discrete smooth interpolation (DSI), based on which a number of orebody modeling works have been published [4,5,6,7]. Additionally, some commercial software [8] has been developed using different technologies including triangulate surfaces, radial functions, and parametric surfaces [9]. However, due to the complex shape of the orebody, much research work still needs to be done in orebody modeling. Our research is dedicated to developing the software for modeling orebody using small triangles based on the Coons surface [10].

We can often obtain the cross-contour polylines of planes and sections of orebodies through geological logging in the geological exploration and production exploration of mines, especially of some precious metal mines with flat orebody shapes. These contour polylines are interlaced with complex shapes. Therefore, if the explicit modeling method based on the contour stitching method is used, the model will be constructed with poor quality and low efficiency [11]. Modeling orebodies with high speed and great quality is of great significance to improving the production efficiency of mines [12]. The orebody modeling method we propose based on Coons surface interpolation can automatically divide the contour polylines into closed loops and interpolate and model the formed n-sided regions; then, the constructed sub-meshes will be combined without manual intervention. Based on the above steps, the complete orebody model can be quickly constructed.

1.1. Related Research Work

Based on the above analysis, we transform the 3D modeling of the orebody into the Coons surface interpolation of the interpreted contour polylines. Therefore, we will briefly summarize the research work on free-form surface reconstruction, contour interpolation modeling, and orebody modeling.

1.1.1. Free-Form Surface Reconstruction

Free-form surface reconstruction is a key technology of orebody 3D modeling [13,14,15]. In practical application, the commonly used surface reconstruction methods include implicit function surface, triangular Bezier surface, B-spline surface, NURBS surface, Coons surface, polygon model, and so on.

The implicit function surface [16,17,18] refers to the surface represented by equation $F\left(x,y,z\right)=0$, which can represent quadric surfaces such as spherical surfaces, cylindrical surfaces, and more complex surfaces. When reconstructing the surface, the implicit function equations are constructed using the existing data, and then the surface model can be generated by obtaining isosurfaces through the Marching Cube method or other methods. Though using the implicit function to reconstruct a surface has the advantages of small calculation and easy solution, there are also some problems such as unclear geometric meaning. In the future, a great amount of research work needs to be carried out to make extensive use of implicit modeling.

In 1982, Farin [19] proposed the method of constructing the triangular Bezier surface based on the triangulation of scattered data. This surface is explicit and $C^1$ continuous with flexible construction. Using the triangular Bezier surface for surface reconstruction requires less manual intervention and is easy to automate. However, it also has some disadvantages, such as the insufficient ability for surface modification and poor controllability.

The B-spline curve and surface [20,21,22] have many excellent properties, such as geometric invariance, convex hull, convexity preservation, and local support, which lead to the powerful function of representing and designing free-form curves and surfaces. It is a commonly used free-form surface suitable for engineering shape design, but it cannot be well applied to represent and design elementary curves and surfaces. Therefore, researchers proposed the NURBS surface method [23,24,25] for surface reconstruction, which both has the powerful function of representing and designing free-form curves and surfaces like B-spline and can accurately represent quadric arcs and quadric surfaces. The unification of analytical geometry and free-form curves and surfaces is well realized by this method. Due to the introduction of a weight factor and a manipulation control vertex, it provides sufficient flexibility for various shape designs. However, the parameterization effect will become very poor with an inappropriate weight factor. In addition, there are still many problems to be further solved in the calculation of NURBS surface intersection.

In 1964, Professor S. A. Coons of the Massachusetts Institute of Technology (MIT) proposed a general method to describe surfaces, the Coons surfaces method [10], which has unique advantages in engineering applications. Given the boundary conditions, Coons patches of the corresponding degree can be constructed. Therefore, in theory, the Coons method can be used to construct arbitrary order patches through boundary constraints [26,27], which makes Coons surfaces show strong advantages in representing and designing surfaces. However, in the Coons surface method, it is difficult to control the surface shape. When the boundary curve is determined, the shape can only be modified by changing the torsion vector [28].

The most commonly used polygon model method is the triangle model method based on triangulation. The triangulation technology in the plane has matured [29], and great research progress has been made in the triangulation of surface data and tetrahedrons of spatial scattered data [30,31]. The characteristics of the polygon model are: (1) adjacent surface patches can only be $C^0$ continuous and (2) data amounts are large. When there needs high accuracy of a model, a large number of polygons need to be used to represent the surface approximately, which has a huge amount of data. As a result, this method is not suitable for constructing complex models with high accuracy. Therefore, in this paper, the polygon model method is used only in the modeling of simple single-sided regions.

1.1.2. Contour Interpolation Modeling

Visualization based on contour interpolation refers to the modeling of three-dimensional entities using a sequence of two-dimensional contours, which is an important research direction in scientific computing visualization. According to the number of contour polylines of two adjacent layers, it can be divided into single contour modeling and multi-contour modeling.

Single contour modeling refers to the three-dimensional model reconstruction of the contour polylines when there is only one contour polyline in each layer on two adjacent planes, which is relatively simple. Researchers have proposed many optimization methods for single contour modeling, such as the minimum surface product method [32], the maximum volume method [33] based on the global search strategy, the shortest diagonal method [34], and the synchronous advance method of adjacent contours [35] based on local calculation and decision.

Multi-contour modeling refers to 3D model reconstruction when there are multiple contour polylines on two adjacent planes, which difficulty is the correspondence and branch processing between contour polylines. In terms of the correspondence of contour polylines, Meyers et al. [36] proposed the minimum spanning tree method, but this method still has some limitations. In some special cases, the contour polylines cannot correctly correspond, and subsequent manual processing is required. For the branching problem between contour lines, Ekoule et al. [37] proposed a method for constructing the middle contour polyline in the middle of the two adjacent planes of contour polylines. This method transforms the branching problem into the connection problem between a series of single contour polylines. However, in some cases, this method will leave holes on the constructed model surface.

Jones et al. [38] proposed an isosurface construction method using volume data, which requires that contour polylines have simple shapes and be completely closed. There is no need to judge the corresponding relationship and branch relationship when using this method for surface reconstruction, and it is suitable for both convex and non-convex contour polylines. However, it also has the disadvantage of a large amount of calculation. Additionally, Zhong et al. [39] proposed a method of reconstructing the 3D orebody by cross-section contour polylines, which allows for adding geometric constraints and can easily control the shape of the model, but many discontinuous parts may be generated after reconstruction if the data are relatively sparse. Moreover, Wu et al. [40] proposed a 3D orebody modeling method based on the normal estimation of cross-contour polylines. First, the normals of the cross-contour polylines will be estimated, based on which the radial basis function will be used for interpolating and modeling. This method improves the modeling automation of contour polylines, but it requires strict intersection between contour polylines, and they can be completely broken at intersections.

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.

• Use approximate planes of each contour polyline to cut all polylines for the closed-loop division.
• Preprocess the formed loops and separately model the processed four-sided regions through Coons surface interpolation.
• 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 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.

To satisfy the slope continuity and the position continuity at the four boundaries, the four boundaries $P\left(u,j\right)$ and $P\left(i,v\right)$ and their cross-border tangent vectors $P_v\left(u,j\right)$ and $P_u\left(i,v\right)$ must be given at the same time, as shown in Figure 2.

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

$$Q(u,v)=-\left[\begin{array}{cc}-1\hfill & F_0\left(u\right)F_1\left(u\right)G_0\left(u\right)G_1\left(u\right)\hfill \end{array}\right]\times \left[\begin{array}{ccccc}0& P(u,0)& P(u,1)& P_v(u,0)& P_v(u,1)\\ P(0,v)& P(0,0)& P(0,1)& P_v(0,0)& P_v(0,1)\\ P(1,v)& P(1,0)& P(1,1)& P_v(1,0)& P_v(1,1)\\ P_u(0,v)& P_u(0,0)& P_u(0,1)& P_{uv}(0,0)& P_{uv}(0,1)\\ P_u(1,v)& P_u(1,0)& P_u(1,1)& P_{uv}(1,0)& P_{uv}(1,1)\end{array}\right]\left[\begin{array}{c}-1\\ F_0\left(v\right)\\ F_1\left(v\right)\\ G_0\left(v\right)\\ G_1\left(v\right)\end{array}\right],u,v\in [0,1]$$

(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_0$ and $G_1$ are mixed functions, which meet the following conditions [10]:

$$\begin{gathered} F_i\left(j\right)=G_i^{\prime }\left(j\right)=\{\begin{array}{c}1,i=j\\ 0,i\ne j\end{array} \\ {F}_i^{\prime }\left(j\right)=G_i\left(j\right)=0,i,j=0,1 \end{gathered}$$

(2)

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

$$\{\begin{array}{c}{F}_0\left(u\right)=2u^3-3u^2+1\\ F_1\left(u\right)=-2u^3+3u^2\\ G_0\left(u\right)=u^3-2u^2+u\\ G_1\left(u\right)=u^3-u^2\end{array}$$

(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_0$ and $G_1$, we can express the four boundaries of the surface patch and their cross-border tangent vectors as follows [10].

$$\left\{\begin{array}{cc}\hfill P(0,v)& =F_0\left(v\right)P(0,0)+F_1\left(v\right)P(0,1)+G_0\left(v\right)P_v(0,0)+G_1\left(v\right)P_v(0,1)\hfill \\ \hfill P(1,v)& =F_0\left(v\right)P(1,0)+F_1\left(v\right)P(1,1)+G_0\left(v\right)P_v(1,0)+G_1\left(v\right)P_v(1,1)\hfill \\ \hfill P(u,0)& =F_0\left(u\right)P(0,0)+F_1\left(u\right)P(0,1)+G_0\left(u\right)P_u(0,0)+G_1\left(u\right)P_u(0,1)\hfill \\ \hfill P(u,1)& =F_0\left(u\right)P(1,0)+F_1\left(u\right)P(1,1)+G_0\left(u\right)P_u(1,0)+G_1\left(u\right)P_u(1,1)\hfill \\ \hfill P_u(u,0)& =F_0\left(u\right)P_u(0,0)+F_1\left(u\right)P_u(0,1)+G_0\left(u\right)P_{uv}(0,0)+G_1\left(u\right)P_{uv}(0,1)\hfill \\ \hfill P_u(u,1)& =F_0\left(u\right)P_u(1,0)+F_1\left(u\right)P_u(1,1)+G_0\left(u\right)P_{uv}(1,0)+G_1\left(u\right)P_{uv}(1,1)\hfill \\ \hfill P_v(u,0)& =F_0\left(u\right)P_v(0,0)+F_1\left(u\right)P_v(0,1)+G_0\left(u\right)P_{uv}(0,0)+G_1\left(u\right)P_{uv}(0,1)\hfill \\ \hfill P_v(u,1)& =F_0\left(u\right)P_v(1,0)+F_1\left(u\right)P_v(1,1)+G_0\left(u\right)P_{uv}(1,0)+G_1\left(u\right)P_{uv}(1,1)\hfill \end{array}\right.$$

(4)

After bringing Equation (2) into Equation (1) and simplifying it, we can obtain the simplified surface equation [10]:

$$\left\{Q(u,v)=\left[F_0\left(u\right)F_1\left(u\right)G_0\left(u\right)G_1\left(u\right)\right]\times \left[\begin{array}{cccc}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{array}\right]\left[\begin{array}{c}{F}_0\left(v\right)\hfill \\ F_1\left(v\right)\hfill \\ G_0\left(v\right)\hfill \\ G_1\left(v\right)\hfill \end{array}\right],u,v\in [0,1]\right.$$

(5)

where $F_0$, $F_1$, $G_0$ and $G_1$ are the same as Equation (3).

The cubic mixed function and Equation (3) can be written as [10]:

$$Q\left(u,v\right)=UMB{M}^T{V}^T$$

(6)

where:

$$U=\left[\begin{array}{ccc}{u}^3& u^2& \begin{array}{cc}u& 1\end{array}\end{array}\right], V=\left[\begin{array}{ccc}{v}^3& v^2& \begin{array}{cc}v& 1\end{array}\end{array}\right],u,v\in \left[0,1\right]$$

(7)

$$M=\left[\begin{array}{cccc}2& -2& 1& 1\\ -3& 3& -2& -1\\ 0& 0& 1& 0\\ 1& 0& 0& 0\end{array}\right]$$

(8)

$$B=\left[\begin{array}{ccc}P\left(0,0\right)& P\left(0,1\right)& \begin{array}{cc}{P}_v\left(0,0\right)& P_v\left(0,1\right)\end{array}\\ P\left(1,0\right)& P\left(1,1\right)& \begin{array}{cc}{P}_v\left(1,0\right)& P_v\left(1,1\right)\end{array}\\ \begin{array}{c}{P}_u\left(0,0\right)\\ P_u\left(1,0\right)\end{array}& \begin{array}{c}{P}_u\left(0,1\right)\\ P_u\left(1,1\right)\end{array}& \begin{array}{cc}\begin{array}{c}{P}_{uv}\left(0,0\right)\\ P_{uv}\left(1,0\right)\end{array}& \begin{array}{c}{P}_{uv}\left(0,1\right)\\ P_{uv}\left(1,1\right)\end{array}\end{array}\end{array}\right]$$

(9)

The parametric equation of the Coons surface is [10]:

$$\{\begin{array}{c}x\left(u,v\right)=UM{B}_x{M}^T{V}^T\\ y\left(u,v\right)=UM{B}_y{M}^T{V}^T\\ z\left(u,v\right)=UM{B}_z{M}^T{V}^T\end{array}$$

(10)

where $B_x, B_y$, and $B_z$ 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 $\alpha ,\beta ,\gamma ,\delta ,\epsilon$ 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 and Figure 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. Initial coding table after cutting polylines with the approximate planes.

	Plane Bit	α	β	γ	δ	ε
Segments
a		1	0	2	1	1
b		1	1	2	1	0
c		1	1	0	1	1
d		1	1	2	0	1

.

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₁: 11,211 and a₂: 12,211. After processing all the contour polylines, the final coding result is obtained as shown in

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

	Plane bit	α	β	γ	δ	ε
Segments
a1		1	1	2	1	1
a2		1	2	2	1	1
b1		1	1	2	1	1
b2		1	1	2	1	2
c1		1	1	1	1	1
c2		1	1	2	1	1
d1		1	1	2	1	1
d2		1	1	2	2	1

.

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_1$ 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.

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.

3.3. Process n-Sided Regions

Using Coons surface for interpolation calculation needs the boundary information of four sides. However, many of the closed loops divided before are non-quadrilateral. Therefore, we need to preprocess the n-sided regions before interpolation calculation.

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.

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.

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_i,i=1,2,\dots ,n$ are the same as $A$. The coordinates of $B_i,i=1,2,\dots ,n$ are the same as $B$.

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

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.

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 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.

3.4. Sided Region Modeling

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\times n$. As shown in Figure 9, set the first group of opposite sides respectively as $u_0$ and $u_1$, the second group of opposite sides respectively as $v_0$ and $v_1$.

Based on Equation (9), the $x, y, \mathrm{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_0$ and $u_1$ as $\left\{T_i|0\le i<m, 0\le T_i\le 1\right\}$, and the $t$ of the second group sides $v_0$ and $v_1$ as $\left\{t_i|0\le i<n, 0\le t_i\le 1\right\}$. Then, set the point interpolated by the $i$ point on the $u_0$ and the $j$ point on the $v_0$ as $A_{ij}$, whose calculation formulas of coordinates are:

$$\begin{array}{cc}\hfill x=u_0{\left(i\right)}_x\times F_0\left(t_j\right)& +v_0{\left(j\right)}_x\times F_0\left(T_i\right)+u_1{\left(i\right)}_x\times F_1\left(t_j\right)+v_1{\left(j\right)}_x\times F_1\left(T_i\right)\hfill \\ & -\left(u_0{\left(0\right)}_x\times F_0\left(T_i\right)+u_0{(m-1)}_x\times F_1\left(T_i\right)\right)\times F_0\left(t_j\right)\hfill \\ & +\left(u_1{\left(0\right)}_x\times F_0\left(T_i\right)+u_1{(m-1)}_x\times F_1\left(T_i\right)\right)\times F_1\left(t_j\right)\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill y=u_0{\left(i\right)}_y\times F_0\left(t_j\right)& +v_0{\left(j\right)}_y\times F_0\left(T_i\right)+u_1{\left(i\right)}_y\times F_1\left(t_j\right)+v_1{\left(j\right)}_y\times F_1\left(T_i\right)\hfill \\ & -(u_0{\left(0\right)}_y\times F_0\left(T_i\right)+u_0{(m-1)}_y\times F_1\left(T_i\right))\times F_0\left(t_j\right)\hfill \\ & +(u_1{\left(0\right)}_y\times F_0\left(T_i\right)+u_1{(m-1)}_y\times F_1\left(T_i\right))\times F_1\left(t_j\right)\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill z=u_0{\left(i\right)}_z\times F_0\left(t_j\right)& +v_0{\left(j\right)}_z\times F_0\left(T_i\right)+u_1{\left(i\right)}_z\times F_1\left(t_j\right)+v_1{\left(j\right)}_z\times F_1\left(T_i\right)\hfill \\ & -\left(u_0{\left(0\right)}_z\times F_0\left(T_i\right)+u_0{(m-1)}_z\times F_1\left(T_i\right)\right)\times F_0\left(t_j\right)\hfill \\ & +\left(u_1{\left(0\right)}_z\times F_0\left(T_i\right)+u_1{(m-1)}_z\times F_1\left(T_i\right)\right)\times F_1\left(t_j\right)\hfill \end{array}$$

(13)

where the harmonic function is:

$$\{\begin{array}{c}{F}_0\left(t\right)=2t^3-3t^2+1\\ F_1\left(t\right)=-2t^3+3t^2\end{array}$$

(14)

$u_0{\left(0\right)}_x,u_0{\left(0\right)}_y,u_0{\left(0\right)}_z$ refers to respectively $x,y,z$ coordinates of the first point on the $u_0$; $u_0{\left(m-1\right)}_x,u_0{\left(m-1\right)}_y,u_0{\left(m-1\right)}_z$ refer to, respectively, $x,y,z$ coordinates of the last point on the $u_0$; $u_1{\left(0\right)}_x,u_1{\left(0\right)}_y,u_1{\left(0\right)}_z$ refers to respectively $x,y,z$ coordinates of the first point on the $u_1$; $u_1{\left(m-1\right)}_x,u_1{\left(m-1\right)}_y,u_1{\left(m-1\right)}_z$ refer to, respectively, $x,y,z$ coordinates of the last point on the $u_1$; $u_0{\left(i\right)}_x,u_0{\left(i\right)}_y,u_0{\left(i\right)}_z$ refer to, respectively, $x,y,z$ coordinates of the point on the $u_0$, whose $t$ is T_i; $v_0{\left(j\right)}_x,v_0{\left(j\right)}_y,v_0{\left(j\right)}_z$ refer to, respectively, $x,y,z$ coordinates of the point on the $v_0$, whose $t$ is $t_j$;$u_1{\left(i\right)}_x,u_1{\left(i\right)}_y,u_1{\left(i\right)}_z$ refer to, respectively, $x,y,z$ coordinates of the point on the $u_1$, whose $t$ is $T_i$; $v_1{\left(j\right)}_x,v_1{\left(j\right)}_y,v_1{\left(j\right)}_z$ refer to, respectively, $x,y,z$ coordinates of the point on the $v_1$, whose $t$ is $t_j$.

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

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.

4. Results and Discussion

In this paper, five groups of different data are used to test the modeling method proposed in this work. These data are from the orebody of a tin mine in Yunnan, China, and are obtained in production and exploration. Based on these data, the original contour polylines can be produced through geological logging. The following experimental figures clearly show the process and the result of single closed-loop modeling and sub-mesh combination.

4.1. Examples

Figure 12 shows the processes of closed-loop division, single closed-loop modeling, and sub-mesh combination of the first group of cross-contour polylines. Figure 12a shows the original cross contour polylines of the orebody, each of which is in the same plane within the tolerance range. Then, the contour polylines are cut with the approximate planes. As shown in Figure 12b, each contour polyline is divided and grouped. After that, single closed-loop modeling is carried out on the contour polylines through the Coons surface interpolation, and the results are shown in Figure 12c. It can be seen that every closed loop has been modeled, but the sub-meshes are separate from each other. Finally, the sub-meshes are combined into a complete orebody model using polygon data merging technology, as shown in Figure 12d.

Figure 13 and Figure 14 show the modeling process on contour polylines with different shapes. It can be seen that the algorithm implemented here can still carry out closed-loop division, single closed-loop modeling, and sub-mesh combination on more complex contour lines with good effects.

Figure 15 shows the modeling process of more complex contour polylines, in which most of the contour polylines do not strictly intersect. Figure 15a shows the original contour polylines. First, search the original intersections of the contour polylines, and the result is shown in Figure 15b. The green points are the original intersections that need no further processing. It can be seen that only a few polylines possess the original intersections. Red points are isolated points that have no other contour polylines within the specified scope. Then, the contour polylines are moved to intersect within the tolerance range. Figure 15c shows the intersections after processing, and most of the contour polylines have intersected with each other. After that, the closed-loop division is carried out based on the intersected contour polylines, and the result is shown in Figure 15d. Based on the closed loops, sub-meshes can be constructed through the Coons surface interpolation. Figure 15e shows the modeling result of every closed loop. The sub-meshes are smooth but split from each other. Finally, the sub-meshes are combined to obtain a complete orebody model, as shown in Figure 15f.

As shown in Figure 16, the method is applied to complex orebodies with different shapes, and good orebody models are obtained, showing that this method has strong adaptability.

4.2. Discussion

The above experimental results show that the proposed method can deal with the contour polylines obtained by geological logging in mine production with good effects. Firstly, the constructed model conforms to the geological trend and is smooth in every mesh. Secondly, the closed-loop division can be carried out automatically, which will greatly improve the efficiency of modeling. Thirdly, since modeling between different regions will not affect each other, it can be carried out at the same time, which will greatly reduce the overall modeling time. However, this method still has some limitations that need to be further studied and solved. Firstly, due to the inherent limitations of the Coons technology, the model constructed by Coons surface interpolation is difficult to update, and it is difficult to take faulting crosscutting into account, which will affect the accuracy of the model. In mine production, many orebodies are cross-cut by faults. For these orebodies, we should choose other modeling methods. Secondly, due to the complexity of the data in actual production, when this method is applied to some complex contour polylines, as shown in Figure 15, it can not complete all closed-loop divisions well, which needs to be specified manually. Finally, before Coons surface interpolation, all regions except simple single-sided regions need to be transformed into four-sided regions, which may affect the final modeling results by introducing numerical instability issues when adding virtual edges and vertexes. Therefore, to further improve the efficiency and accuracy of the modeling of cross-contour polylines, it is necessary to develop a closed-loop division method with higher precision and an n-sided regions expansion method with less effect on interpolation modeling in the future.

5. Conclusions

We propose a modeling method based on bicubic Coons surface interpolation, which needs the contour polyline to be in an approximate plane within the tolerance range. The innovation of this paper is that we proposed a creative method of automatically dividing the closed loops, which will greatly reduce the modeling time. Additionally, the constructed sub-meshes can be combined without manual intervention. This method solves the problem that it is difficult to divide closed loops based on contour polylines with complex shapes, and it greatly improves the efficiency of modeling based on complex cross-contour polylines. Firstly, the contour polylines are cut by the approximate planes to divide the subspaces. Then, the divided polylines are grouped and connected according to the specified tolerance to form closed loops. Secondly, the Coons surface interpolation is used for modeling based on the closed-loop information. Finally, the constructed sub-meshes are combined to form a complete orebody model. Through experiments with multiple groups of cross-contour polylines, the results show that the proposed method can accurately divide closed loops and model the contour polylines with good effects.