Research and application of Boolean operation for triangular mesh model of underground space engineering—Boolean operation for triangular mesh model

This study is aiming at the problems of slow collision detection speed of bounding box, low efficiency of triangle intersection detection, and triangulation in regions where complex intersecting lines are segmented that exist in Boolean operation of triangular mesh model in underground space engineering. The orientation calculation method of underground space engineering is put forward, and the precise definition of bounding box is realized. A simple method to realize the fast intersection detection algorithm between triangles is presented, the process of transforming the projection of intersecting lines into scalar and then judging is eliminated, and the efficiency of triangle intersection detection is improved. On the basis of Delaunay triangulation robust algorithm, numbering method of intersection position is improved. Central axis region division method and hierarchical subdivision triangulation method are proposed, which can divide all the regions of mixed complex intersections in one traversal. A subdivision triangulation method by comparing the cosine value of internal tension angle is proposed, and the subdivision triangulation of polygon is optimized. The feasibility of the algorithm is verified by practical engineering application.

characteristics. It can approximate the surface of a complex body with multiple patches and is easy to handle. Therefore, it is widely used in surveying and mapping engineering, geological engineering, three-dimensional animation, CAD/ CAM, virtual reality, forward/reverse engineering, and other fields.
So far, many scholars have done researches on Boolean operations. Shiqi Ou 3 introduces a mechanism to register intersecting points in the process of grid intersection, classifies the topological relations of intersecting points, and then realizes the Boolean operation between grid models. A CAD modeling method based on inferential Boolean operation is proposed by Sun 4 to construct heterogeneous material objects. Miyazaki 5 uses Boolean operations to identify buildings in different GIS data maps. Ming Chen 6 uses hierarchical depth image technology and vertex connection propagation rules to determine the internal and external relations of triangular mesh units, and the amount of triangular mesh classification calculation was reduced. Sato, Andr Kubagawa 7 uses irregular Boolean operations to determine the collision-free region of material packaging to achieve the optimal spatial layout of irregular objects. Zhanglong Yang 8 gets the point cloud model of ray segment by discrete sampling of the grid model and transforms the three-dimensional Boolean operation of triangular meshes into one-dimensional ray segments, so as to accurately solve the intersection point and interpolate around the intersection point. Gao Yi 9 has constructed an efficient data structure of interoperability between CPU and GPU, and proposed a fast contour extraction method based on GPU to realize the Boolean operation of polygon. Although the above research of Boolean operation has made positive progress, there are still some problems, including slow speed of bounding box collision detection, low efficiency of triangle intersection detection, and complex triangulation of intersection segmentation region. To solve these problems, a faster detection method of bounding box collision and a simpler algorithm for fast intersection detection between triangles are proposed. Based on the robust Delaunay triangulation algorithm, numbering method of intersection line position is improved, and the central axis region division method and the hierarchical subdivision method of triangulation are proposed.

| ALGORITHM PRINCIPLE
The key points of Boolean operation of triangular mesh model are as follows: (a) fast collision intersection detection between triangular elements; (b) intersection judgment and generation between triangles, according to which the main triangle is divided into regions and retriangulated; and (c) accurate judgment of the triangle inside or outside another model.
The Boolean operation of triangular mesh model is as follows: 1. Orientation of the engineering body is calculated automatically according to the triangle mesh model and constructing the triangle bounding box along the strike direction. Triangle outside the bounding box is excluded and classified as an external triangle. 2. Each triangle's AABB bounding box is constructed along the strike direction, and collision detection is carried out by combining the space division of a balanced binary tree. The intersection judgment of triangles is carried out when bounding boxes intersect. 3. The intersection of space triangle is judged by the criterion of positive and negative area of the plane and whether the intersection of the triangle and plane is overlapped, and then extracting and generating intersection rings or chains of intersection triangle. 4. Triangle is divided into polygon regions by intersection chain and ring, and the regions are divided into triangles through the subdivision triangulation method again. 5. The triangles are internally and externally discriminated according to the intersection normal vector. 6. Classify the remaining undistinguishable triangles into internal or external according to the principle of adjacent triangles.
Assume that triangulation models participating in Boolean operations are A and B, respectively, triangulation A is the main triangulation model, and B is tangent triangulation model. The basic idea of a Boolean operation is to use B to segment A into two parts, including A in , inside B, and A ex , outside B. Similarly, B as the main triangle and A as tangent triangle, B is divided into B in and B ex . Boolean operations achieve intersection A∩B = A in + B in ; union A∪B = A ex + B ex ; subtraction A -B = A ex + B in ; and subtraction B -A = A in + B ex . Based on the above operational thinking, it is only necessary to study the perfect triangulation of the main triangle. Boolean operation flow of triangular mesh model is shown in Figure 1:

| Collision detection of AABB bounding box along the strike direction
In order to quickly find the intersection between the triangular meshes, that is, to realize the collision detection between geometric models, the bounding box of the geometric model is often used for collision detection in advance. The commonly used bounding box types are as follows: AABB (axis-aligned bounding box), OBB (oriented bounding box), k-DOps (k-sided discretely oriented polytopes), CBB (cylindrical bounding box), Bounding Sphere, and Convex envelope. Tightness of CBB and Bounding Sphere is poor, and Convex envelope is complex to construct. Bounding box is F I G U R E 1 Flowchart of Boolean operation among models the ideal choice. Therefore, it is preferred to choose a tight bounding box. AABB 10 is a relatively simple type of bounding box. It has simplicity intersection test and good tightness. But for the elongated object placed diagonally, the tightness is poor. OBB 11 has better tightness than AABB, but its structure and intersection test are more complex.
For underground engineering with obvious azimuth characteristics, if the AABB is constructed with the engineering trend as the axis, it can effectively compensate for the tight defects.

| Automatic calculation of engineering strike direction
At present, there are few studies on the automatic calculation of engineering strike direction by using 3D laser scanning data. Two kinds of strike direction methods are proposed through research: (a) Projecting the point cloud data onto the XOY plane, and the direction of the longest two-point connection is the strike direction of the engineering body; (b) on the basic of the triangulation model, triangle normal vector is projected on the XOY plane, and the four quadrants projection normal vectors are respectively accumulated. The vertical direction of the accumulated normal vector which has the maximum mode is the strike direction of the engineering body. As shown in picture 2, the method (1) has a case where the rectangular diagonal is larger than the rectangular side, the orientation is shifted, and the direction A1 is calculated by the method (1); in the method (2), the triangular normal vector is equivalent to the triangular cumulative area, the exposed area of the engineering side along the strike is largest, and the calculated orientation is more accurate. In Figure 2, N1 and N2 are the normal vectors of triangle 1 and triangle 2, respectively, N = N1 + N2, and A2 perpendicular to N is the strike direction of engineering.
The algorithm flowchart of the preferred method (2) is shown in Figure 3, where A is the direction vector, Z is the Z-axis normal vector, and N max is the maximum value of the normal vectors sum in the quadrant.

| AABB intersection test
The intersection test between bounding boxes is simpler than that between triangles. Firstly, AABB of tangent triangles are constructed to test the intersection, and the triangles of the main triangle outside the bounding box are quickly eliminated. Then, the triangular network is constructed by a top-down method to construct its hierarchical bounding box (AABB tree), traversing the leaf nodes in the bounding box tree of the main triangular model. Finally, the triangles contained in the leaf nodes are accurately intersected, and if they intersect, the intersection segments are extracted.
The efficiency of bounding box intersection test determines the time cost of a Boolean operation. The intersection condition of two AABB bounding boxes is that they have overlapping regions on three coordinate axes. In order to improve efficiency, BSP, 12 Q-Tree, 13 and Octree 14 can be used to construct hierarchical bounding boxes. Since the triangulation model of the underground space engineering is a hollow model, some triangles may span multiple boundaries of octree subdivision cube or quadtree subdivision rectangle, resulting in complicated division and search judgment. In order to minimize the crossover of triangle and improve the efficiency, a balanced binary tree along strike is used to construct the hierarchical bounding box.
Traditional binary tree indexes have the following disadvantages: (a) Entity can only be stored in a leaf node. The intermediate node and the root node do not store entity information, resulting in a binary deep layer and low query efficiency. (b) The same entity may be stored in multiple nodes during splitting process of binary tree, resulting in waste of Due to the uneven distribution of objects, the tree structure is extremely unbalanced, and the storage space is wasted.
By improving the binary tree, the entity is stored in the smallest rectangular node which contains it completely, not in its leaf node. Each entity is stored only once in the tree to avoid the waste of storage space. First, a full binary tree is generated to avoid redistributing memory allocation when the entity is inserted and speeding up the insertion, and finally, releasing the memory space occupied by the empty node. The improved binary tree structure is shown in Figure 4. The depth of the binary tree is set between 4 ~ 10. The bounding box intersection test flow is shown in Figure 5.

| Triangle intersection judgment
After the bounding box collision detection, intersection judgment between triangles is performed. There are several cases where two space triangles do not intersect: (a) Planes of the two triangles are parallel but not coplanar; (b) vertices of the triangle are located on the plane side of the other triangle; and (c) vertices of the triangle are located on both sides of the plane of the other triangle but do not intersect. Excluding the three cases, the two triangles intersect. If two triangles are coplanar, the problem is converted to intersect of line segments on the plane.  15 conducted a comparative study on typical spatial triangle collision detection algorithms and concluded that: Devillers algorithm is generally optimal. Tropp algorithm is most suitable for high intersection rates. Shen algorithm is more suitable for medium and low intersection rates. Although overall performance of Möller algorithm is slightly weak, it is the basis of other algorithms and is widely used.

F I G U R E 4
The common points of the above algorithms is to determine whether the three vertices of a triangle are on the same side through a plane. In the intersecting triangle, the intersection point of the triangle's edge and the plane where the other triangle is located is needed. A simpler method is implemented to realize the fast intersection detection algorithm between triangles referring to the above algorithm. For convenience, it is assumed that there are triangles TraA and TraB, vertices are A1, A2, and A3 of TraA and B1, B2, and B3 of TraB, respectively. N is the normal vector, and the plane of TraB is Π2. Sort by the hour hand, vertices of triangle are ranked by counterclockwise.

1.
Taking any vertex (B1) of triangle TraB as the starting point, the three vertex vectors V1, V2, and V3 to form B1 to triangle TraA are calculated, respectively.

2.
Calculating the cosine values cos1, cos2, and cos3 of the angle between vectors V1, V2, and V3 and normal vector N, respectively; 3. As shown in Figure 6, if cos1 and Cos2 are not greater than 0 or less than 0, A1 and A2 are on the opposite side of 2, respectively. Similarly, cos1 and cos3, and cos2 and cos3 are judged. If any two cosine values equal to 0, the two triangles are coplanar. Firstly, it is judged whether there are intersections between three sides of a triangle and three sides of another triangle. If there are intersections, the two triangles intersect. Otherwise, it is further judged whether there are intersections. To determine whether two triangles contain each other, it is necessary to determine whether the vertex of one triangle is in another triangle. According to the vector cross-product, whether the point is included in the triangle is judged. If any point is contained in another triangle, the two triangles intersect. Otherwise, the two triangles do not intersect. If the three groups of cosine values do not have an opposite sign, it is determined that the two triangles do not intersect. Otherwise, go to the next step. 4. Determine the three vertices of the triangle TraB according to the steps (a) to (c). If the three vertices are on the same side of the plane where TraA is located, the two triangles do not intersect. Otherwise, go to the next step. 5. The triangle intersection judgment and the intersection line extraction are performed according to the overlapping segments, and the method is schematically shown in Figure 7. Let intersection points of triangle TraA and plane Π2 be P1 and P2, and the intersection points of triangle TraB and plane Π1 be P3 and P4. Lines P1P2 and P3P4 are on the intersection line L of the two planes, and the triangle intersecting condition is that P1P2 and P3P4 have overlapping segments. Overlapping segments can be judged and extracted through the intersection spacing criterion. The steps are as follows: a. Calculate the distance between the intersection points, which are D12, D34, D13, D14, D23, and D24. is in the middle of P3P4, then P2P3 is the intersection of the triangle TraA and TraB and is saved to the intersection container.
The intersection judgment and extraction of overlapping line segment can be directly performed by intersection distance criterion, which eliminates the process of transforming the intersection line into scalar and then judging, and the efficiency is improved.

| Intersection generation
Triangular TraA intersects with triangular network B to generate intersecting line segments. One or more intersection line is generated by searching for the beginning and ending of the intersection segments, and triangle TraA is divided into a plurality of closed regions.
According to the relationship between intersecting lines in the triangle, intersecting lines are divided into three basic types 16 : cross-angle type, cross-border type, and ring type, which are shown in Figure 8. Complex intersection line is a combination of the three basic types.
The intersection line divides the triangle into two or more regions. The next step is to discuss how to divide and triangulate the regions accurately to determine whether the newly generated triangles belong to internal or external objects. Intersection of triangles is expressed by extending the point structure with normal vectors. The normal vector of the intersection expressed by the extension point structure is used as the normal vector of the opposite triangle. Through this normal vector, internal and external relations of the newly generated triangle can be judged.

| Triangulation method selection
In the literature, 17 taking the vertices and edges of a triangle as reference, vertices of the triangle are sorted counterclockwise. Combining the Euclidean distance between intersection point and the triangle vertex, position of the intersection points is numbered. When dividing regions, according to the location number of intersection lines, the regions are divided into three types: cross-border type, cross-angle type, and mixed type. This method has a good reference for solving the triangle internal region segmentation, but it needs to distinguish the cross-border type and the cross-angle type separately.
In the literature, 18 According to the topological relationship between the triangle edge and the intersection line, the Delaunay triangulation robust algorithm is used to realize the triangulation inside the triangle. This method does not require special regions to generate new triangles, as shown in Figure 9. This method is difficult to implement in case of complex constraints.
Through analysis and comparison, based on the above triangle internal region segmentation method, intersection line location numbering method is improved, and the central axis region division method and hierarchical subdivision triangulation method are proposed.

| Central axis region division method
The central axis region division method takes triangle center point and normal vector as the central axis, and takes the connection of the first point and central point of the triangle as the initial line. The cosine of azimuth angle between intersection point and the initial line is calculated by referring to clockwise direction of the central axis. B, the internal region, is automatically divided by sorting the cosine value. To illustrate the algorithm, the case of mixed complex intersecting line division regions is selected. As shown in Figure 10, which includes cross-edge type, cross-angle type, and cross-line common point type, all the regions can be automatically partitioned by the central axis division method at once.
Cosine calculation of the azimuth angle in the axis area division method is illustrated in Figure 11. In the figure, OA is the initial line. According to the counterclockwise order, the orientation is divided into four areas: Quad1, Quad2, Quad3, and Quad4. Points B, C, D, and E are respectively located in the four areas, and the azimuth cosines are cos B, cos C, cos D, and cos E, respectively. Steps of the central axis division method are as follows. Figure 12 is the method illustration.  otherwise, the intersection point is reversely sorted and updated. Through this step, the order of all the intersection points is adjusted. 2. Store the vertices A1, A2, and A3 as intersection lines in the intersection container. The intersection lines are sorted from the largest to the smallest according to azimuth angle of its first point. If the first points of the two intersection lines coincide, the larger the cosine value of the tail point is, the higher the ranking is. 3. Extract regional boundaries: By Figure 13 Similarly, Q 3 , Q 4 , and Q 5 are divided. When A 3 as the initial boundary, there is no point whose azimuth cosine is less than that of A 3 , so the region is not divided. All the regions of mixed complex intersection can be completely divided at once.

| Subdivision triangulation of region polygon
The divided areas are plane polygon, and there are many mature methods for triangulation of planar polygon, 19 such as ear-cutting method, node connection method, template method, topological subdivision method, grid method, geometric decomposition method, and Delaunay triangulation algorithm. After triangulation, the fewer elongated triangles, the more acute triangles, and the more uniform the triangular mesh distribution, the better the triangulation results. In order to make the subdivided triangle get close to the optimal triangle, a subdivided triangulation method by comparing the cosine values of the internal tension angle is proposed.
The method is as follows: The tension angles of each edge with all points are calculated separately, and the edge and point with the largest tension angle are extracted and formed a triangle. So the polygon is subdivided into two parts, until all polygons are subdivided into triangles. The method is shown in Figure 14.

| Internal and external relation judgment of triangle
For entities involved in a Boolean operation, triangular meshes can be classified into three categories: within, outside, and intersecting with another entity. The intersecting triangular meshes are divided into subdividing triangles by regional division. It is only necessary to judge the internal and external relationship between the subdividing triangles and another entity. The external triangles are judged by whether it is co-edge with the seed triangle. The internal triangles are judged by whether the angle between the regional normal vector and the intersection vector is less than 90 degrees.

| ALGORITHM APPLICATION
Two measured goafs, one roadway and one ore pass, are selected as objects. The number of triangles of goafs is 108 580   Figure 15 is triangulation models of objects. Figure 16 is a schematic diagram of the intersecting position of objects. Figure 17 is the calculation result of the Boolean operation. Table 1 compares the time consumed by different algorithms among the results of Figure 17.

| DISCUSSION
For collision test between models, the traditional AABB bounding box has poor compactness for slender objects placed along oblique diagonal direction. In this paper, the strike direction of slender objects is obtained firstly. Then, the strike direction is deduced to construct their bounding box, which solves the problem of poor compactness. Compared with the reference 20, for slender objects, the bounding box constructed by the improved method is significantly smaller than that of traditional method. The efficiency of collision detection is improved.
For the extraction of intersection lines, intersection judgment and extraction can be directly performed by the intersection spacing criterion. The process of transforming the projection of intersecting lines into scalar and then judging is eliminated. Reference 21 proposes a BSP to B-rep conversion algorithm which can accurately generate boundaries with sharp and small features. Compared with the reference 21, our method is faster in extracting intersection lines, but reference 21 is better in dealing with sharp and minute feature structures.
When reconstructing the surface in the intersecting region, our method can accurately capture the geometric details, while the uncrossed region remains the same. Unlike the previous works of adaptive contouring (reference 22,23) which reconstruct the whole mesh surface, we take a similar strategy as references 24 and 25. Reference 26 generates a new triangulation only in a close vicinity around the intersections of the input meshes and thus preserve as much of the original mesh structure as possible. It has very high accuracy, but it reduces the efficiency of intersection detection. Reference 27 introduces a novel solid modeling framework taking advantage of the architecture of parallel computing on modern graphics hardware. It has a relatively high computational efficiency. However, this method is memory occupied. Reference 28 is designed for sculpting, which is used for BSP. However, most of the systems are likely to use B-reps, not BSP trees. The described input and output paths provide a rudimentary sketch of the representation conversions needed.
As can be seen from Table 1, when A∪B is carried out, collision tests are mainly carried out, and the retriangulation area is very small. The advantage of our method in collision detection is brought into play, and the calculation speed is relatively fast. When A∩B is carried out, the area of subdivision enlarges, reference 25 has obvious advantages, and the calculation speed is faster. When A-B and B-A are carried out, our method is close to reference 25, and the efficiency is relatively high. Reference 27 has high accuracy in dealing with details and reduces the computational efficiency. In reference 28, the calculation time is increased when BSP is converted to B-reps. So the methods in references 27 and 28 take a relatively larger computing time.

1.
A method for calculating the orientation of the underground space engineering is proposed. Based on the established triangular mesh model, the normal vectors of four quadrants projected on the XOY plane are accumulated separately. The direction perpendicular to the accumulated normal vector which has the maximum vector mode is the orientation of the engineering.

2.
A simple method to realize the fast intersection detection algorithm between triangles is presented. Intersection judgment and extraction can be directly performed by the intersection spacing criterion. The process of transforming the projection of intersecting lines into scalar and then judging is eliminated, and the efficiency of intersection detection is improved. 3. On the basis of the Delaunay triangulation robust algorithm, numbering method of intersection position is improved. Central axis region division method and hierarchical subdivision triangulation method are proposed, which can divide all the regions of mixed complex intersections in one traversal. 4. In order to make the subdivided triangle get close to the optimal triangle, a subdivided triangulation method by comparing the cosine values of the internal tension angle is proposed. Polygon subdivision triangulation is realized.