A Universal Generating Algorithm of the Polyhedral Discrete Grid Based on Unit Duplication

Abstract

Based on the analysis of the problems in the generation algorithm of discrete grid systems domestically and abroad, a new universal algorithm for the unit duplication of a polyhedral discrete grid is proposed, and its core is “simple unit replication + effective region restriction”. First, the grid coordinate system and the corresponding spatial rectangular coordinate system are established to determine the rectangular coordinates of any grid cell node. Then, the type of the subdivision grid system to be calculated is determined to identify the three key factors affecting the grid types, which are the position of the starting point, the length of the starting edge, and the direction of the starting edge. On this basis, the effective boundary of a multiscale grid can be determined and the grid coordinates of a multiscale grid can be obtained. A one-to-one correspondence between the multiscale grids and subdivision types can be established. Through the appropriate rotation, translation and scaling of the multiscale grid, the node coordinates of a single triangular grid system are calculated, and the relationships between the nodes of different levels are established. Finally, this paper takes a hexagonal grid as an example to carry out the experiment verifications by converting a single triangular grid system (plane) directly to a single triangular grid with a positive icosahedral surface to generate a positive icosahedral surface grid. The experimental results show that the algorithm has good universality and can generate the multiscale grid of an arbitrary grid configuration by adjusting the corresponding starting transformation parameters.

1. Introduction

The Discrete Global Grid System (DGGS) is a new kind of organization, management, and application of the spatial data model, which uses specific methods to evenly discretize the Earth into a seamless, non-overlapping multi-resolution grid hierarchy [1,2]. The cells of DGGS may be triangular, quadrilateral (squares or diamonds), or hexagonal. Dutton [3,4], Goodchild [5], Szalay [6], and so on, studied the triangulation based on the regular octahedron, as shown in Figure 1a. Alborzi and Samet et al. [7] studied the cube-based quadrangle subdivision, as shown in Figure 1b. White et al. [8] studied the diamond subdivision based on the regular icosahedron, as shown in Figure 1c. Sahr [9], Peterson [10], Vince [11], Zheng [12], and Tong et al. [13] studied the global grid of hexagon subdivision on the regular icosahedron, as shown in Figure 1d. DGGS has been widely applied in many fields of Geoscience, such as Earth system modeling [14], global spatial data indexing [15,16,17,18], global environmental monitoring [8,19], global meteorological simulation [20], global navigating [21], ocean path planning [22], optimal path determination [23], cellular automata [24], and so on.

In recent years, different types of normal polyhedral grid generation methods have been proposed at home and abroad, such as generation methods for Quaternary Triangular Mesh (QTM) [25], the hierarchical grid generation algorithm for diamond grids [26], and the regular polyhedron generation algorithm for hexagonal grids [9]. Among these research results, there are mainly two ideas regarding the construction of grid units of different levels in the generation of a polyhedral grid: One is to subdivide using recursion layer by layer and the other is to subdivide using the hierarchical arrangement unit by unit. While these methods can generate different kinds of triangular, quadrilateral, and hexagon grids, there are two kinds of problems in the traditional multi-scale grids generation methods:

• It is necessary to customize the corresponding grid generation algorithm according to the shape (triangle, quadrilateral, or hexagon) of the target grid, the aperture (3 or 4 aperture) [9], the prime factor, the duality, the centricity of the subdivision frequency [1], and other factors [8]. The recursive subdivision methods of different partition topologies are shown in Figure 2.
• In the process of the recursion grid generation method, an effective distributed algorithm with high algorithmic complexity needs to be designed for the generation of high-level grids.

In view of the above questions, this paper proposes a universal global discrete grid generation algorithm based on the unit replication method by analyzing different grid cell types and different grid levels. This method adopts the “simple unit replication + effective region restriction” method, which can solve the universal generation problem of different grid forms and different grid levels. Compared with the recursive splitting method, its significant advantages are as follows: (1) In the process of generating ultra-high-level grid cells, the complexity of the algorithm is reduced, and the grid generation efficiency is greatly improved; and (2) it highlights the potential links between different types of grid cells, and promotes the theoretical research on communication, sharing and interoperability of global discrete grid systems. This paper elaborates the generation algorithm based on the hexagonal grid of a regular icosahedron, which can be adopted by generating other grid types. The generation method of the global discrete grid adopts the method shown in Figure 1 in Reference [27]. This paper mainly focuses on the construction algorithm of multi-scale grid systems.

2. Algorithm Thoughts

Since the regular icosahedron is centrally symmetric, we only consider the grid subdivision of a triangular surface, which is defined as the single triangular grid system. After determining the basic configuration of the grid, we analyzed the existing single-triangular grid systems with various subtypes and found that the essential differences lie in the number of subgrids that are contained in the single-triangular grid system and their positional attitude relative to the triangular grid. Therefore, from another perspective, we fix the basic subdivision form of the grid and extend it to an infinite plane to make it seamless and overlapping, which we call the “base extended planar grid system.” The basic idea of the algorithm will be introduced as follows.

In the base extended planar grid system, an appropriate discrete grid skew coordinate system is established. First, fix the coordinates of the vertices of an initial equilateral triangle, which is called the initial single triangle. By appropriately rotating, translating, and scaling the coordinates of the boundary vertices of the initial single triangle, the initial simple triangle can be grouped with some grid cells of the base extended planar grid. In this way, different combinations of transformation parameters are used to establish a one-to-one corresponding mathematical relationship with a single triangular grid system of specific subdivision types and levels. By changing the configurations of the grid and calculating the transformation parameters, a single triangular grid system with arbitrary cross-sections of triangular, quadrilateral, and hexagonal arbitrary segmentation types can be generated.

In summary, the algorithm is divided into the following steps:

(1)

establish the skew coordinate system of the base expanded grid,

(2)

determine the type of the grid system to be calculated,

(3)

determine the effective control boundaries of the grids of different levels,

(4)

calculate the node coordinates of the grid cells of different levels,

(5)

generate the single triangular grid system coordinates,

(6)

establish the indexing relationship between the grids of different levels.

The specific algorithm flowchart is shown in Figure 3:

3. Methodology

3.1. Establishment of the Skew Coordinate System of the Base Expanded Grid

The base extended planar grid coordinate system is a skew coordinate system that is established on the infinite plane discrete grid, and its components are called grid coordinates in the following context. It is mainly used to determine the central position of the grid cells in a single triangular grid system and to describe the regular arrangement of grid types. The coordinate values are all integers. For the convenience of this research, this paper takes the hexagonal grid as an example (the grid that is mentioned in the following refers to a hexagonal grid). The coordinate origin is defined in the center of any grid cell in the infinite plane discrete grid. The horizontal adjacent grid centerline is used as the coordinate system’s I axis, and the coordinate system’s J axis is rotated 120° counterclockwise in the direction of the I axis.

In addition, the auxiliary spatial Cartesian coordinates should also be established based on the origin of the coordinate system. The origin and the X axis of the coordinate system coincide with the origin of the base extended grid oblique coordinate system and the I axis. The Y-axis is rotated 90° counterclockwise on the X-axis, and the Z-axis follows the right-handed spiral rule, the components are hereinafter referred to as Cartesian coordinates.

In the extended grid skew coordinate system, the radius of a grid cell is defined as R, and the vertex grid coordinates of the initial single triangle are $A_0\left(0,0\right)$, $B_0\left(1,0\right)$, and $C_0\left(1,1\right)$, respectively.

Of the seven nodes of a grid cell with the grid coordinates (i,j), the central node is identified as $O_{\left(i,j\right)}$. The six boundary nodes are identified in the counterclockwise direction from the upper right corner as $A_{\left(i,j\right)}$, $B_{\left(i,j\right)}$, $C_{\left(i,j\right)}$, $D_{\left(i,j\right)}$, $E_{\left(i,j\right)}$, and $F_{\left(i,j\right)}$. The specific situation is shown in Figure 4.

Through the geometrical transformation between the extended grid skew coordinate system and the spatial Cartesian coordinate system, the spatial rectangular coordinates of the seven nodes of any grid unit $\left(i,j\right)$ in the extended grid skew coordinate system are obtained, as in Formula (1).

$$\{\begin{array}{c}{A}_{\left(i,j\right)}\left(\frac{2}{\sqrt{3}}R+i·\sqrt{3}R-\frac{1}{2}\times (\frac{1}{\sqrt{3}}R+j·\sqrt{3}R),\frac{\sqrt{3}}{2}\times (\frac{1}{\sqrt{3}}R+j·\sqrt{3}R),0\right)\\ B_{\left(i,j\right)}\left(\frac{1}{\sqrt{3}}R+i·\sqrt{3}R-\frac{1}{2}\times (\frac{2}{\sqrt{3}}R+j·\sqrt{3}R),\frac{\sqrt{3}}{2}\times (\frac{2}{\sqrt{3}}R+j·\sqrt{3}R),0\right)\\ C_{\left(i,j\right)}\left(\frac{-1}{\sqrt{3}}R+i·\sqrt{3}R-\frac{1}{2}\times (\frac{1}{\sqrt{3}}R+j·\sqrt{3}R),\frac{\sqrt{3}}{2}\times (\frac{1}{\sqrt{3}}R+j·\sqrt{3}R),0\right)\\ D_{\left(i,j\right)}\left(\frac{-2}{\sqrt{3}}R+i·\sqrt{3}R-\frac{1}{2}\times (\frac{-1}{\sqrt{3}}R+j·\sqrt{3}R),\frac{\sqrt{3}}{2}\times (\frac{-1}{\sqrt{3}}R+j·\sqrt{3}R),0\right)\\ E_{\left(i,j\right)}\left(\frac{-1}{\sqrt{3}}R+i·\sqrt{3}R-\frac{1}{2}\times (\frac{-2}{\sqrt{3}}R+j·\sqrt{3}R),\frac{\sqrt{3}}{2}\times (\frac{-2}{\sqrt{3}}R+j·\sqrt{3}R),0\right)\\ F_{\left(i,j\right)}\left(\frac{1}{\sqrt{3}}R+i·\sqrt{3}R-\frac{1}{2}\times (\frac{-1}{\sqrt{3}}R+j·\sqrt{3}R),\frac{\sqrt{3}}{2}\times (\frac{-1}{\sqrt{3}}R+j·\sqrt{3}R),0\right)\\ O_{\left(i,j\right)}\left(0+i·\sqrt{3}R-\frac{1}{2}\times (0+j·\sqrt{3}R),\frac{\sqrt{3}}{2}\times (0+j·\sqrt{3}R),0\right)\end{array}$$

(1)

3.2. Determination the Type of the Grid System to Be Calculated

This paper takes the hexagonal grid system as an example. The hexagonal grid is more complicated. At present, the aperture 3 subdivision and aperture 4 subdivision are widely used in the hexagonal subdivisions. There is one common type of aperture 3 subdivision and four common types of aperture 4 subdivision, which are shown in Figure 5.

Based on the analysis of the morphology of the various single triangular grid systems of the common hexagonal grid system, it is found that the geometries of the grids are the same. The only difference is that the numbers of grid cells that are contained in the triangulation planes are different from the positions and attitudes of these grids relative to the triangulation surface, and the specific situations are similar for different levels of grids and different types of grids. Therefore, the morphology of a single triangular grid in the extended grid skew coordinate system is mainly determined by three factors: The position of the starting point S = (i₀, j₀), the length of the starting edge L_N, and the angle θ between the starting edge and the I axial direction.

There are only two cases for the starting point location, which are starting from the grid center node and starting from a grid boundary node. Different values of the length of the starting edge can generate grid systems of different scales. θ is the angle of the starting edge and the positive direction of the X axis, which has three conditions: 0°, 30°, and 60°. The detailed subdivision of a hexagon grid in the extended grid skew coordinate system is shown in Figure 6.

Through the comparative analysis of the subdivision types of Figure 5 and Figure 6, and based on the grid morphology at the triangle boundary (the specific situation has been marked in the gray area of Figure 5), a conclusion can be obtained that the subdivision method in Figure 6a corresponds to the common subdivision of the class I aperture 4 grid and class II grid and the even level of the aperture 3 grid, Figure 6b corresponds to the class III aperture 4 grid and the odd level of the aperture 3 grid, and Figure 6c corresponds to the class IV aperture 4 grid. The classification method in Figure 6d does not have an obvious geometric symmetry relationship between the grid form at the triangular boundary and the grid cells, which leads to the extremely complicated ownership problem of the cross-surface elements between adjacent triangular surfaces [25]. Therefore, the paper does not extensively address this kind of subdivision method.

3.3. Determination of the Effective Control Boundary of the Grids of Different Levels

To illustrate the problem, this paper takes the subdivision method in Figure 6a as an example. It is found that the key factor of the multiscale grid generation depends on the length $L_N$ of the starting edge, by analyzing the single triangular facet grid system of this type of subdivision in the extended grid skew coordinate system. In other words, different starting edge lengths can be used to produce grid systems with different levels under this subdivision method. Therefore, it is necessary to study the mathematical relationship between different starting edge lengths and grid levels.

The 20 vertices of an icosahedron have the characteristics of symmetry, equality and reflexivity. Therefore, with the fixed subdivision method, the grid morphology at the vertex of any triangular boundary must be consistent. With this consideration, we find that the value of the starting edge length $L_N$ can be expressed as $L_N=M·\sqrt{3}R$. In the formula, $M\in Z^{+}$, where $R$ is the radius of the hexagonal grid unit.

To reflect the relationship between the different levels of the grid system, through the study of the hierarchical grid, the relationship between the specific value $L_N$ of the starting edge and the grid level is expressed as follows: when $L_N=3·2^{N-1}·\sqrt{3}R$, take the appropriate value of N corresponding to the length of the grid boundary of level N.

Therefore, the grid unit where the starting node is located is defined as $S=i_0,j_0$, the initial transformation parameter is changed to $F=(S,L_N,\theta )=\left(\right(i_0,j_0),3·2^{N-1}·\sqrt{3}R,{60}^{°})$, and the grid coordinates of the three vertices of the triangle of layer N are $A_N(i_0,j_0)$, $B_N(i_0+3·2^{N-1},j_0)$, and $C_N(i_0+3·2^{N-1},j_0+3·2^{N-1})$. Then, the coordinates of the boundary vertex of the multiscale grid in the spatial Cartesian coordinate system can be rapidly determined by Formula (1). The specific result is shown in Formula (2).

$$\{\begin{array}{c}{A}_N\left(\left(i_0-\frac{1}{2}{j}_0\right)·\sqrt{3}R,\frac{\sqrt{3}}{2}{j}_0·\sqrt{3}R,0\right) \\ B_N\left(\left(i_0-\frac{1}{2}{j}_0+3·2^{N-1}\right)·\sqrt{3}R,\frac{\sqrt{3}}{2}{j}_0·\sqrt{3}R,0\right) \\ C_N\left(\left(i_0-\frac{1}{2}{j}_0+3·2^{N-2}\right)·\sqrt{3}R,\frac{\sqrt{3}}{2}\left(j_0+3·2^{N-1}\right)·\sqrt{3}R,0\right)\end{array}$$

(2)

3.4. Calculation of the Node Coordinates of Grid Cells of Different Levels

According to the above research, given the grid coordinates $S=(i_0,j_0)$ and the grid division level N of any starting node, in the base extended grid skew coordinate system, the grid coordinates $\left(i,j\right)$ of any grid cells in the grid division can be rapidly calculated. The specific values of the grid coordinates are calculated by Formula (3). The details of the multiscale grid boundary in the Nth layer in the extended grid skew coordinate system are shown in Figure 7. Therefore, the spatial rectangular coordinates of the seven nodes of any grid cells in the Nth layer can be obtained by Formula (1).

$$\{\begin{array}{c}{i}_0\le i\le i_0+3·2^{N-1}\\ j_0\le j\le j_0+3·2^{N-1}\\ j\le i-\left(i_0-j_0\right) \end{array}$$

(3)

3.5. Generation of the Single Triangular Grid System Coordinates

In the Nth layer, when the grid coordinates $\left(i,j\right)$ satisfy Formula (3), the components of the center node of the seven nodes that are mentioned in Section 3.1 for the spatial Cartesian coordinate system are expressed as $O_{\left(i,j\right)}^X$, $O_{\left(i,j\right)}^Y$ and $O_{\left(i,j\right)}^Z$ for the X-axis, Y-axis and Z-axis, respectively. The components of the six boundary nodes, which are defined counterclockwise from the upper right, on the X-axis, Y-axis and Z-axis in the space Cartesian coordinate system are expressed as follows: $A_{\left(i,j\right)}^X$, $A_{\left(i,j\right)}^Y$ and $A_{\left(i,j\right)}^Z$; $B_{\left(i,j\right)}^X$, $B_{\left(i,j\right)}^Y$ and $B_{\left(i,j\right)}^Z$; $C_{\left(i,j\right)}^X$, $C_{\left(i,j\right)}^Y$ and $C_{\left(i,j\right)}^Z$; $D_{\left(i,j\right)}^X$, $D_{\left(i,j\right)}^Y$ and $D_{\left(i,j\right)}^Z;$ $E_{\left(i,j\right)}^X$ $E_{\left(i,j\right)}^Y$ and $E_{\left(i,j\right)}^Z$; $F_{\left(i,j\right)}^X$, $F_{\left(i,j\right)}^Y$, and $F_{\left(i,j\right)}^Z$.

Through the geometric relationship between the multiscale grid and the boundary vertices of the initial single triangle in the spatial Cartesian coordinate system, the following coordinate transformation relationship is given:

$$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=K·H·\left[\begin{array}{c}x\\ y\\ z\end{array}\right]+\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]$$

(4)

In Formula (4), $K$ represents the scaling coefficient, $H$ represents the rotation matrix, ${\left[\begin{array}{ccc}{x}_0& y_0& z_0\end{array}\right]}^T$ represents the translation vector, and ${\left[\begin{array}{ccc}x& y& z\end{array}\right]}^T$ and ${\left[\begin{array}{ccc}X& Y& Z\end{array}\right]}^T$ represent the coordinate values before and after the transformation, respectively. There are two specific situations.

• When $Z=0$, there is a only single triangular grid system unit (plane) and H is a two-dimensional rotation matrix.
• When $Z\ne 0$, a transformation from the initial single triangular face (plane) to the spatial triangular face is required, through which the grid of the icosahedral surface can be directly generated.

In view of the abovementioned subdivision method, this study shows that given the initial transformation parameter $F=(S,L_N,\theta )=\left(\right(i_0,j_0),3·2^{N-1}·\sqrt{3}R,{60}^{°})$, through the similar scaling transformation, where $K=\frac{1}{3·2^{N-1}}$, $H=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$, and $\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]=-\frac{1}{3·2^N}·\left[\begin{array}{c}{O}_{\left(i_0,j_0\right)}^X\\ O_{\left(i_0,j_0\right)}^Y\\ O_{\left(i_0,j_0\right)}^Z\end{array}\right]$, when the coordinate of the grid (i,j) in the Nth layer of the extended grid skew coordinate system satisfies Formula (3), the rectangular coordinates of the seven nodes of the corresponding grid unit in the single triangular grid system can be calculated as follows:

$$\{\begin{array}{c}{A}_{\left(i,j\right)}^{\prime }\left(\frac{1}{3·2^{N-1}}\left(A_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{3·2^{N-1}}\left(A_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{3·2^{N-1}}\left(A_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ B_{\left(i,j\right)}^{\prime }\left(\frac{1}{3·2^{N-1}}\left(B_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{3·2^{N-1}}\left(B_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{3·2^{N-1}}\left(B_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ C_{\left(i,j\right)}^{\prime }\left(\frac{1}{3·2^{N-1}}\left(C_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{3·2^{N-1}}\left(C_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{3·2^{N-1}}\left(C_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ D_{\left(i,j\right)}^{\prime }\left(\frac{1}{3·2^{N-1}}\left(D_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{3·2^{N-1}}\left(D_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{3·2^{N-1}}\left(D_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ E_{\left(i,j\right)}^{\prime }\left(\frac{1}{3·2^{N-1}}\left(E_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{3·2^{N-1}}\left(E_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{3·2^{N-1}}\left(E_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ F_{\left(i,j\right)}^{\prime }\left(\frac{1}{3·2^{N-1}}\left(F_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{3·2^{N-1}}\left(F_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{3·2^{N-1}}\left(F_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ O_{\left(i,j\right)}^{\prime }\left(\frac{1}{3·2^{N-1}}\left(O_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{3·2^{N-1}}\left(O_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{3·2^{N-1}}\left(O_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\end{array}$$

(5)

3.6. Establishment of the Relationship between the Grids of Different Levels

The establishment of the relationship between the grids of different levels is conducive to the rapid retrieval of the spatial data and the efficient calculation of the application analysis. The indexing relationship here refers to the correspondence between the grid nodes of different levels. The goal is to establish the correlation of the nodes of different levels, that is, to calculate the specific grid coordinates and node locations of the seven nodes in the (N + 1)th layer grid in the grid unit of the Nth grid coordinates $\left(I,J\right)$.

In the extended grid skew coordinate system shown in Figure 4, the positions of the seven nodes of the grid $\left(I,J\right)$ in the Nth layer of the single triangular grid system are expressed as $A_{\left(N,I,J\right)}$, $B_{\left(N,I,J\right)}$, $C_{\left(N,I,J\right)}$, $D_{\left(N,I,J\right)}$, $E_{\left(N,I,J\right)}$, $F_{\left(N,I,J\right)}$, and $O_{\left(N,I,J\right)}$. The arrangement of the positions between the seven nodes follows the aforementioned principles. By analyzing the geometric positioning relationship of the subunits and parent units of the adjacent hierarchical grids, the following conclusions are drawn:

$$\{\begin{array}{c}{A}_{\left(N,I,J\right)}~B_{\left(N+1,2I+1,2J\right)} \\ B_{\left(N,I,J\right)}~C_{\left(N+1,2I+1,2J+1\right)} \\ C_{\left(N,I,J\right)}~D_{\left(N+1,2I,2J+1\right)} \\ D_{\left(N,I,J\right)}~E_{\left(N+1,2I-1,2J\right)} \\ E_{\left(N,I,J\right)}~F_{\left(N+1,2I-1,2J-1\right)} \\ F_{\left(N,I,J\right)}~A_{\left(N+1,2I,2J-1\right)} \\ O_{\left(N,I,J\right)}~O_{\left(N+1,2I,2J\right)} \end{array}$$

(6)

In Formula (6), $A_{\left(N,I,J\right)}~B_{\left(N+1,2I+1,2J\right)}$ indicates that boundary node A of the grid $\left(I,J\right)$ in the Nth layer corresponds to boundary node B of the grid $\left(2I+1,2J\right)$ in the (N + 1)-th layer. The other cases in Formula (6) are similar.

3.7. Some Key Parameters

According to the previous algorithm, the key parameters in Section 3.1, Section 3.2, Section 3.3, Section 3.4, Section 3.5 and Section 3.6 in the outlined single triangular grid system are calculated, as shown in

Table 1. The initial parameters of the various split grid types, the effective control of the boundary, and the unit node spatial Cartesian coordinates

	Grid parameter	Grid Type	Starting Position	$\mathbf{Relationship}\mathbf{between}\mathbf{Starting}\mathbf{Edge}\mathbf{Length}{\mathit{L}}_{\mathit{N}}$$\mathbf{and}\mathbf{Level}\mathit{N}$	$\varnothing$	Grid Node Coordinates of the Nth Level	Grid Coordinate Range of the Nth Level	Rectangular Coordinates of Unit Nodes in the Single Triangular Grid System
Grid Shape at Triangular Boundaries
		class I aperture 4 grid	$O_{\left(i_0,j_0\right)}$	$L_N=3·2^{N-1}·\sqrt{3}R$	$60°$	$\{\begin{array}{c}{A}_N\left(i_0,j_0\right) \\ B_N\left(i_0+3·2^{N-1},j_0\right) \\ C_N\left(i_0+3·2^{N-1},j_0+3·2^{N-1}\right)\end{array}$	$\{\begin{array}{c}{i}_0\le i\le i_0+3·2^{N-1}\\ j_0\le j\le j_0+3·2^{N-1}\\ j\le i-\left(i_0-j_0\right)\end{array}$	Formula (5)
		class II aperture 4 grid	$O_{\left(i_0,j_0\right)}$	$L_N=2^N·\sqrt{3}R$	$60°$	$\{\begin{array}{c}{A}_N\left(i_0,j_0\right) \\ B_N\left(i_0+2^N,j_0\right) \\ C_N\left(i_0+2^N,j_0+2^N\right)\end{array}$	$\{\begin{array}{c}{i}_0\le i\le i_0+2^N\\ j_0\le j\le j_0+2^N\\ j\le i-\left(i_0-j_0\right)\end{array}$	Formula (7)
		aperture 3 even grid	$O_{\left(i_0,j_0\right)}$	$L_N={\left(\sqrt{3}\right)}^N·\sqrt{3}R$	$60°$	$\{\begin{array}{c}{A}_N\left(i_0,j_0\right) \\ B_N\left(i_0+{\left(\sqrt{3}\right)}^N,j_0\right) \\ C_N\left(i_0+{\left(\sqrt{3}\right)}^N,j_0+{\left(\sqrt{3}\right)}^N\right)\end{array}$	$\{\begin{array}{c}{i}_0\le i\le i_0+{\left(\sqrt{3}\right)}^N\\ j_0\le j\le j_0+{\left(\sqrt{3}\right)}^N\\ j\le i-\left(i_0-j_0\right) \end{array}$	Formula (8)
		aperture 3 odd grid	$O_{\left(i_0,j_0\right)}$	$L_N={\left(\sqrt{3}\right)}^{N+1}·R$	$30°$	$\{\begin{array}{c}{A}_N\left(i_0,j_0\right) \\ B_N\left(i_0+{\left(\sqrt{3}\right)}^{N-1},j_0-{\left(\sqrt{3}\right)}^{N-1}\right) \\ C_N\left(i_0+2·{\left(\sqrt{3}\right)}^{N-1},j_0+{\left(\sqrt{3}\right)}^{N-1}\right)\end{array}$	$\{\begin{array}{c}i+j\ge i_0+j_0 \\ i-2j\ge i_0-2j_0 \\ 2i-j\le 2i_0-j_0+{\left(\sqrt{3}\right)}^{N+1}\end{array}$	Formula (9)
		class III aperture 4 grid	$O_{\left(i_0,j_0\right)}$	$L_N=6·2^{N-1}·R$	$30°$	$\{\begin{array}{c}{A}_N\left(i_0,j_0\right) \\ B_N\left(i_0+2^N,j_0-2^N\right) \\ C_N\left(i_0+2^{N+1},j_0+2^N\right)\end{array}$	$\{\begin{array}{c}i+j\ge i_0+j_0 \\ i-2j\ge i_0-2j_0 \\ 2i-j\le 2i_0-j_0+6·2^{N-1}\end{array}$	Formula (10)
		class IV aperture 4 grid	$B_{\left(i_0,j_0\right)}$	$L_N=3·2^{N-1}·R$	$30°$	$\{\begin{array}{c}{A}_N\left(i_0,j_0\right) \\ B_N\left(i_0+2^{N-1},j_0-2^{N-1}\right)\\ C_N\left(i_0+2^N,j_0+2^{N-1}\right) \end{array}$	$\{\begin{array}{c}i+j\ge i_0+j_0 \\ i-2j\ge i_0-2j_0 \\ 2i-j\le 2i_0-j_0+3·2^{N-1}\end{array}$	Formula (11)

Note: Formula (5) in the table is shown in Section 3.5. Formulas (7)–(11) in the table are shown below.

.

$$\{\begin{array}{c}{A}_{\left(i,j\right)}^{\prime }\left(\frac{1}{2^N}\left(A_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{2^N}\left(A_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{2^N}\left(A_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ B_{\left(i,j\right)}^{\prime }\left(\frac{1}{2^N}\left(B_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{2^N}\left(B_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{2^N}\left(B_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ C_{\left(i,j\right)}^{\prime }\left(\frac{1}{2^N}\left(C_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{2^N}\left(C_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{2^N}\left(C_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ D_{\left(i,j\right)}^{\prime }\left(\frac{1}{2^N}\left(D_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{2^N}\left(D_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{2^N}\left(D_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ E_{\left(i,j\right)}^{\prime }\left(\frac{1}{2^N}\left(E_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{2^N}\left(E_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{2^N}\left(E_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ F_{\left(i,j\right)}^{\prime }\left(\frac{1}{2^N}\left(F_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{2^N}\left(F_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{2^N}\left(F_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ O_{\left(i,j\right)}^{\prime }\left(\frac{1}{2^N}\left(O_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{2^N}\left(O_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{2^N}\left(O_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\end{array}$$

(7)

$$\{\begin{array}{c}{A}_{\left(i,j\right)}^{\prime }\left(\frac{1}{{\left(\sqrt{3}\right)}^N}\left(A_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{{\left(\sqrt{3}\right)}^N}\left(A_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{{\left(\sqrt{3}\right)}^N}\left(A_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ B_{\left(i,j\right)}^{\prime }\left(\frac{1}{{\left(\sqrt{3}\right)}^N}\left(B_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{{\left(\sqrt{3}\right)}^N}\left(B_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{{\left(\sqrt{3}\right)}^N}\left(B_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ C_{\left(i,j\right)}^{\prime }\left(\frac{1}{{\left(\sqrt{3}\right)}^N}\left(C_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{{\left(\sqrt{3}\right)}^N}\left(C_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{{\left(\sqrt{3}\right)}^N}\left(C_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ D_{\left(i,j\right)}^{\prime }\left(\frac{1}{{\left(\sqrt{3}\right)}^N}\left(D_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{{\left(\sqrt{3}\right)}^N}\left(D_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{{\left(\sqrt{3}\right)}^N}\left(D_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ E_{\left(i,j\right)}^{\prime }\left(\frac{1}{{\left(\sqrt{3}\right)}^N}\left(E_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{{\left(\sqrt{3}\right)}^N}\left(E_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{{\left(\sqrt{3}\right)}^N}\left(E_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ F_{\left(i,j\right)}^{\prime }\left(\frac{1}{{\left(\sqrt{3}\right)}^N}\left(F_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{{\left(\sqrt{3}\right)}^N}\left(F_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{{\left(\sqrt{3}\right)}^N}\left(F_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\\ O_{\left(i,j\right)}^{\prime }\left(\frac{1}{{\left(\sqrt{3}\right)}^N}\left(O_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\right),\frac{1}{{\left(\sqrt{3}\right)}^N}\left(O_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\right),\frac{1}{{\left(\sqrt{3}\right)}^N}\left(O_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\right)\right)\end{array}$$

(8)

$$\left[\begin{array}{c}{O}_{\left(i,j\right)}^{X^{\prime }}\\ O_{\left(i,j\right)}^{Y^{\prime }}\\ O_{\left(i,j\right)}^{Z^{\prime }}\end{array}\right]=\frac{1}{{\left(\sqrt{3}\right)}^N}·\left[\begin{array}{ccc}\frac{\sqrt{3}}{2}& \frac{-1}{2}& 0\\ \frac{1}{2}& \frac{\sqrt{3}}{2}& 0\\ 0& 0& 1\end{array}\right]·\left[\begin{array}{c}{O}_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\\ O_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\\ O_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\end{array}\right]$$

(9)

$$\left[\begin{array}{c}{O}_{\left(i,j\right)}^{X^{\prime }}\\ O_{\left(i,j\right)}^{Y^{\prime }}\\ O_{\left(i,j\right)}^{Z^{\prime }}\end{array}\right]=\frac{1}{\sqrt{3}·2^N}·\left[\begin{array}{ccc}\frac{\sqrt{3}}{2}& \frac{-1}{2}& 0\\ \frac{1}{2}& \frac{\sqrt{3}}{2}& 0\\ 0& 0& 1\end{array}\right]·\left[\begin{array}{c}{O}_{\left(i,j\right)}^X-O_{\left(i_0,j_0\right)}^X\\ O_{\left(i,j\right)}^Y-O_{\left(i_0,j_0\right)}^Y\\ O_{\left(i,j\right)}^Z-O_{\left(i_0,j_0\right)}^Z\end{array}\right]$$

(10)

$$\left[\begin{array}{c}{O}_{\left(i,j\right)}^{X^{\prime }}\\ O_{\left(i,j\right)}^{Y^{\prime }}\\ O_{\left(i,j\right)}^{Z^{\prime }}\end{array}\right]=\frac{1}{\sqrt{3}·2^{N-1}}·\left[\begin{array}{ccc}\frac{\sqrt{3}}{2}& \frac{-1}{2}& 0\\ \frac{1}{2}& \frac{\sqrt{3}}{2}& 0\\ 0& 0& 1\end{array}\right]·\left[\begin{array}{c}{O}_{\left(i,j\right)}^X-B_{\left(i_0,j_0\right)}^X\\ O_{\left(i,j\right)}^Y-B_{\left(i_0,j_0\right)}^Y\\ O_{\left(i,j\right)}^Z-B_{\left(i_0,j_0\right)}^Z\end{array}\right]$$

(11)

In Formulas (9)–(11), the rectangular coordinates of the other six boundary nodes are similar to those of the center node O. By replacing the node letter O of the grid (i,j) unit with the letters A, B, C, D, E, and F, the respective coordinates can be obtained.

In addition, this paper also gives a direct conversion to the regular polyhedral surface grid through this method. It only needs to make Z > 0 in Equation (4), and then combine the geometric property information of the regular polyhedron in order to establish a suitable mathematical transformation relationship. The transformation from the initial single triangular face (plane) to the spatial triangular face can be done to directly generate a grid of regular polyhedral surfaces from the base extended planar grid system.

4. Experiments and Discussion

4.1. Experiments

This paper mainly carried out the following two experiments:

Experiment 1. 

According to the above algorithm, this paper takes the class I aperture 4 grid as an example. By querying the relevant parameters of the class I aperture 4 grid in Table 1 and combining Formula (4) and the spatial geometry relationship of the regular icosahedron, the multiscale grid of this kind on the regular icosahedron can be directly generated.

Experimental environment: ASUS All Series desktop computer, Win8.1 Enterprise 64 bit operating system, Intel fourth generation Core i5-4460 3.20 GHz quad core processor, 8 GB RAM, 7200 rpm 500 GB hard disk drive, and AMD Radeon R7 200 Series graphics card. Compilation environment: Visual Studio 2013, Release version, C++, and OSG.

The specific experimental results are shown in Figure 8.

Experiment 2. 

According to the above algorithm, the comparative experiments are carried out in this paper with four representative subdivision types, which are the class I aperture 4, the class II aperture 4, the class III aperture 4 and the aperture 3 subdivision grid. Referring to the grid subdivision data in Table 1, the Earth is approximately represented by a regular icosahedron, and the radius of the Earth is the average radius of the WGS-84 reference ellipsoid [1], i.e., R = 6371008771 m. Taking Snyder’s iso-product polyhedron projection as the basic projection type, the mapping relationship between the plane and the sphere is established, and various spherical hexagonal multiscale grid systems are generated. The generation process of the global discrete grid system is shown in Figure 1 of Reference [28].

The relevant experimental environment and compilation environment are the same as in Experiment 1. The results of the specific experimental procedures are shown in Figure 9.

4.2. Discussion

On the basis of previous studies, this paper makes an in-depth analysis of the production methods of multiscale grid elements. It is found that the geometric structure of a single triangular mesh system with different hexagonal triangulation types is different, but its common feature is that the shape of the mesh is the same. The only difference is that the numbers of meshes that are contained in the triangles are different with respect to the positions and attitudes of these meshes relative to the triangles, and the different levels of meshes are different. The grids and the different mesh types are similar. Based on the above characteristics, a unified base extended grid oblique coordinate system is established, and a multiscale single triangular grid system of an arbitrary hexagonal section is generated by changing the initial parameters. In addition, the method is equally applicable to the twelve pentagons at the vertex of an icosahedron. Since the method of this paper is to use a region to restrict the range of grids, and then to make spatial similarity transformation and stitching according to the coordinates of icosahedron vertexes, so each vertex is composed of five triangles because of the characteristics of the icosahedron itself. Therefore, 1/6 of a hexagonal grid at the vertices constrained by triangles (as shown in Table 1) automatically form a pentagonal grid, and 12 vertices form 12 pentagonal grids.

The advantages of the algorithm are as follows:

(1)

Compared with the existing hierarchical recursive subdivision, it is not necessary to customize the exclusive generation algorithm for a specific aperture, subdivision frequency, and eccentricity, which provides convenience for different split type grid generation algorithms. Furthermore, it highlights the potential links between the various types of grids and strengthens the theoretical research on the communication, sharing and interoperability among the global discrete grid systems.

(2)

Compared with hierarchical recursive subdivision, in the process of high-level grid generation, it is not necessary to design an effective distributed algorithm, which reduces the complexity of the algorithm.

(3)

In fact, the algorithm is also applicable to the generation of triangular and quadrilateral global grid systems. It only needs to change the basic shape of the mesh in the base extended grid coordinate system. In essence, the algorithm truly achieves the generalized generation of multiscale meshes of each split type.

(4)

Compared with the current single-level grid generation method, the algorithm is more concise, more flexible and more operable.

While the algorithm has many advantages, there are also some shortcomings.

(1)

The algorithm only considers one triangle of a regular polyhedron and does not involve the processing of boundary trans-plane elements in terms of the global grid elements.

(2)

There is a problem of repeated rendering at the triangle boundary and the vertex of a boundary grid, which affects the efficiency of the global grid generation to some extent.

5. Conclusions

In this paper, a new algorithm for generic grid system generation based on cell replication is proposed. The hexagonal discrete grid is taken as an example to verify the algorithm, and the experiment results show that the method has good universality and is suitable to be reused in the multitype grid generation process. The core of the algorithm is “simple cell replication + effective region control”. By changing the initial transformation parameters, any shape and any type of multiscale grid can be generated, which avoids the limitations of traditional customization algorithms and achieves the goal of the unified generation of all kinds of spherical grid systems. In conclusion, this method is superior to the traditional methods in the following aspects: Firstly, in the traditional multi-scale grids generation methods, users need to customize the corresponding grids generation algorithm according to the factors of grids aperture, the prime factor of subdivision frequency, duality, eccentricity, etc. While adopting this method, users only need to replace the basic parameters and rules, and can use the same program to generate different types of grids, which is easier to implement. Secondly, with this universal generation method, due to the use of a unified benchmark, the potential links among different grid cells are highlighted. Not only grids with the same shape, but also those having different shapes only need to correspond to the nodes on the basis of the plane coordinates of Figure 4, so that different types of grids can be rapidly correlated. Such a tool will promote the development of theoretical research on communication, sharing, and interoperability of multi-type discrete grid systems.