A GPU-based method for determining whether a point is within a spherical polygon

Abstract

The present invention discloses a method for determining whether a point is located within a spherical polygon based on a GPU. The method of the present invention is as follows: 1) Initializing an adaptive hexagonal hierarchical grid structure (abbreviated as AHHG structure); 2) Starting from the first layer of grid, the AHHG structure is created layer by layer by iteratively and parallelly subdividing the cells containing polygon edges in the current layer of grid; when the cutoff subdivision condition is met, the subdivision is stopped, and then the edges of the polygon are distributed in parallel to the bottom layer of grid cells to complete the creation of the AHHG structure; 3) The leaf cell attributes of the bottom layer of grid are determined in parallel, and then the leaf cell attributes of other layers of grid are determined in parallel from bottom to top; 4) For a given point, according to the AHHG structure created in the above steps, the leaf cell into which the measured point falls is searched in parallel, and whether the point is located within the spherical polygon is determined based on the leaf cell attributes. The present invention can efficiently complete the determination of whether a point is located within a dynamically changing spherical polygon.

Description

Technical Field

[0001] The invention belongs to the technical field of geographic information and computer graphics, and relates to a method for determining whether a point is located within a spherical polygon based on a GPU. Background Art

[0002] Determining whether a point is located on a spherical polygon is a basic operation in spherical information processing. Compared with the operation of determining whether a point is located inside a planar polygon, this problem has been less studied. Existing methods are divided into two categories: conversion to plane space processing and direct processing in spherical space. The first type of method usually first projects the spherical polygon onto a plane using some projection method, and then uses a method to determine whether the point is located inside the planar polygon (National Aeronautics and Space Administration, 2020. General mission analysis tool (gmat). Available from: https: / / sourceforge.net / projects / gmat / files / GMAT / GMAT-R2020a / .). These methods are easy to implement, but because projecting a spherical arc segment onto a straight line on a plane will cause deformation, the judgment result may be incorrect.

[0003] Regarding the second type of method, Bevis et al. (Bevis, M. and Chate Lain, J., 1989. Locating a point on a spherical surface relative to an aspherical polygon of an arbitrary shape. Mathematical Geology, 21, 811–828.) proposed performing raycrossing calculations directly on a spherical surface. This method first finds a point X inside the spherical polygon, then connects it to the measured point Q to form a spherical arc segment QX. The number of intersections between this arc segment and the polygon edge is calculated. If there is an odd number of intersections, Q is inside the polygon; otherwise, it is outside the polygon. This method requires checking all polygon edges and requires a large number of trigonometric functions to determine whether two arc segments intersect, making it inefficient.

[0004] Another approach (Li, J., Zhang, H., and Wang, W., 2017. Fastand robust point-in-spherical-polygon tests using multilevel spherical grids. In: J. Chang, J. J. Zhang, N. Magnenat Thalmann, S. M. Hu, R. Tong, and W. Wang, eds. Next Generation Computer Animation Techniques, Cham. Springer International Publishing, 56–66) proposes creating a latitude and longitude grid by evenly dividing the longitude and latitude grids and precalculating the attributes of the grid cell center points (inside / outside the polygon). This allows for rapid identification of the cell in which a point is detected. An arc is then formed connecting the point to the cell center point, and the detection result is then determined using local ray-crossing calculations along this arc. This approach performs arc intersection calculations in Cartesian coordinate space, thus avoiding trigonometric calculations and resulting in high speed. However, the size of its grid cells changes with latitude, which may cause some cells to contain a large number of polygon edges, resulting in excessive intersection calculations and reduced judgment efficiency.

[0005] Ketzner et al. (Kerzner et al. 2022) improved the method of (Bevis, M. and Chate Lain, J., 1989. Locating a point on a spherical surface relative to an aspherical polygon of an arbitrary shape. Mathematical Geology, 21, 811–828.). They adopted the idea of ​​(Li, J., Zhang, H., and Wang, W., 2017. Fastand robust point-in-spherical-polygon tests using multilevel spherical grids. In: J. Chang, J. J. Zhang, N. Magnenat Thalmann, S. M. Hu, R. Tong and W. Wang, eds. Next Generation Computer Animation Techniques, Cham. Springer International Publishing, 56–66.) to determine whether two arc segments intersect in Cartesian space, thereby significantly improving the efficiency of the method of (Bevis, M. and Chate Lain, J., 1989. Locating a point on a spherical surface relative to an aspherical polygon of an arbitrary shape. Mathematical Geology, 21, 811–828.) Based on this, they also proposed a method that requires pre-computation to handle the situation where a large number of measured points need to be determined for the same polygon. It first identifies a known reference point within the polygon and connects it to each vertex of the polygon, forming an uneven partition of the spherical longitude, called a lune. The polygon's edges are then assigned to these lunes, and a binary search tree is created to manage these lunes. During the decision process, the lune containing the measured point is first found through a binary search tree. Then, local ray-crossing is performed based on the polygon edges within that lune to determine the decision result. Because this method partitions the space based on polygon vertices only in the longitude direction, some lunes may contain too many edges, increasing computational overhead and affecting decision efficiency.

[0006] In 2022, Li Jing et al. (Li, J. and Wang, W., 2022. Adaptive hexagonal hierarchical grid for point-in-spherical-polygon tests. Journal of Software, 33, 3485–3497.) proposed a determination method based on an adaptive hexagonal hierarchical grid, referred to as the AHHG method. This method first creates a regular icosahedron structure that fits the sphere and generates a first-level global hexagonal grid. Then, using a 1-to-4 subdivision method, cells containing polygon edges are adaptively subdivided based on the polygon edges, forming the next level of subdivided grids containing hexagons, pentagons, and quadrilaterals. This subdivision is iteratively performed layer by layer, forming a multi-layer adaptive hexagonal grid structure. The leaf cell center attributes (inside / outside the polygon) in the AHHG structure are then precalculated. When determining a measured point, the leaf cell containing the point is first quickly found. Then, the measured point attributes are determined using a local ray-crossing method based on the leaf cell center attributes. Because hexagonal grid cell size does not vary with latitude, it effectively addresses the uneven cell size issues encountered in methods based on latitude and longitude grids (Li, J., Zhang, H., and Wang, W., 2017. Fastand robust point-in-spherical-polygon tests using multilevel spherical grids. In: J. Chang, J. J. Zhang, N. Magnenat Thalmann, S. M. Hu, R. Tong, and W. Wang, eds. Next Generation Computer Animation Techniques, Cham. Springer International Publishing, 56–66). However, despite the high efficiency of the AHHG method, creating the AHHG structure is time-consuming. This is due to the complex structure for recording the adjacency relationships between cells at each level and the non-trivial process of pre-calculating leaf cell attributes. As we all know, there are many applications that require real-time determination of dynamic polygons. However, when a polygon changes, the AHHG structure must be recreated. This time-consuming creation process prevents the AHHG method from achieving real-time determination. Summary of the Invention

[0007] In order to solve the above technical problems, the present invention proposes a GPU-based method for determining whether a point is located within a spherical polygon. The method performs preprocessing (including creating an AHHG structure and determining leaf unit attributes) and determination operations in parallel on the GPU.

[0008] The technical solution of the present invention is:

[0009] A method for determining whether a point is located within a spherical polygon based on a GPU, comprising the following steps:

[0010] 1) Pre-allocate storage space for an adaptive hexagonal hierarchical grid structure, referred to as an AHHG structure, on the GPU and CPU sides; create the first-level grid of the AHHG structure in a serial manner on the CPU side; after the creation is completed, copy the first-level grid data from the CPU side to the GPU side.

[0011] 2) Starting from the first layer of mesh, the AHHG structure is created layer by layer by iteratively and parallelly subdividing the cells in the current layer of mesh that contain polygon edges, referred to as non-empty cells. When the cutoff subdivision condition is met, the subdivision is stopped, and then the polygon edges are distributed in parallel to the bottom layer of mesh cells, finally completing the creation of the AHHG structure.

[0012] 3) Determine the attributes of all leaf cells of the bottom grid of the AHHG structure in parallel, and then determine the attributes of the leaf cells of each grid layer in parallel starting from the grid layer above the bottom grid to the first grid layer.

[0013] 4) Based on the AHHG structure, determine in parallel whether the measured point is located within the spherical polygon.

[0014] According to another aspect of the present invention, a computer storage medium is provided, wherein the storage medium stores at least one executable instruction, and the executable instruction enables a processor to execute operations corresponding to the estimation method described in the present invention.

[0015] The beneficial effects of the present invention are as follows:

[0016] The AHHG method can effectively handle the problem of determining whether a point is located within a polygon by creating an AHHG structure. However, the expensive creation overhead of the AHHG structure hinders its application in practice, especially when the polygon changes dynamically. The present invention effectively avoids the problem of inconsistent data access during parallel operations by rationally dividing the task into several parallel subtasks that are executed serially in sequence. In this way, the parallel computing capabilities of the GPU can be fully developed, effectively accelerating the creation and processing of the AHHG structure. Therefore, the present invention can efficiently complete the determination of whether a point is located within a dynamically changing spherical polygon. Experiments show that for a dynamic polygon containing 100,000 edges, the present invention can determine 1,000,000 measured points in real time on a PC, which is impossible with existing methods. Thanks to this, as a basic operation, the present invention can improve many problems in global information processing. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] In order to more clearly illustrate the technical solutions of the embodiments of the present application, the following is a brief introduction to the drawings required for use in the embodiments of the present application. For those skilled in the art, other drawings can be obtained based on these drawings without any creative work.

[0018] Figure 1 is a flow chart of a GPU-based method for determining whether a point is within a spherical polygon;

[0019] Figure 2 This is a diagram illustrating the storage structure on the GPU side and the CPU side;

[0020] Figure 3 It is a schematic diagram for calculating the number of new geometric elements of the boundary subunit;

[0021] Figure 4 This is a diagram of how the polygon edges are assigned to the bottom-level grid;

[0022] Figure 5 It is a schematic diagram for determining the properties of empty cells adjacent to subdivided cells. DETAILED DESCRIPTION

[0023] The following describes a method for determining whether a point is located within a spherical polygon based on a GPU provided by the present invention in conjunction with the accompanying drawings and embodiments. The flowchart of the method for determining whether a point is located within a spherical polygon based on a GPU is shown as follows: Figure 1 As shown, an introduction is given.

[0024] The present invention comprises four phases: the first is an initialization phase, which involves pre-allocating a storage structure and serially creating the first-level AHHG grid on the CPU and copying it to the GPU. The second phase involves creating the remaining levels of the AHHG structure in parallel on the GPU. The third phase involves determining the leaf cell attributes of the AHHG structure in parallel on the GPU. The fourth phase involves determining whether a measured point lies within a spherical polygon based on the AHHG structure in parallel on the GPU.

[0025] Before creating an AHHG, we read polygon data from a polygon file and copy it to the GPU. Polygon data consists of a polygon vertex structure array (PV) and a polygon edge structure array (PE). PV stores the Cartesian coordinates of the polygon vertices. PE stores the index of each edge's starting point and ending point in PV, in the order of the polygon's edge chain. For clarity, we will refer to AHHG edges as mesh edges and spherical polygon edges as polygon edges. The specific steps of the present invention are detailed below.

[0026] (1) Initialization

[0027] (1.1) Pre-allocate the storage space of the AHHG structure on the GPU and CPU

[0028] In this step, storage space for the GPU-side AHHG structure is allocated, including: a vertex structure array V, a mesh edge structure array B, a cell structure array C, and a polygon edge index array R. The vertices of all levels are stored in array V in order of level, that is, the vertices of the first mesh level are stored first, then the vertices of the second mesh level, and so on. Similarly, the mesh edges of all levels are stored in array B in order of level, and the cells of all levels are stored in array C in order of level. For elements in arrays V, B, and C, we can use indexes or level-relative indexes to search. The index refers to its serial number in the entire array element, while the level-relative index refers to its serial number in the array element of this level. The vertex structure stores the vertex information of the mesh cell, that is, the Cartesian coordinates of the vertex. The mesh edge structure stores the information of the directed edges of the mesh cell; the cell structure stores the mesh cell information. The polygon edge index array R stores the index of the polygon edge contained in the cell in the polygon edge array PE, which is called the edge pointer. The CPU's storage space includes a grid structure array G, which records the starting storage location and number of geometric elements (vertices, grid edges, and cells) of all levels in the corresponding array. Figure 2 This is a diagram illustrating the storage structure on the GPU side and the CPU side.

[0029] (1.2) Create the first layer grid

[0030] In this step, we create the first-level mesh of AHHG in a serial manner on the CPU side. The creation method is the same as that in the literature (Li, J. and Wang, W., 2022. Adaptive hexagonal hierarchical grid for point-in-spherical-polygon tests. Journal of Software, 33, 3485–3497.). The first-level mesh created covers the entire sphere and includes 80 vertices, 240 directed edges, and 42 cells (12 pentagonal cells and 30 hexagonal cells). These data are first stored in a temporary array on the CPU side. After creation, they are copied from the CPU side to the pre-allocated V, B, and C arrays on the GPU side. Then, the first-level mesh is used as the current level mesh to facilitate subsequent iterative subdivision. For convenience of description, we use l to represent the current grid level, and for the first-level mesh, l = 0.

[0031] (2) Create the entire grid structure

[0032] In this step, the present invention begins with the first-level mesh and creates meshes layer by layer from top to bottom by parallel iterative subdivision of the current level. When the subdivision cutoff condition is met, subdivision stops, and the polygon edges are then distributed in parallel to the cells of the bottom-level mesh, ultimately completing the construction of the AHHG structure. The specific steps are as follows.

[0033] (2.1) Parallel marking of non-empty cells

[0034] We first mark all cells in the current layer that need to be subdivided, that is, non-empty cells containing polygon edges. Subsequent subdivision operations are performed only on these non-empty cells. Here, we track and mark the grid cells in the current layer that the polygon edges traverse. We assign a thread to each polygon edge. Each thread first locates the current layer cell where the polygon edge starts, then sequentially searches for the cell edge it exits, and finds the next cell to enter based on the adjacency relationship between cell edges, until it reaches the cell where the polygon edge ends. During this process, marks are set in all traversed cell structures. Since multiple edges may enter the same cell, multiple threads may write to the same cell structure. However, since the cell is marked, there is no need to handle write conflicts. These marked cells are the cells that need to be subdivided.

[0035] (2.2) Determine whether to end the subdivision

[0036] Count the total number of marked cells, then calculate the average number of polygon edges γ per non-empty cell in the current mesh layer (i.e., layer l), where γ = total number of polygon edges / total number of marked cells. There are two conditions for stopping the creation of the next mesh layer: first, the number of layers has reached the user-defined maximum number of layers; second, γ is less than a given threshold. If either of these conditions is met, subdivision ends and proceeds to step (2.6). Otherwise, proceed to step (2.3).

[0037] (2.3) Generate central subunit

[0038] In this step, we create the central sub-cell of the current layer (i.e., layer l), that is, we generate new geometric elements belonging to the central sub-cell, including new vertices, new mesh edges, and new cells, and write them to the corresponding positions in the l+1 layer in arrays V, B, and C. This step is divided into three stages: calculating the number of new geometric elements, calculating the write positions of the new geometric elements based on the number, and writing the new geometric elements according to the write positions.

[0039] (2.3.1) Calculate the number of new geometric elements

[0040] In the AHHG structure, each non-empty cell, after being subdivided, will generate a central sub-cell, new vertices with the same number of vertices as the sub-cell, and new mesh edges with the same number of vertices as the sub-cell. To this end, we use two buffer arrays, buf0 and buf1, to record the number of new vertices and new cells for each central sub-cell, respectively. Since the number of new mesh edges is the same as the number of new vertices, buf0 can be shared. We assign a thread to each current layer cell. For a thread i, if its corresponding cell is a marked cell, let buf0[i] be the number of vertices of the cell, and buf1[i] be 1; otherwise, the values ​​of buf0[i] and buf1[i] are both 0.

[0041] (2.3.2) Calculate the write position of the new geometric element

[0042] Calculate the preorder sum of buf0 and buf1 in parallel and store the result back in buf0 and buf1. At this point, buf0 stores the relative index of each new vertex in V, which is also the relative index of each new mesh edge in B, and buf1 stores the relative index of each new cell in C.

[0043] (2.3.3) Writing new geometric elements

[0044] We assign a thread to each current layer unit. If the current unit is a marked unit, new geometric element data is generated and written to the corresponding positions in the V, B, and C arrays according to the values ​​in buf0 and buf1. The specific process is to set the unit corresponding to thread i to c. If c is a marked unit and it contains k mesh edges, then use the method in (Li, J. and Wang, W., 2022. Adaptive hexagonal hierarchical grid for point-in-spherical-polygon tests. Journal of Software, 33, 3485–3497.) to generate a total of k new vertices, k new mesh edges, and 1 new unit for the central subunit of c; let the next layer of the grid of the current layer be The starting storage position of the vertices of the next layer (i.e., the l+1 layer) in V is vstart, then k new vertex data are written to V in sequence starting from vstart+buf0[i]; let the starting storage position of the mesh edges of the next layer of the current layer (i.e., the l+1 layer) in B be bstart, then k new mesh edges are written to B in sequence starting from bstart+buf0[i]; let the starting storage position of the cells of the next layer of the current layer (i.e., the l+1 layer) in C be cstart, then new cells are written to the position of cstart+buf1[i] in C.

[0045] (2.4) Generate boundary subunits

[0046] In this step, we create the boundary sub-cell of the current layer (i.e., layer l), that is, generate new geometric elements belonging to the boundary sub-cell, including new mesh edges and new cells, and write them to the corresponding positions in arrays B and C. This step is also divided into three stages: calculating the number of new geometric elements, calculating the write positions of the new geometric elements based on the number, and writing the new geometric elements according to the write positions.

[0047] (2.4.1) Calculate the number of new geometric elements

[0048] In the AHHG structure, after each non-empty cell is subdivided, each of its edges will generate a corresponding boundary subcell (except the longest edge of the quadrilateral cell), and it is adjacent to the center subcell of the subdivided cell. Therefore, we can create adjacent boundary subcells starting from the edges of the center subcell. This is more parallel than creating all its boundary subcells sequentially starting from the current layer cell. To this end, we assign a thread to each edge of the center subcell generated in step (2.3) to parallelly calculate the number of new geometric elements from the boundary subcell. Here, we still use two buffer arrays, buf0 and buf1, to record the number of new mesh edges and new cells from each boundary subcell, respectively. For a thread i, if its corresponding edge b belongs to the center subcell of the quadrilateral and is the longest edge in the cell, then this edge has no corresponding boundary subcell, and let buf0[i] = 0 and buf1[i] = 0. Otherwise, search for the upper unit of the central sub-unit to which b belongs (called the parent unit), the edge in the parent unit that is parallel to b (called the parent edge), the partner edge of the parent edge (called the parent partner edge), and the unit to which the parent partner edge belongs (called the adjacent parent unit). If the adjacent parent unit does not exist or the adjacent parent unit has no central sub-unit, the number of edges of the boundary sub-unit is 4, let buf0[i]=4, buf1[i]=1. If the adjacent parent unit exists and has a central sub-unit, compare the unit indexes of the parent unit and the adjacent parent unit. If the index of the parent unit is larger, the boundary sub-unit is a hexagonal unit, let buf0[i]=6, buf1[i]=1. Figure 3 For example, the mesh edge V3V4 belongs to the central subunit The parent unit is The edge parallel to edge V3V4, that is, the parent edge of V3V4 is V0V1, the parent partner edge of V3V4 is V1V0, and the adjacent parent unit of V3V4 is unit And unit There is a central sub-unit, so the number of edges of the boundary sub-unit generated corresponding to V3V4 is 6.

[0049] (2.4.2) Calculate the write position of the new geometric element

[0050] New geometric elements from the boundary subcells are written after the new geometric elements from the center subcells. To do this, we first set bstart = the index of the last new mesh edge of the center subcell of the current layer (i.e., layer l) + 1, and cstart = the index of the last new cell of the center subcell of the current layer (i.e., layer l) + 1. Then, we calculate the pre-order sum of buf0 and buf1 in parallel and store the results back in buf0 and buf1, respectively. At this point, buf0 stores the relative index of each new mesh edge in B starting from bstart, and buf1 stores the relative index of each new cell in C starting from cstart.

[0051] (2.4.3) Writing new geometric elements

[0052] We again assign a thread to each edge of the central subcell generated in step (2.3) and repeat the process of step (2.4.1), but replace the counting operation with writing new geometric elements. For a thread i, if its corresponding edge b belongs to a quadrilateral central subcell and is the longest edge in the cell, then this edge has no corresponding border subcell and no operation is performed. Otherwise, the parent cell of the central subcell to which b belongs, b's parent edge, b's parent partner edge, and b's adjacent parent cell are searched in sequence. If the adjacent parent cell does not exist or has no central child cell, a quadrilateral boundary child cell is generated according to the method in (Li, J. and Wang, W., 2022. Adaptive hexagonal hierarchical grid for point-in-spherical-polygon tests. Journal of Software, 33, 3485–3497.), and the four new grid edges are sequentially written to B starting from bstart+buf0[i], and a new cell is written to C at cstart+buf1[i]. If the adjacent parent cell exists and has a central child cell, the cell indexes of the parent cell and the adjacent parent cell are compared. If the index of the parent cell is larger, a hexagonal boundary subcell is generated according to the method in (Li, J. and Wang, W., 2022. Adaptive hexagonal hierarchical grid for point-in-spherical-polygon tests. Journal of Software, 33, 3485–3497.), and the six new grid edges are written sequentially to the position starting from bstart+buf0[i] in B, and one new cell is written to the position cstart+buf1[i] in C.

[0053] (2.5) Set the adjacency relationship between newly created subunits

[0054] In this step, we set the adjacency relationship between each newly created sub-unit in parallel. The adjacency relationship between sub-units is achieved by setting the partner edge relationship, which can be achieved with the help of the topological relationship with the upper-level unit. In order to facilitate the search for the upper-level unit, we choose to process the sub-units in parallel. Since the boundary sub-unit must be adjacent to the central sub-unit, we only need to process one of the two. Since there are more boundary sub-units than central sub-units, we choose to process the boundary sub-units in parallel to obtain a higher degree of parallelism. That is, we assign a thread to each boundary sub-unit, and for each boundary sub-unit a, process each of its mesh edges in turn. For each mesh edge u of unit a, first find the parent unit a' of unit a, then find the mesh edge u1' in unit a' whose corresponding boundary sub-unit is a, then find the mesh edge u2' in unit a' connected to the mesh edge u1', and then find the boundary sub-unit b corresponding to the mesh edge u2', find the mesh edge w with the same vertex as the mesh edge u from unit b, and set the mesh edge u and the mesh edge w as partner edges to each other. With Figure 3 For example, let us be the boundary subunit To find the sub-units adjacent to it through its edge V1V4, we first find the unit The parent unit Then find the unit The corresponding boundary subunit is Find the edge V0V1, then find the edge V1V2 connected to V0V1, and then find the boundary subunit corresponding to the edge V1V2 Finally found the unit The middle vertex is the edge V4V1 between V1 and V4. Set V4V1 and V1V4 as partner edges to complete the boundary subunit. and Setting of neighbor relationships.

[0055] After completing the setting of the adjacent relationship, we take the layer obtained by subdivision as the current layer, that is, let l = l + 1, and then return to step (2.1) to start the next round of subdivision.

[0056] (2.6) Assign the polygon edges to the bottom grid

[0057] After completing the iterative subdivision, we also need to assign the polygon edges to the bottom-level cells to finally complete the creation of the AHHG structure. This step searches in parallel for the polygon edges contained in each non-empty cell of the bottom-level grid, and stores their indexes in the polygon edge array, called edge pointers, in the array R in the order of the indexes of their corresponding non-empty cells. The starting storage position and number of polygon edge pointers in R are recorded in the cell. The essence of this step is to find a (cell index, edge pointer) pair, where the cell index is the index of the bottom-level grid cell containing the polygon edge, and the edge pointer is the index of the polygon edge located in the cell. If the edges of the polygons that intersect with the cell are directly searched for, each cell needs to traverse all polygon edges. When there are many polygon edges, the overhead is very high. For this reason, we still adopt the method of starting from the polygon edges, that is, tracking the polygon edges, and all cells that are traversed are the cells that intersect with the edge. The result obtained by this method is a (edge ​​pointer, cell index) pair. Since the (edge ​​pointer, cell index) pair is equivalent to the (cell index, edge pointer) pair, we can use the buffer array to convert the (edge ​​pointer, cell index) pair into the (cell index, edge pointer) pair.

[0058] exist Figure 4 In (a), the three edges of a spherical polygon intersect with a total of 7 grid cells, corresponding to 10 (edge ​​pointer, cell index) pairs, arranged in order of polygon edge index. After the transformation, the content of each pair remains unchanged, but is changed to be arranged in order of cell index. The edge pointer sequence is stored in the array R, and the starting position and number of the edge pointer sequence of each cell in R are stored in the cell structure. For example, the cell The starting position of is 0 and the number is 2. The starting position of unit ⑥ is 8 and the number is 2. For clarity, Figure 4 The cell index in is represented by a circled number.

[0059] In summary, we divide this step into two steps: obtaining (edge ​​pointer, cell index) pairs and converting them into (cell index, edge pointer) pairs.

[0060] (2.6.1) Get (edge ​​pointer, cell index) pair

[0061] We set up four temporary integer arrays, buf0, buf1, buf2, and buf3, for buffer storage and initialize them to all zeros. First, we parallelize the number of cells traversed by each polygon edge and write it into buf0. Then, we calculate the pre-order sum of buf0 and, based on the value of buf0, write the edge pointers in the (edge ​​pointer, cell index) pair into buf1 and the cell indexes in the (edge ​​pointer, cell index) pair into buf2.

[0062] For example, Figure 4In (b), for the polygon edge with index 0, buf0[0] is equal to 0, so the cell it passes through ①,② are written into the 3 positions starting from 0 in buf1 in sequence, and the edge pointer,0, is written into the 3 positions starting from 0 in buf2.

[0063] (2.6.2) is converted into (cell index, polygon edge) pairs

[0064] First, based on the value of buf1, the number of polygon edges contained in each cell is calculated in parallel and written into the array buf3. Then the pre-order of buf3 is calculated and written into buf0. Finally, based on the values ​​of buf1, buf2, and buf3, the edge pointers are written into the R array in parallel with the cell index as the main order. At this point, R stores the edges of the polygons contained in the cells in sequence with the cell index as the main order. For example, in 4(b), for thread 0, according to buf1[0] and buf2[0], it can be seen that the (edge ​​pointer, cell index) pair corresponding to thread 0 is ( 0), buf3[0] is in the storage unit The edge pointer of is at the starting storage position in R, so the edge pointer 0 is written to the position of buf0[0] in R (i.e. R[0]). For thread 9, the corresponding edge pointer pair is ( 2), edge 2 is also a unit Therefore, edge 2 is written to the position after buf0[0] in R (i.e., R[1]). Finally, based on the values ​​of buf3 and buf0, the starting position and number of edge pointers in each cell in R are set in parallel. At this point, the edges of the polygon are assigned to each cell, and the entire AHHG structure is created.

[0065] (3) Determine leaf unit attributes

[0066] After the AHHG structure is created, this step determines the attributes of all leaf cells in the AHHG structure relative to the polygon (inside / outside the polygon) to accelerate the determination of whether a point is within the spherical polygon. Leaf cells are divided into two types: empty leaf cells, which do not contain polygon edges, and non-empty leaf cells, which contain polygon edges. Empty leaf cells have clear attributes: inside or outside the polygon. However, since non-empty leaf cells occupy both inside and outside the polygon, their classification is unclear. To this end, we precompute the internal and external attributes of their center points and use them as cell attributes. The center point attributes serve as important information for determining the attributes of empty leaf cells and measured points. This step determines all leaf cells (including empty and non-empty cells) in parallel from the bottom up. Determining each grid cell requires at most three parallel processing operations for all cells, resulting in high efficiency. The entire process is divided into two steps: determining the attributes of leaf cells in the bottom-level grid and determining the attributes of leaf cells in non-bottom-level grids.

[0067] (3.1) Determine the leaf unit attributes of the bottom grid

[0068] The leaf cells of the bottom-level grid are divided into three categories: non-empty leaf cells, cells adjacent to non-empty cells (called directly empty cells), and cells adjacent only to empty cells (called indirect empty cells). Due to the characteristics of the AHHG structure, the attribute determination of these three types of cells is dependent on each other: the determination of directly empty cells depends on non-empty cells, and the determination of indirect empty cells depends on directly empty cells. To this end, in this step, we perform three parallel processing operations on all cells, each time determining a cell type, as follows.

[0069] First, the attributes of the center points of all non-empty cells in the bottom-level grid are determined in parallel and used as the non-empty cell attributes. Specifically, a thread is assigned to each bottom-level cell. If the cell corresponding to the thread is non-empty, the center point attributes are set using the method described in (Li, J. and Wang, W., 2022. Adaptive hexagonal hierarchical grid for point-in-spherical-polygon tests. Journal of Software, 33, 3485–3497.). This involves finding the edge of the polygon closest to the center point within the cell and determining whether the center point lies to the left of the plane formed by that edge and the sphere center. If so, the center point is inside the polygon; otherwise, it is outside the polygon.

[0070] Then, the attributes of directly empty cells in the bottom-level grid are determined in parallel. Specifically, a thread is assigned to each bottom-level cell. If the cell corresponding to the thread is a directly empty cell with no attributes set, the cell attributes are determined according to the method described in (Li, J. and Wang, W., 2022. Adaptive hexagonal hierarchical grid for point-in-spherical-polygon tests. Journal of Software, 33, 3485–3497.). Specifically, the midpoint of an edge adjacent to a non-empty cell is taken, and a line is drawn connecting this midpoint with the center point of the non-empty cell. The number of intersections between this line segment and the polygon edge in the non-empty cell is calculated. If the number of intersections is even, the attributes of the directly empty cell are the same as those of the center point of the non-empty cell; otherwise, they are different.

[0071] Finally, the attributes of the non-directly empty cells in the bottom-level grid are determined in parallel. This means that a thread is assigned to each bottom-level cell. If the cell corresponding to the thread is a non-directly empty cell with no attributes set, the attributes of the adjacent directly empty cell are used as the attributes of the non-directly empty cell.

[0072] (3.2) Determine the leaf unit attributes of non-bottom-level grids

[0073] After completing the classification of the bottom-level grid cells, we then proceed from bottom to top to determine the leaf cells of the non-bottom-level grid. Leaf cells of the non-bottom-level grid are divided into two categories: empty cells adjacent to subdivided cells and empty cells not adjacent to subdivided cells. Due to the characteristics of the AHHG structure, the attributes of empty cells not adjacent to subdivided cells depend on the attributes of empty cells adjacent to subdivided cells. To this end, for each layer of the non-bottom-level grid, we perform two parallel processing operations on all cells, determining one type of cell each time. The specific steps are as follows.

[0074] Starting from the bottom grid to the first grid (i.e., the 0th grid), perform the following operations on each grid layer.

[0075] First, if the previous layer of grid of the current layer grid does not exist, then end; otherwise, the previous layer of grid of the current layer grid is used as the current layer grid.

[0076] Then, we parallelly determine the properties of the empty cells adjacent to the subdivided cells in the current layer grid. That is, we assign a thread to each current layer cell. If the cell corresponding to the thread is an empty cell adjacent to the subdivided cell with no attributes set, we determine the cell properties as follows. Figure 5 For example, an empty cell and subdivided units adjacent, for them The shared edge V0V1 in the image must correspond to a quadrilateral element, and V0V1 is its longest edge. Since the midpoint U of V0V1 is the longest edge of the element have the same attributes, so we can just determine the attributes of U. Here, there are three possibilities for quadrilateral cells: empty cells, non-empty cells, and subdivided cells. If the quadrilateral cell is an empty cell, then the attributes of U are the same as the attributes of the cell, and the attributes of the cell have been determined as the leaf cells of the lower grid in the bottom-up determination process. If the quadrilateral cell is a non-empty cell, since the center point of the cell has been determined in the bottom-up process, the method in step (3.1) can be used to determine the attributes of U based on the attributes of the center point, that is, draw a line segment connecting the center point and U, and calculate the number of intersections of the line segment with the polygon edges inside the quadrilateral. If it is an even number, then the attributes of U and the center point are the same, otherwise they are different. If the quadrilateral cell is a subdivided cell, then we iterate downward from it to find the central subcell of the current cell until we reach the leaf cell. The leaf cell is still a quadrilateral cell containing the midpoint. Therefore, if the leaf cell is an empty cell, then U has the same attributes as it. If the leaf unit is not an empty unit, then make a line segment connecting the center point of the leaf unit and U, and calculate the number of intersections between the line segment and the polygon edge inside the quadrilateral. If it is an even number, then the attributes of U and the center point are the same, otherwise they are different. Then use the attributes of U as the unit Attributes.

[0077] Finally, the attributes of empty cells in the current grid layer that are not adjacent to the subdivided cells are determined in parallel. Specifically, a thread is assigned to each cell in the current grid layer. If the cell corresponding to the thread is an empty cell that is not adjacent to the subdivided cell and has no attributes set, the attributes of the adjacent empty cell that is adjacent to the subdivided cell are used as the attributes of the thread.

[0078] (4) Determine whether the measured point is located within the spherical polygon

[0079] Based on the AHHG structure, which predetermines leaf cell attributes, we can achieve rapid, parallel determination of a point's location within a spherical polygon. Specifically, we assign a thread to each point being measured. For each point, we search down through the grid, starting from the first layer, until we find a cell containing the point. We then determine the point's attributes based on that leaf cell. Each thread performs the following operations:

[0080] (4.1) Find the leaf cell containing the measured point

[0081] The first-level grid cells are traversed sequentially until a cell containing the measured point is found and used as the current cell. If the current cell is not a leaf cell, the process starts from the current cell and searches downwards for a leaf cell containing the measured point. The process is to traverse the center subcell of the current cell and all its boundary subcells until a subcell containing the measured point is found. If the subcell is not a leaf cell, it is used as the current cell and the process is repeated.

[0082] (4.2) Determine the attributes of the measured point based on the leaf unit attributes

[0083] If the leaf cell is an empty cell, the attributes of the measured point are the same as those of the leaf cell; if the leaf cell is a non-empty cell, take the center point of the leaf cell, make a connecting line segment between the center point and the measured point, and calculate the number of intersections between the line segment and the polygon edge in the current cell; if it is an even number, the attributes of the measured point are the same as those of the center point of the leaf cell; if it is an odd number, the attributes of the measured point are different from those of the center point of the leaf cell.

[0084] An embodiment of the present application also provides a computer-readable storage medium, on which a computer program is stored. When the computer program is executed by a processor, the various processes of the embodiment of the above-mentioned GPU-based method for determining whether a point is located within a polygon are implemented, and the same technical effect can be achieved. To avoid repetition, it will not be repeated here.

[0085] The various embodiments in this specification are described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the various embodiments can be referenced to each other.

[0086] Those skilled in the art will appreciate that the embodiments of the present application can be provided as methods, devices, or computer storage media. Therefore, the present application can adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Moreover, the present application can adopt the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to magnetic disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0087] The embodiments of the present application are described with reference to the flowcharts and / or block diagrams of the methods, terminal devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing terminal device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing terminal device generate instructions for implementing the steps in the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0088] These computer program instructions may also be stored in a computer readable memory that can direct a computer or other programmable data processing terminal device to operate in a specific manner, so that the instructions stored in the computer readable memory produce a manufactured product including an instruction device, which implements the process Figure 1 a process or multiple processes and / or boxes Figure 1 The function specified in one or more boxes.

[0089] These computer program instructions can also be loaded onto a computer or other programmable data processing terminal device so that a series of operating steps are executed on the computer or other programmable terminal device to produce a computer-implemented process, thereby providing instructions for executing on the computer or other programmable terminal device to implement the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0090] Although preferred embodiments of the present invention have been described, those skilled in the art may make additional changes and modifications to these embodiments once they become aware of the basic inventive concepts. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments and all changes and modifications that fall within the scope of the embodiments of the present invention.

[0091] This document uses specific examples to illustrate the principles and implementation methods of this application. The description of the above embodiments is only used to help understand the method and core ideas of this application. At the same time, for those skilled in the art, based on the ideas of this application, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting this application.

Claims

1. A GPU-based method for determining whether a point is located within a spherical polygon, characterized in that: The following steps are involved: (1) Pre-allocate the adaptive hexagonal hierarchical grid AHHG structure and allocate a vertex structure array on the GPU side V , an array of grid edge structures B , a cell structure array C , an array of polygon edge indices R， On the CPU side, a grid structure array G that records the storage location of geometric elements is pre-allocated; on the CPU side, an AHHG first-layer grid covering the entire sphere is created in a serial manner, and the first-layer grid data is copied to the V, B, and C arrays on the GPU side in hierarchical order; (2) Starting from the first layer of grid, the AHHG structure is created layer by layer by iteratively and parallelly subdividing the non-empty cells containing polygon edges in the current layer of grid; when the cutoff subdivision condition is met, the subdivision is stopped, and then the edges of the polygons are distributed in parallel to the bottom layer of grid cells, and finally the creation of the AHHG structure is completed; (3) Parallel determination of the attributes of all leaf cells in the bottom grid of the AHHG structure, and then parallel determination of the attributes of leaf cells in each grid layer starting from the grid layer above the bottom grid to the first grid layer; (4) Based on the AHHG structure, it is determined in parallel whether the measured point is located within the spherical polygon; The creation of the AHHG structure in (2) includes: assigning a thread to each polygon edge, tracking the cells that the edge passes through and marking non-empty cells; Calculate the average number of polygon edges contained in each non-empty cell in the current layer. If the number is less than a given threshold or the current number of layers is equal to the maximum number of layers set by the user, stop subdividing. Otherwise, create the center sub-cell and the boundary sub-cell in parallel. Assign threads to each boundary sub-unit and set the sub-unit adjacency relationship through the partner edge relationship; Distribute polygon edges to the bottom-level grid cells in parallel: first obtain (edge ​​pointer, cell index) pairs in edge pointer-major order, then convert them into (cell index, edge pointer) pairs in cell index-major order and store them in the R array.

2. The method according to claim 1, wherein In step (1), all vertices of the layers are stored in the array in order of layers. V , that is, first the vertices of the first layer of mesh, then the vertices of the second layer of mesh, and so on; similarly, the mesh edges of all levels are stored in the array in order of level B , all levels of cells are stored in the array in order of level C ; Array V , B and C The elements in the grid are searched using an index or a level-relative index; the index refers to its sequence number in the entire array element, while the level-relative index refers to its sequence number in the array element of this layer; the vertex structure stores the grid unit vertex information, that is, the Cartesian coordinates of the vertex; the grid edge structure stores the information of the grid unit directed edges; the cell structure stores the grid unit information; the polygon edge index array R The index of the polygon edge contained in the storage unit in the polygon edge array is called the edge pointer; the CPU side storage space includes a grid structure array G , which records the starting storage position and number of geometric elements of all levels, including vertices, mesh edges and cells, in the corresponding arrays.

3. The method according to claim 1, wherein The characteristic of the first-layer grid of AHHG created in a serial manner on the CPU side in step (1) is that it covers the entire sphere.

4. The method according to claim 3, wherein The first-level grid of the AHHG includes 80 vertices, 240 directed edges and 42 units, including 12 pentagonal units and 30 hexagonal units.

5. The method according to claim 1, wherein The method of step (2) creating the AHHG structure layer by layer by iteratively subdividing the cells containing polygon edges in the current layer grid starting from the first layer grid is as follows: marking all non-empty cells in the current layer in parallel; determining whether to end the subdivision, if so, allocating the polygon edges to the current layer grid in parallel; if not, creating all the center sub-cells of the current layer grid in parallel; creating all the boundary sub-cells of the current layer grid in parallel; and setting the adjacent relationship between all the newly generated sub-cells in parallel; Use the new layer of mesh generated by subdivision as the current layer of mesh.

6. The method according to claim 5, wherein The method for parallel marking all non-empty cells in the current layer is: a thread is assigned to each polygon edge, each thread first locates the current layer cell where the starting point of the polygon edge is located, and then tracks the polygon edge, that is, first finds the cell edge that the polygon edge passes through, and finds the next cell that the polygon edge enters based on the adjacent relationship between the cell edges, until reaching the cell where the end point of the polygon edge is located; during this process, marks are set in the structures of all traversed cells.

7. The method according to claim 5, wherein The method for creating all the central sub-units of the current layer grid in parallel is: calculating the number of new geometric elements from the central sub-units in parallel; calculating the writing positions of the new geometric elements in parallel; and writing the new geometric elements in parallel.

8. The method according to claim 7, wherein The method of parallelizing the number of new geometric elements from the central subcell is to use two buffer arrays, buf 0 and buf 1, to record the number of new vertices and new units of each central subunit respectively; Assign a thread to each current layer unit; for a thread i , if the corresponding unit is a marked unit, then let buf 0[ i ] is the number of vertices of the unit, buf 1[ i ] is 1; otherwise , buf 0[ i ]and buf 1[ i ] are all 0.

9. The method according to claim 7, wherein The method of parallel calculation of the writing position of the new geometric elements is: parallel calculation buf 0 and buf 1's pre-order sum and store the result back buf 0 and buf 1.

10. The method according to claim 9, wherein The method of writing new geometric elements in parallel is: assign a thread to each current layer unit; if the current unit is a marked unit, generate new geometric element data and write it according to buf 0 and buf The value in 1 is written to V , B , C The corresponding position in the array, the specific process is: set the thread i The corresponding unit is c ,like c is the marked unit, k The grid edges are composed of c The total number of central subunits k A new vertex, k New mesh edges and 1 new unit; let the vertices of the next layer of mesh in the current layer be V The starting storage location in is vstart , then from vstart + buf 0[ i ]Starting position in turn to V Write in k New vertex data; Let the mesh edge of the next layer of mesh in the current layer be B The starting storage location in is bstart , then from bstart + buf 0[ i ]Starting position in turn to B Write in k new mesh edges; Let the cells of the next layer of grid in the current layer be C The starting storage location in is cstart , then to C middle cstart + buf 1[ i ] position to write the new unit.

11. The method according to claim 5, wherein The method for creating all boundary sub-units of the current layer grid in parallel is: calculating the number of new geometric elements from the boundary sub-units in parallel; calculating the writing positions of the new geometric elements in parallel; and writing the new geometric elements in parallel.

12. The method according to claim 11, wherein The method of parallelizing the number of new geometric elements from the boundary subcells is to use two buffer arrays, buf 0 and buf 1, to record the number of new mesh edges and new cells from each boundary sub-cell respectively; assign one thread to each edge of the newly generated central sub-cell, and assign one thread to each edge of the newly generated central sub-cell. i , if the corresponding edge b Belongs to the center sub-unit of the quadrilateral and is the longest side in the unit, then this side has no corresponding boundary sub-unit, let buf 0[ i ]=0, buf 1[ i ]=0; Otherwise, first look for b The upper unit of the central sub-unit is called the parent unit, and then find the parent unit with the same b The parallel edge is called the parent edge, and the partner edge of the parent edge is called the parent partner edge. Finally, the cell to which the parent partner edge belongs is called the adjacent parent cell. If the adjacent parent cell does not exist or the adjacent parent cell has no central child cell, the number of edges of the boundary child cell is 4. Let buf 0[ i ]=4, buf 1[ i ]=1; if the adjacent parent cell exists and has a central child cell, compare the cell indexes of the parent cell and the adjacent parent cell; if the parent cell index is larger, the boundary child cell is a hexagonal cell, let buf 0[ i ]=6, buf 1[ i ]=1.

13. The method according to claim 11, wherein The method of parallel calculation of the writing position of the new geometric elements is: parallel calculation buf 0 and buf 1's pre-order sum, and store the results back separately buf 0 and buf 1.

14. The method according to claim 11, wherein The method of writing new geometric elements in parallel is: first, bstart =The last grid edge of the center sub-unit of the current layer grid is B Index +1 in, cstart =The last cell of the center subcell of the current layer grid is in C Index in +1; Then assign a thread to each edge of the newly generated central sub-unit, and assign a thread to each edge of the newly generated central sub-unit. i , if the corresponding edge b If it belongs to a quadrilateral center sub-unit and is the longest side in the unit, then this side has no corresponding boundary sub-unit and no operation is performed; Otherwise, search for b The parent unit of the central child unit, b On the father's side, b The father partner Bian He b If the adjacent parent cell does not exist or the adjacent parent cell has no center child cell, a quadrilateral boundary child cell is generated and the four new mesh edges are written into the B Zhongcong bstart+buf 0[ i ] At the beginning of the C middle cstart+buf 1[ i ]; if the adjacent parent cell exists and has a center child cell, compare the cell index of the parent cell and the adjacent parent cell; if the parent cell index is larger, generate a hexagonal boundary child cell and write the 6 new grid edges into B Zhongcong bstart+buf 0[ i ] At the beginning of the C middle cstart+buf 1[ i ] position.

15. The method according to claim 6, wherein The method of setting the adjacent relationship between all newly generated sub-units in parallel is: assigning a thread to each boundary sub-unit, a , process each of its mesh edges in turn; for the unit a Each mesh edge u , first find the unit a The parent unit a’ , then find a’ The corresponding boundary subunit is a Mesh edge u 1', then find the unit a 'Middle and grid edge u 1' connected mesh edges u 2 ’ , and then find the grid edge u 2 ’ Corresponding boundary subunit b , from the unit b Find the edges of the mesh u Mesh edges with the same vertices w , the mesh edges u and mesh edges w Set each other as partner sides.

16. The method according to claim 1, wherein The method to obtain (edge ​​pointer, unit) pairs in parallel with edge pointer as the main order is: set up four temporary integer arrays, buf 0, buf 1, buf 2 and buf 3. Used as buffer storage and initialized to all zeros; Parallel calculation of the number of cells each polygon edge passes through and writes it to buf 0 in parallel computing buf 0's pre-order sum, and store it back buf 0; According to buf The value of 0 is buf 1 writes the edge pointer in the (edge ​​pointer, cell index) pair to buf Write the cell index in the (edge ​​pointer, cell index) pair in 2.

17. The method according to claim 1, wherein The method of converting (edge ​​pointer, cell index) pairs into (cell index, edge pointer) pairs in parallel with cell index as the main order is: first, according to buf A value of 1 calculates the number of polygons contained in each cell in parallel and writes them into an array buf 3; then calculate buf 3's precedence and write to buf0 ,according to buf 1, buf 2 and buf The value of 3 writes the edge pointer to the R In the array; according to buf 3 and buf0 The value of , and set the edge pointer of each unit in parallel R The starting position and quantity in .

18. The method according to claim 1, wherein The step (3) of determining the attributes of all leaf cells of the bottom grid of the AHHG structure in parallel is as follows: determining the attributes of all non-empty cells of the bottom grid in parallel; determining the attributes of all empty cells adjacent to non-empty cells in the bottom grid, which are called direct empty cells, in parallel; In parallel, the properties of all empty cells in the bottom grid that are not adjacent to non-empty cells are called non-direct empty cells.

19. The method according to claim 18, wherein The method for parallel determination of the non-empty cell attributes of the bottom-level grid is as follows: a thread is assigned to each bottom-level grid cell; if the cell corresponding to the thread is a non-empty cell, its center point attribute is set, that is, the edge of the polygon closest to the center point in the cell is found, and it is determined whether the center point is located to the left of the plane formed by the edge and the center of the sphere; if so, the center point is inside the polygon, otherwise it is outside the polygon.

20. The method of claim 18, wherein: The method for parallelly determining the attributes of all empty cells adjacent to non-empty cells in the bottom-level grid as direct empty cells is as follows: assign a thread to each bottom-level cell; if the cell corresponding to the thread is a direct empty cell with no attributes set, take the edge adjacent to the non-empty cell, take the midpoint on the edge, connect it with the center point of the non-empty cell, and calculate the number of intersections between the line and the polygonal edge in the non-empty cell; if the number of intersections is an even number, then the attributes of the empty cell corresponding to the thread are the same as those of the center point, otherwise they are different.

21. The method of claim 18, wherein: The method for parallelly determining the attributes of all empty cells in the bottom-level grid that are not adjacent to non-empty cells as non-direct empty cells is as follows: a thread is assigned to each bottom-level cell; if the cell corresponding to the thread is a non-direct empty cell with no attributes set, the attributes of the direct empty cells adjacent to the thread are used as the attributes of the non-direct empty cell.

22. The method of claim 1, wherein The method for determining the leaf unit attributes of each layer of grid in parallel in step (3) is: if the previous layer of grid of the current layer does not exist, then end; Otherwise, the previous grid of the current grid is used as the current grid; the properties of all empty cells of the current grid that are adjacent to the subdivided cells are determined in parallel; the properties of all empty cells of the current grid that are not adjacent to the subdivided cells are determined in parallel.

23. The method according to claim 22, wherein The method of determining the attributes of all empty cells adjacent to the subdivided cells in the current layer grid in parallel is as follows: a thread is assigned to each current layer cell; if the cell corresponding to the thread is an empty cell adjacent to the subdivided cell with no attribute set, c , which is subdivided into units c 'Adjacent, then look for c 'Zhong and c Adjacent mesh edges b ,Pick b midpoint U， Then find b Corresponding boundary subunit k； If the unit k If it is an empty cell, U The properties of the unit are the same as the properties of the unit; if the unit k If it is a non-empty unit, it is a connecting unit. k The center point and U The line segment is calculated by the intersection of the line segment and the unit k If the number of intersections of the inner polygon edges is even, then U Same as the properties of the center point, otherwise different; if the unit k If it is a subdivided unit, it will iterate downward from it to find the central subunit of the current unit until it reaches the leaf unit; if the leaf unit is an empty unit, then U With the same properties, if the leaf unit is not an empty unit, then connect the center point of the leaf unit and U , calculate the number of intersections between the line segment and the polygon edge in the leaf unit. If it is an even number, then U Same as the attributes of the center point, otherwise different; finally, U Attributes as units c Attributes.

24. The method of claim 22, wherein: The method for parallelly determining the attributes of all empty cells in the current layer grid that are not adjacent to the subdivided cells is: assigning a thread to each current layer cell; if the cell corresponding to the thread is an empty cell that is not adjacent to the subdivided cell and has no attributes set, then the attributes of the cells adjacent to the subdivided cell are used as the attributes of the cell corresponding to the thread.

25. The method of claim 1, wherein The method for determining in parallel whether the measured point is located within the spherical polygon based on the AHHG structure in step (4) is as follows: assigning a thread to each measured point; for each measured point, searching for a cell containing the measured point layer by layer starting from the first layer of the grid downward until the cell containing the measured point is a leaf cell; and determining the attributes of the measured point based on the leaf cell.

26. The method of claim 25, wherein: The method of searching for the unit containing the measured point layer by layer starting from the first layer of grid until the unit containing the measured point is a leaf unit is: traversing the first layer of grid units in sequence until the unit containing the measured point is found and taking it as the current unit; if the current unit is not a leaf unit, searching for the leaf unit containing the measured point layer by layer starting from the current unit.

27. The method according to claim 26, wherein The method of searching for the leaf unit containing the measured point layer by layer starting from the current unit is: traversing the central sub-unit of the current unit and all its boundary sub-units in turn until the sub-unit containing the measured point is found. If the unit is not a leaf unit, it is used as the current unit and this step is repeated.

28. The method of claim 25, wherein: The method for determining the attribute of the measured point based on the leaf cell is as follows: if the leaf cell is an empty cell, the attribute of the measured point is the same as the attribute of the leaf cell; if the leaf cell is not an empty cell, the center point of the leaf cell is taken, a line segment is made connecting the center point and the measured point, and the number of intersections of the line segment and the polygon edge in the leaf cell is calculated; If it is an even number, the attributes of the measured point are the same as those of the center point of the leaf unit; if it is an odd number, the attributes of the measured point are different from those of the center point of the leaf unit.

29. A computer storage medium, characterized in that The computer storage medium stores computer program instructions, which, when executed by a processor, implement a GPU-based method for determining whether a point is located within a spherical polygon as described in any one of claims 1 to 28.

Images