Geometric model quadrilateral mesh generation method and system based on GPU parallelism

Abstract

The invention discloses a geometric model quadrilateral mesh generation method and system based on GPU parallelism, and belongs to the technical field of computer graphics. The method comprises the following steps: for a geometric model formed by a plurality of surfaces, obtaining parameterized coordinates of surface patches in parallel by a GPU thread, and performing spatial transformation on the parameterized coordinates by establishing a bounding box of the surface patches to obtain a surface patch splicing result; generating an original triangular mesh model based on a patch splicing result, and processing triangular units in the original triangular mesh model in parallel by the GPU thread to obtain a triangular mesh after cavity filtering and network geometric information of the triangular mesh; the GPU thread solves the network topology structure of the triangular mesh in parallel; the GPU thread performs parallel solution of the frame field based on the network geometric information and the network topology information of the triangular mesh; and generating a quadrilateral grid structure of the geometric model in combination with the frame field. According to the method, the grid structure generation efficiency is improved by improving the calculation performance of the curved surface frame field.

Description

Technical Field

[0001] The present invention relates to the field of computer graphics technology, and in particular to a method and system for generating quadrilateral meshes of geometric models based on GPU parallelism. Background Art

[0002] Frame fields are widely used in computer graphics. Frame field calculation is a crucial step in generating high-quality quadrilateral meshes for curved surfaces. It provides clear direction for mesh generation and, through its singularity control mechanism, effectively regulates key nodes in mesh generation, thereby aiding in the generation of high-quality quadrilateral meshes and significantly improving the accuracy and efficiency of quadrilateral mesh generation for complex surfaces.

[0003] Mainstream methods for calculating surface frame fields are divided into two stages: initial surface triangulation and frame field solution on the triangular mesh. The mainstream approach for initial surface triangulation is to generate a triangular mesh in a 2D parameterized space using the Delaunay algorithm, which is then projected onto a 3D surface. Solving the frame field is performed on the edges of the triangular mesh. This requires first calculating the topological relationships of the triangular mesh, then generating a system of linear equations based on these topological relationships to obtain the frame field.

[0004] 1. Surface triangulation.

[0005] Surface triangulation, the most mature technique in current surface polygon meshing research, is widely used in geographic modeling, 3D reconstruction, and graphics rendering. In 2008, Rong et al. proposed a 2D Delaunay triangulation method based on GPU-CPU collaboration. The method's main steps include parallel construction of discrete Voronoi diagrams on the GPU, followed by refinement and correction of the Delaunay structure on the CPU. This method encompasses ten steps: texture mapping, a skip flooding algorithm, initial mesh construction on the GPU, convex hull correction on the CPU, missing point insertion, and edge flipping optimization. It achieves a 53% speedup over traditional CPU algorithms for evenly distributed point sets of millions. Subsequently, Qi et al. extended this method to a full GPU implementation in 2011, supporting parallel processing of constraint insertion and edge flipping operations. This method achieved speedups of up to 31x on synthetic data and approximately 4x on real GIS data, significantly enhancing its applicability for processing large point sets. The GPU-DT algorithm builds on this work, leveraging CUDA parallelism to significantly improve the speed of triangulating large point sets.

[0006] While these algorithms have demonstrated significant advantages across a wide range of applications, they still face numerous limitations. A typical issue is the use of GPU Delaunay triangulation for each facet when processing multi-faceted surface models. This unit-based triangulation approach not only fails to fully utilize the parallel computing capabilities of the GPU, but also incurs high I / O overhead due to frequent data transfers, significantly impacting overall triangulation efficiency.

[0007] 2. Frame field The frame field is a four-fold symmetric direction field defined based on the quotient group space , with four-fold symmetry. The development of frame fields has gone through several important stages. In 2000, Hertzman et al. first proposed the concept of frame fields (cross fields) for texture rendering, but this technology did not receive widespread attention until Palacios et al. introduced the rotationally symmetric field (rotor field) in 2007. On this basis, Panozzo et al. proposed a frame field generation method based on telescopic orthogonal fields, which significantly improved the universality of the frame field. Jiang et al. reduced the variable dimension by introducing constraints and combined it with a customized Riemannian metric algorithm to further expand the application scope of surface frame fields. In 2017, Beaufort et al. innovatively introduced the Ginzburg-Landau theory into the frame field solution, providing an efficient numerical solution for this problem. The latest progress came from Reberol et al. in 2021. They constructed a quadrilateral mesh generation framework for CAD models based on the frame field, providing an important tool for the field of industrial design.

[0008] Computing frame fields on curved triangular meshes involves calculating mesh topological relationships and solving linear equations. The complexity of this problem increases dramatically with increasing mesh size. Generating a large number of topological relationships between mesh cells and nodes is fundamental to frame field calculations, and the storage and computational overhead of this process increases quadratically with increasing mesh resolution. Constructing and solving a system of linear equations based on the Ginzburg-Landau theory is a key step in generating the frame field. As the mesh size increases, the coefficient matrix dimension increases significantly, leading to a significant increase in the storage resources used and the solution time for the system. Given these issues, computing high-resolution frame fields requires significant storage and computational resources, making it difficult to implement. Eigen is a C++ linear algebra library widely used in numerical computing, and most frame field solutions on CPUs are implemented using Eigen.

[0009] In summary, existing algorithms perform serial calculations on the CPU, and face obvious bottlenecks in computational efficiency when faced with the demand for high-precision frame field calculations on complex surfaces: for the initial triangulation process of the surface, the mainstream Delaunay triangulation method consumes a large amount of time for computationally intensive operations such as topological queries and hole detection as the cell size increases; solving the frame field on the triangular mesh requires constructing a global sparse linear system of equations, which is inefficient in matrix assembly and difficult to solve. Direct methods are limited by memory usage and cannot handle large-scale problems, and iterative methods are affected by convergence speed and difficult to guarantee solution performance. Summary of the Invention

[0010] The present invention discloses a method and system for generating quadrilateral meshes of geometric models based on GPU parallelism, which increases the efficiency of generating mesh structures of geometric models by improving the computing performance of curved surface frame fields.

[0011] To achieve the above objectives, the technical solution of the present invention includes the following contents.

[0012] A method for generating quadrilateral meshes of geometric models based on GPU parallelism, the method comprising: For a model composed of multiple faces, each GPU thread obtains the parametric coordinates of the face in the UV coordinate system in parallel, and performs spatial transformation on the parametric coordinates of the face by establishing the face's bounding box to obtain the face splicing result; Generate an original triangular mesh model based on the patch stitching result, and use each GPU thread to process each triangular unit in the original triangular mesh model in parallel to obtain a cavity-filtered triangular mesh and network geometry information of the triangular mesh, wherein the network geometry information includes vertices and triangles; Each GPU thread solves the network topology structure of the triangular mesh in parallel; wherein the network topology structure includes: the topological composition of the edge, the adjacent triangles of the edge, and the adjacent edges of the vertex; Each GPU thread performs parallel solution of the frame field based on the network geometry information and network topology information of the triangular mesh; The frame field is combined to generate a quadrilateral grid structure of a geometric model.

[0013] Furthermore, each GPU thread obtains the parametric coordinates of the patch in the UV coordinate system in parallel, and performs spatial transformation on the parametric coordinates of the patch by establishing a bounding box of the patch to obtain the patch stitching result, including: Each GPU thread obtains the parameterized coordinates of all points in the patch in the UV coordinate system in parallel , Indicates the parameterized coordinate value of a point in the patch in the horizontal direction in the UV coordinate system. Indicates the parametric coordinate value of a point in the patch in the vertical direction in the UV coordinate system; For each face of the model, a bounding box is used to enclose the face; wherein the bounding box is the minimum rectangular frame that contains the parameterized plane of the face, and the extreme value coordinates of the bounding box include: horizontal extreme value U_min, horizontal extreme value U_max, vertical extreme value V_min and vertical extreme value Vmax; Each GPU thread performs vertical offset preprocessing based on the vertical extreme value Vmax and the patch processing order index in parallel on a patch basis to obtain the vertical reference value of the patch. Each GPU thread obtains the position of the corresponding patch in the patch queue according to the patch processing order index, and obtains the total offset of the preceding patches of the patch based on the vertical reference of all patches before the patch. Each GPU thread is based on the total offset of the previous patch of the patch, the parameterized coordinate value Perform vertical spatial transformation on all points in the patch and the vertical extreme value V_min to obtain the transformed parameterized coordinate values ​​of the corresponding points ; Get the parameterized coordinate values ​​of all points in the patch after spatial transformation .

[0014] Furthermore, the method for generating the original triangular mesh model based on the patch stitching result includes: GPU threads executing the Delaunay decomposition algorithm in parallel.

[0015] Furthermore, each GPU thread is used to process each triangle unit in the original triangle mesh model in parallel to obtain a cavity-filtered triangle mesh and network geometry information of the triangle mesh, including: Each GPU thread calculates the centroid coordinates of different triangle units in parallel and constructs a centroid-triangle unit index mapping table; Construct a first ray along the extended diagonal direction of the bounding box, and select an external point as the end point of the first ray; Generate a second ray starting from the centroid of the triangular element and ending at the external point; Establishing a boundary edge quadtree space partitioning structure, and combining the boundary edge quadtree space partitioning structure with calculating the intersection detection between the second ray and the edge of the triangle unit to obtain an intersection count of the triangle unit; Perform Boolean judgment on the intersection count of each triangle unit to obtain the reserved flag bitmap of the triangle unit; Compressing the original triangular mesh model based on the retention flag bitmap of the triangular unit, and after searching the corresponding triangle index using the centroid-triangle index mapping table, deleting the triangles with the mark bit of 0 in the retention flag bitmap to obtain a triangular mesh after cavity filtering; Extract the network geometry information of the triangular mesh ;in, represents a vertex, Represents a triangle.

[0016] Furthermore, each GPU thread solves the network topology structure of the triangular mesh in parallel, including: For each vertex in the triangle mesh , traverse the triangle in parallel through GPU threads To generate a P_Tris array, the P_Tris array is used to store the adjacent triangle index set of the vertex; Determine vertices based on the P_Tris array Is there an adjacent edge at this vertex? If there are adjacent edges, perform parallel prefix sum calculations based on the P_Tris array to generate the P_Edges_Num_Prefix_Sum array, which is used to store the prefix sums of the adjacent edges of the vertex; Allocate storage space for the P_Edges array based on the prefix sum of the adjacent edges, and write the vertices through multi-threaded writing. The adjacent edge index of the vertex is written into the P_Edges array, and the P_Edges array is used to store the adjacent edge index set of the vertex; Remove duplicates from the adjacent edges in the P_Edges array and write the deduplicated adjacent edge index set into the Point_Start_Edge_Num array; Perform parallel prefix sum calculation based on the Point_Start_Edge_Num array to generate the P_Start_Edge_Num_Prefix_Sum array, where the P_Start_Edge_Num_Prefix_Sum array is used to store the prefix sums of the endpoints of the adjacent edges; Based on the prefix sum of the adjacent edge endpoints, the adjacent edge endpoints are written into the E_Fir_Point array and the E_Sec_Point array respectively; wherein, for any adjacent edge, the index value of the adjacent edge endpoint in the E_Fir_Point array is less than the index value of the adjacent edge endpoint in the E_Sec_Point array; Perform parallel edge traversal on the triangle mesh and perform a parallel edge traversal on each edge. , retrieve the vertex The adjacent triangle set P_Tris[ ];in, Represents an edge Another vertex of Based on the adjacent triangle set P_Tris[ ] Perform triangle containment detection, and based on the detection results, store the number of adjacent triangles of the edge and the index of the adjacent triangles of the edge in the E_Tris_Num array and the E_Tris array respectively; The topological composition of the edge, the adjacent triangles of the edge, and the adjacent edges of the vertex are obtained through the E_Fir_Point array and E_Sec_Point array, the E_Tris_Num_Prefix_Sum array and E_Tris array, and the P_Edges_Num_Prefix_Sum array and P_Edges array respectively.

[0017] Furthermore, the vertex is realized by the following steps and vertices The network topology query operation between: According to the vertex index value, the vertex and vertices Perform sequence normalization; Get the vertex through the P_Start_Edge_Num_Prefix_Sum array As the storage starting position of the edge head node, start; According to Point_Start_Edge_Num[ ] Calculate the end position end = start + out-degree, lock the candidate edge set E_Sec_Point[start : end]; Starting from the starting position start, access the vertices one by one in memory order Associated edge, read the secondary endpoint index current_ = E_Sec_Point[idx], if current_ == : Immediately returns the current global index idx and ends the search when traversal reaches position end; if no match is found, returns -1.

[0018] Furthermore, each GPU thread performs parallel solution of the frame field based on the network geometry and topology information of the triangular mesh, including: Construct the optimization problem: ; Among them, the first component , the second component , represents the frame field, represents the surface, Indicates adjustable parameters; Using the heat diffusion system to derive the optimization problem, we obtain the equation: ;in, represents the number of iterations, Represents a constant related to the diagonal length and the shortest side length of the surface bounding box. is a diagonal matrix, the elements in the diagonal matrix , 、 Represents edges The area of ​​an adjacent triangle of Represents the stiffness matrix, condition matrix is A single-column matrix, Indicates the The first component in the iteration The Quantity, Indicates the The first component in the iteration The Components, initial condition matrix The construction follows the principle of boundary distinction; Calculate the diagonal matrix based on the mesh topology information and mesh geometry information and the stiffness matrix The value of each element in; The diagonal matrix and the stiffness matrix Substitute the values ​​of each element in into the equation to solve the conditional frame field.

[0019] Furthermore, based on the grid topology information and grid geometry information, the diagonal matrix is ​​calculated and the stiffness matrix The values ​​of each element in include: Based on the adjacent triangles of the edge, parallel threads are synchronized to obtain the diagonal matrix The value of each element in; Traverse each edge and find the adjacent triangle of the current edge based on the E_Tris array and adjacent triangles ; In adjacent triangles and adjacent triangles Find the four adjacent edges of the current edge and calculate the angle between the four adjacent edges and the current edge , angle , angle and angle Then, based on the angle , angle , angle and angle Calculate the weight of the current edge; Calculate adjacent triangles and adjacent triangles The local coordinate system of the current edge is calculated based on the local coordinate system. ; Based on the angle of the current edge and weights, calculate the stiffness matrix The value of the corresponding element in .

[0020] Furthermore, the diagonal matrix and the stiffness matrix Substitute the values ​​of each element in into the equation to solve the conditional frame field, including: Step 1: Calculate the matrix ; Step 2: Based on the diagonal length of the bounding box and the shortest side length of the grid , calculate the initial angle and the final angle XX ; Step 3: Set the total number of loops , current number of cycles , and execute Step 4-Step 8, each cycle ; Step 4: Based on the initial angle , final angle , Current cycle number and the total number of cycles , calculate the current angle ; Step 5: Based on the current angle , diagonal matrix and the stiffness matrix Parallel matrix updates , and let the number of iterations ; Step 6: Based on the current angle and the conditional matrix , calculate the matrix ; Step 7: Solve using the parallel preconditioned conjugate gradient method ; Step 8: Parallel Update , ; Step 9: If , then go to Step 7; where c represents a set threshold; Step 10: If , then calculate the edge frame field.

[0021] A GPU-based quadrilateral mesh generation system for geometric models, comprising: The UV stitching module is used to obtain the parametric coordinates of the facets in the UV coordinate system in parallel for models composed of multiple faces. It then performs spatial transformation on the parametric coordinates of the facets by establishing a bounding box for the facets to obtain the facet stitching result. A triangular mesh generation module is used to generate an original triangular mesh model based on the facet splicing result, and use each GPU thread to process each triangular unit in the original triangular mesh model in parallel to obtain a triangular mesh after cavity filtering and network geometry information of the triangular mesh, wherein the network geometry information includes vertices and triangles; A topology structure generation module is used for each GPU thread to solve the network topology structure of the triangular mesh in parallel; wherein the network topology structure includes: the topological composition of the edge, the adjacent triangles of the edge, and the adjacent edges of the vertex; The frame field generation module is used by each GPU thread to perform parallel solution of the frame field based on the network geometry information and network topology information of the triangular mesh; The quadrilateral mesh generation module is used to generate the quadrilateral mesh structure of the geometric model in combination with the frame field.

[0022] Compared with the prior art, the present invention has at least the following beneficial effects.

[0023] 1) The UV stitching algorithm of the present invention maps the parameterized coordinates of multiple surfaces into a unified texture space to construct a virtual stitching point cloud. While achieving an order of magnitude increase in the data scale of concurrent computing, it also merges data communication batches between the CPU and GPU through the coordinate mapping mechanism, significantly reducing data transmission overhead.

[0024] 2) This paper designs a parallel topology solver with an excellent space complexity of only O(N) (N is the number of points), significantly reducing memory usage compared to existing methods. When solving linear equations, it employs a parallel preconditioned conjugate gradient method. Compared to direct methods, this method avoids matrix decomposition, significantly reducing memory consumption and making it suitable for computing high-resolution frame fields.

[0025] 3) The process of solving the surface frame field in the present invention is all completed on the GPU, which can improve the computing performance by 1 to 2 orders of magnitude compared with existing methods. BRIEF DESCRIPTION OF THE DRAWINGS

[0026] Figure 1 Schematic diagram of the frame field of the present invention.

[0027] Figure 2 This is a flowchart of the GPU-based parallel fast calculation method for curved frame fields of the present invention.

[0028] Figure 3A Schematic diagram of UV occlusion of two surfaces.

[0029] Figure 3B Schematic diagram of incorrect mesh topology caused by Delaunay after UV occlusion.

[0030] Figure 4 The following are example images of different faces of the model. Different colors represent different faces.

[0031] Figure 5 Figure 1 shows an example of how the unmodified Delaunay algorithm cannot correctly handle cavities. (a) shows the original model, and (b) shows the model without the cavity.

[0032] Figure 6 Schematic diagram of UV stitching algorithm.

[0033] Figure 7A Unified parameterization after UV splicing.

[0034] Figure 7B This is the Delaunay effect after UV stitching.

[0035] Figure 8 Schematic diagram of the cavity detection algorithm principle.

[0036] Figure 9 This is an example diagram after cavity processing of the present invention.

[0037] Figure 10 is a matrix Schematic diagram of the calculation. The red edge is the current edge.

[0038] Figure 11 for Schematic diagram of the calculation. Among them, Z is the normal vector of triangle t1, and XYZ forms a right-handed direct coordinate system. For the edge Projection on the XY plane, for The angle of rotation about the X axis.

[0039] Figure 12 This is the Delaunay mesh obtained using GPU-DT and the method of the present invention.

[0040] Figure 13This is an example of the frame field calculation results obtained by the present invention in the experiment on the Mambo dataset. DETAILED DESCRIPTION

[0041] In the following description, the present invention is further described through specific embodiments so that those skilled in the art can have a more thorough understanding of the features and advantages of the present invention.

[0042] In this embodiment, the frame field of the present invention is defined on the edge of the triangular mesh. After the local coordinate system on the edge is determined, Figure 1 As shown, the local coordinate system can be used express.

[0043] like Figure 2 As shown, the present invention can be divided into two major stages: surface triangulation and frame field solution. In the surface triangulation stage, the present invention implements GPU-parallel constrained Delaunay triangulation and designs UV stitching and cavity processing operations to improve computational efficiency and correctly handle models with cavities. In the frame field solution stage, for high-resolution frame field GPU parallel solution, the present invention designs a grid topology data structure with a spatial complexity of O(N) and an access time complexity of O(1), and implements a rapidly convergent parallel preconditioned conjugate gradient method to solve linear equations.

[0044] 1. UV stitching.

[0045] For models composed of multiple faces, the two-dimensional parameter spaces of the multiple faces often block each other. If Delaunay meshing is performed directly in this parameterized space, the points between different faces will be connected to each other, resulting in incorrect mesh topology. For example, when the UV coordinates of the facets on the two-dimensional plane are as follows Figure 3A When the faces overlap, the traditional Delaunay triangulation algorithm will mistakenly connect the sampling points of different faces into cross-face triangles ( Figure 3B This cross-surface connection destroys the geometric continuity of the surface, resulting in wrinkles or tears when mapped back to the 3D model.

[0046] In order to allow the existing GPU Delaunay algorithm to perform Delaunay on multiple faces at the same time, this paper proposes a UV splicing strategy that can solve the UV occlusion problem. Figure 4For each face of the model, the present invention first calculates its parameterized bounding box, which is the smallest rectangular box that can completely accommodate the parameterized plane of the face, and it provides a quantifiable boundary range for subsequent operations. Then, through carefully designed translation operations, the parameterized bounding boxes of each face are spatially transformed. The purpose of these operations is to readjust the position and posture of the parameterization of each face in the parameter space so that there is no mutual occlusion between the bounding boxes. In this way, the parameterization results of multiple faces can be smoothly mapped to Figure 6 In the unified texture space shown in . Figure 7A This is a schematic diagram of unified parameterization after UV splicing.

[0047] During the entire process, since the limiting edges of each face have not changed, this ensures that no triangles crossing the face will be generated during the subsequent Delaunay triangulation process.

[0048] Input: Independent parameterized coordinates for each patch.

[0049] Output: Global UV layout with no overlap.

[0050] Step 1: Each GPU thread performs bounding box feature extraction in parallel on a facet basis, calculates the axis-aligned bounding box (AABB) for each facet, and records the bounding box's extreme coordinates U_min / U_max (horizontal extreme value) and V_min / V_max (vertical extreme value).

[0051] Step 2: Each GPU thread performs vertical offset preprocessing in parallel on a face-by-face basis, generating a vertical reference for each face: V_base = V_max + face_index × ε (ε = 1e-6 to prevent numerical conflicts, face_index is the processing order index).

[0052] Step 3: Spatial accumulation calculation: calculate the prefix sum array of the vertical reference in the processing order, use the classic prefix sum calculation algorithm on the GPU to obtain the total offset of the first i patches S[i] = Σ(V_base[0…i]), and build a cumulative height reference system.

[0053] Step 4: Each GPU thread performs coordinate space transformation in parallel on a patch basis, performing a vertical translation on the i-th patch: v' = v + (S[i-1] - V_min) (S[i-1] represents the total offset of the previous patch, and V_min is the original lower boundary of the current patch), keeping the horizontal U coordinate unchanged. v represents the coordinate of the point on the patch in the V direction in the UV parameter space.

[0054] Experimental verification shows that when processing tens of millions of points, the computational efficiency advantage of this method becomes increasingly pronounced as the number of faces increases. This is primarily because, without the Delaunay algorithm for UV stitching, GPU utilization is significantly reduced due to the small size of individual face data, which prevents full utilization of the GPU's computing cores. Furthermore, the time consumed by I / O operations resulting from multiple face-by-face data transfers to the GPU increases.

[0055] 2. Parallel processing of cavities.

[0056] 2.1 Cavity.

[0057] When the existing Delaunay triangulation algorithm on GPU processes surface models with cavities, due to the lack of effective cavity identification and boundary constraint mechanism, it is inevitable that the following problems will occur at the boundaries of adjacent faces: Figure 7B The extra triangles shown in the blue box of . These triangles can be regarded as triangles existing in the cavity. For special surfaces with internal cavities, this result is not correct, such as Figure 5 Based on parallel Delaunay triangulation, the present invention proposes a method for parallel removal of triangles in cavities, which can not only accurately generate internal edges but also eliminate cross-face triangles generated by UV stitching.

[0058] 2.2 Cavity parallel processing algorithm.

[0059] The current mainstream cavity detection algorithm mainly implements its core judgment logic based on the CPU, that is, based on the ray casting principle, such as Figure 8 As shown, any point inside the triangle to be tested is selected, and a ray pointing to the infinite point p is constructed. The cavity is determined by counting the parity of the number N of intersections between the ray and the model boundary (when N is an even number, it is determined to be a triangle inside the cavity, and when N is an odd number, it is a normal surface. However, the traditional implementation scheme requires all the geometric data generated by the GPU to be returned to the CPU for processing. Not only does it generate significant I / O overhead due to frequent CPU-GPU data transmission, but it also wastes hardware computing potential due to the inability to effectively utilize the GPU's stream processor architecture, which seriously restricts the algorithm execution efficiency. To address this bottleneck, the present invention proposes a cavity parallel processing framework based on the hardware characteristics of the GPU. By reconstructing the data structure and calculation paradigm of the ray intersection test, the full-process GPU-side calculation of cavity determination is realized, giving full play to the large-scale parallel computing advantages of modern graphics processors.

[0060] The specific algorithm steps are as follows: Input: original triangle mesh model, boundary edge set.

[0061] Output: Clean mesh after cavity filtering.

[0062] Step 1: Each GPU thread generates the barycenter coordinates in units of triangle units, synchronously calculates the centroid coordinates for each triangle unit, and constructs a centroid-triangle index mapping table.

[0063] Step 2: Dynamically construct the far end point, calculate the mesh axis-aligned bounding box (AABB), and select the external point along the diagonal extension direction of the bounding box as the ray end point.

[0064] Step 3: Accelerate the structure-assisted intersection test, establish a boundary edge quadtree space partitioning structure, and perform ray-edge intersection detection in parallel: (1) Among them, Q is the quadtree query result, is the intersection indicator function, Indicates the starting point , the end point is rays.

[0065] Step 4: Parallel determination of odd and even marks, and implementation of Boolean determination based on intersection count: (2) Generates a triangle preservation flag bitmap.

[0066] Step 5: Mesh reconstruction and parallel reduction operation. Based on the flag bitmap compression storage structure, the corresponding triangle index is found using the centroid-triangle index mapping table, and the triangles with the flag bit set to 0 are deleted to output a topologically coherent cavity filter mesh. Figure 9 (b) in the equation is correct. Figure 9 (a) is the effect diagram obtained by executing the cavity detection algorithm of the present invention on the original model.

[0067] 3. Parallel solution of topological structures.

[0068] When processing mesh data, complete topological structure information is crucial. Because the definition of the frame field relies on the geometry and connectivity of the mesh edges, key topological data such as the topological composition of the edges, adjacent triangle information, and vertex adjacent edges must be obtained. However, the mesh data structure generated by the commonly used Delaunay triangulation algorithm is relatively basic, typically only recording basic information such as triangle vertex indices. It neither explicitly stores edge adjacency relationships nor effectively identifies boundary edges, making it difficult to directly provide the topological relationship support required for frame field construction.

[0069] Key-value map data structures enable efficient CPU-based solution of network topology, conveniently storing and querying relationships between vertices and triangles, and easily constructing complete topological information. However, when migrating to the GPU for processing, the situation becomes significantly more complex. GPU architectures are designed with a focus on parallel computing, and their hardware characteristics impose certain limitations on data structure processing. Currently, the industry generally uses two-dimensional arrays to simulate map structures, simulating the storage and access of key-value pairs through cleverly designed array indexes. However, this approach has significant drawbacks and consumes a large amount of storage space. GPU memory resources are relatively expensive and have limited capacity. Such a large-scale waste of storage space is undoubtedly unreasonable. This not only increases hardware costs but can also lead to inability to operate properly when processing large-scale mesh data due to insufficient video memory, significantly limiting the application scope and efficiency of GPU-based mesh processing algorithms.

[0070] The present invention implements the topology solution process on the GPU, wherein Table 1 contains the data structures used in the topology solution process. The method performs well in terms of space complexity, which is only O(N) (N is the number of points), greatly optimizing the video memory usage. During the calculation process, the multi-threaded parallel mechanism of the GPU is fully utilized to obtain results at an extremely fast speed. Under normal circumstances, the time consumed by the topology query operation is stable at O(1). Even in extreme cases, the complexity of some query operations is only O(N) (N is the number of points). The method has good stability and efficiency. Experiments show that the method can achieve an acceleration ratio of 2-3 orders of magnitude compared to traditional solutions.

[0071] According to Table 1, the topological structure of the edge can be obtained through E_Fir_Point and E_Sec_Point, the adjacent triangles of the edge can be obtained through E_Tris_Num_Prefix_Sum and E_Tris, and the adjacent edges of the vertex can be obtained through P_Edges_Num_Prefix_Sum and P_Edges. The specific algorithm is shown in 3.2.

[0072] 3.1: UV topology relationship construction process.

[0073] Input: Triangle mesh .

[0074] Output: Complete topology diagram.

[0075] Phase 1: Vertex-triangle adjacency construction.

[0076] Step 1: Parallelism statistics.

[0077] Step 1.1: Allocate the P_Tris_Num array and traverse the triangles in parallel using GPU threads .

[0078] Step 1.2: The three vertices of the vertices perform the atomic addition operation atomicAdd(&P_Tris_Num[v_i], 1).

[0079] Step 2: Parallel prefix sum calculation generates the P_Tris_Num_Prefix_Sum array.

[0080] Step 3: Fill in the adjacency relationship.

[0081] Step 3.1: Allocate P_Tris storage space based on the prefix sum.

[0082] Step 3.2: Multithreaded writing of triangle indices to ensure vertex-level data locality.

[0083] Phase 2: Edge topology derivation.

[0084] Step 1: Edge existence judgment, vertex The number of adjacent edges of satisfies: (3) Step 2: Parallel prefix sum calculation generates the P_Edges_Num_Prefix_Sum array.

[0085] Step 3: Fill in the adjacency relationship.

[0086] Step 3.1: Allocate P_Edges storage space based on the prefix sum.

[0087] Step 3.2: Multi-threaded writing to ensure vertex-level data locality.

[0088] Step 4: Compress and store undirected edges.

[0089] Step 4.1: Remove duplicate points from adjacent edges (keep edge), calculate Point_Start_Edge_Num.

[0090] Step 4.2: Parallel prefix sum calculation generates the P_Start_Edge_Num_Prefix_Sum array.

[0091] Step 4.3: Multi-threaded writing E_Fir_Point, E_Sec_Point.

[0092] Stage 3: Edge-triangle adjacency derivation.

[0093] Step 1: Parallel edge traversal, that is, for each edge , search The adjacent triangle set P_Tris[ ].

[0094] Step 2: Triangle inclusion detection.

[0095] Step 2.1: Perform bit mask detection. Step 2.2: Count the number of triangles that meet the conditions (1 / 2).

[0096] Stpe3: Topology label storage: Fill E_Tris_Num according to the detection results.

[0097] Step 4: Fill with E_Tris.

[0098] Table 1 Data structures used in the topology solution process 3.2 Topology query acceleration system.

[0099] The query operation has a unified paradigm. The edge existence query example is as follows Input: vertex v1, vertex v2.

[0100] Output: edge index / -1.

[0101] Step 1: Normalize the vertex order, that is, if the v1 index is greater than v2, swap the two values.

[0102] Step 2: Positioning candidate edge intervals.

[0103] Step 2.1: Prefix and fast addressing, that is, through the P_Start_Edge_Num_Prefix_Sum array, directly obtain the vertex v1 as the storage starting position start of the edge head node.

[0104] Step 2.2: Determine the range based on the out-degree. Calculate the end position end = start + out-degree according to the value of Point_Start_Edge_Num[v1] and lock the candidate edge set E_Sec_Point[start : end].

[0105] Step 3: Sequential search execution Starting at the starting position start, access the edges associated with vertex v1 one by one in memory order, read the secondary endpoint index of the current edge current_v2 = E_Sec_Point[idx], if current_v2 == v2: immediately return the current global index idx, and terminate the search when the traversal reaches position end. If no match is found, return -1.

[0106] 4. Parallel solution of frame field.

[0107] Since the frame field has four-fold symmetry, in order to ensure numerical continuity, we can use (4) Representing the frame field, this ensures numerical continuity while strictly meeting the physical requirements of four-fold symmetry. To further analyze the transformation characteristics of the field, the present invention decouples the frame field into two independent components: (5) in, (6) In order to make the frame field transformation smoother, we need to find a condition that satisfies, (7) of and Where S represents the surface. and , direct solution cannot guarantee (8) So the above formula (7) is adjusted to (9) By adjusting an adjustable parameter To adjust the extent of the requirements and The square of is 1. Represents an adjustable parameter, Represents a surface.

[0108] Reberol et al. used the heat diffusion system to derive equation (9) as (10) in, is the number of iterations, is a constant, which is related to the diagonal length and the shortest side length of the surface bounding box. is a diagonal matrix, represents the identity matrix, which can be expressed as (11) , For the edge The area of ​​the adjacent triangle, is the Crouzeix-Raviart stiffness matrix (thestiffness Crouzeix-Raviart matrices), It can be expressed as, (12) Where k represents the kth component. The initial condition matrix The construction follows the principle of boundary distinction. (13) where index .

[0109] In order to ensure the equation (8), each time we solve the equation (10), we need to , Mapped to the unit circle, the solution of Equation (10) often needs to be performed multiple times, which consumes the most time in the entire pricing field solution process.

[0110] 4.1: and Parallel computing.

[0111] According to formula (11), the matrix The construction of is closely related to the triangle area. Thanks to the GPU's advantage in parallel computing of geometric quantities, the efficient solution of the triangle area can be directly realized through GPU acceleration. Combined with the above-mentioned GPU-based edge adjacency triangle topology information, the matrix The individual elements of can be computed synchronously by parallel threads.

[0112] matrix The dimension of is 2N, where N is the number of edges. Figure 10 As shown, taking the calculation of a single edge as an example, the edge With edge In the matrix The elements in can be expressed as, (14) in, (15) like Figure 11 As shown, the edge can be calculated similarly side ,side side ,side side In the matrix Elements in .

[0113] Due to the different edge weights and There is no interdependence between the calculations, so the matrix The calculation steps are as follows: Input: mesh topology information, mesh geometry information V Output: Matrix Step 1: Traverse each edge in parallel and calculate .

[0114] Step 1.1: Find the adjacent triangles t1 and t2 of the current edge in E_Tris.

[0115] Step 1.2: In t1, t2, find the four adjacent edges of the current edge and calculate the angle .

[0116] Step 1.3: According to formula 15, we get .

[0117] Step 2: Traverse each edge in parallel and calculate Step 2.1: Calculate the local coordinate system in t1 and t2.

[0118] Step 2.2: According to Figure 11 , calculate the angle .

[0119] Step 3: According to formula 14, calculate 4.2: Parallel solution of systems of equations.

[0120] The solution of Equation 10 involves the solution of a system of linear equations. Reberol et al. The invariant properties can be solved using the direct method. However, when the matrix size is large, the direct method consumes a lot of memory, which is not very friendly for GPUs with limited video memory.

[0121] When solving sparse linear systems, the choice between direct and iterative methods requires a comprehensive consideration of matrix size and hardware characteristics. For small matrices (<100,000 degrees of freedom), direct methods (such as LU decomposition) provide exact solutions, and CPUs generally outperform due to their memory capacity advantages. However, when matrices scale to 100,000-1,000,000 degrees of freedom, GPU memory limitations (typically less than 48GB) prevent the intermediate storage required by direct methods from being fully resident, making CPU memory access efficiency more competitive. When the data size exceeds CPU memory capacity, direct methods become ineffective due to memory limitations. At this point, iterative methods (such as the preconditioned conjugate gradient method) become the only viable solution due to their O(N) memory complexity. Through optimized preconditioning techniques, iterative methods can achieve speedups of over 10 times compared to CPUs on GPU parallel architectures. This combination of memory efficiency and parallel computing advantages makes iterative methods a core technology path for GPU acceleration in large-scale frame field calculations.

[0122] This paper uses the parallel preconditioned conjugate gradient iterative method to solve Equation (10). The specific algorithm steps are as follows.

[0123] Input: triangular mesh and its topological relationship after segmentation.

[0124] Output: The intersection field on each edge.

[0125] Step 1: Calculate the matrix , .

[0126] Step 2: Calculate the matrix .

[0127] Step 3: Calculation , ,in is the diagonal length of the grid bounding box, is the shortest side length of the grid.

[0128] Step 4: Set the total number of loops N and the current number of loops , execute Step5-Step9 N times, each loop .

[0129] Step 5: Calculation .

[0130] Step 6: Parallel Update ,make .

[0131] Step 7: Parallel Computing .

[0132] Step 8: Solve using the parallel preconditioned conjugate gradient method .

[0133] Step 9: Parallel Update , ,if , then go to Step 7.

[0134] Step 10: Calculate the edge frame field.

[0135] 5. Model grid structure generation

[0136] After the frame field is acquired, the grid structure of the model can be generated based on the frame field. The method for generating the network structure can adopt any related technology in the prior art, and the present invention is not limited thereto.

[0137] Next, a specific experiment is used to illustrate the fast calculation of the surface frame field in the GPU-based parallel model grid structure generation method provided by the present invention.

[0138] This paper conducted experiments on computers primarily equipped with an Intel Core i7-10700K and an RTX 3080 GPU. Using surface models from the Mambo dataset, we tested the performance of three key steps: initial surface triangulation, parallel topology solution, and parallel frame field solution. First, we sampled 100,000, 1,000, and 1,000,000 points from the model and partitioned them into approximately 10, 100, and 1,000 faces, respectively. We compared the results using the CPU Delaunay method (CDT library), the GPU Delaunay method without UV stitching (GPU-DT library), and our proposed method. The results are shown in Table 2. When the number of points per face is small and the number of faces increases, the data transfer from the CPU (CDT) to the GPU consumes the most time, resulting in the GPU (GPU-DT) performing slower than the CPU.

[0139] Table 2 Delaunay time consumption of different methods The Delaunay mesh obtained by using GPU-DT and the method of the present invention is as follows Figure 12 As shown in the figure, (a) is the original model, (b) is the Delaunay mesh obtained using the GPU-DT method, and (c) is the Delaunay mesh obtained using the method of the present invention. Because GPU-DT does not handle cavities, the resulting mesh is inconsistent with the original surface model.

[0140] The topology of the triangular mesh after Delaunay is calculated using the proposed method and the traditional CPU-based method. The time consumption comparison results are shown in Table 3.

[0141] Table 3 Topology construction time comparison chart In the final solution phase, we tested models with approximately 100,000, 1 million, and 1,000 million edges. When the number of edges reached 100,000, the Eigen library's preconditioned conjugate gradient solver took over an hour to iterate and produce a result. The direct method on the CPU exceeded the device memory limit when the data size exceeded 1 million. The final results are shown in Table 4.

[0142] Table 4 Comparison of cross-field solution time Experiments were conducted using several models in the Mambo dataset to obtain some frame field calculation results, such as Figure 13 shown.

[0143] This paper successfully achieves efficient solution of surface model frame fields on a GPU architecture. By innovatively employing a parallel Delaunay algorithm based on UV stitching and a GPU parallel topology algorithm, the method achieves a speedup of one to two orders of magnitude compared to traditional methods. This achievement not only strengthens the technical foundation for parallel solution of surface model frame fields on GPUs but also provides a practical solution for migrating many surface model algorithms that previously ran only on CPUs to GPU platforms, significantly boosting computational efficiency in related fields.

[0144] The above embodiments are merely illustrative of the principles and effects of the present invention and are not intended to limit the present invention. Any skilled artisan may modify or alter the above embodiments without departing from the spirit and scope of the present invention. Therefore, the scope of protection of the present invention shall be as set forth in the appended claims.

Claims

1. A method for generating quadrilateral meshes of geometric models based on GPU parallelism, characterized in that: The method comprises: For a model composed of multiple faces, each GPU thread obtains the parametric coordinates of the face in the UV coordinate system in parallel, and performs spatial transformation on the parametric coordinates of the face by establishing the face's bounding box to obtain the face splicing result; Generate an original triangular mesh model based on the patch stitching result, and use each GPU thread to process each triangular unit in the original triangular mesh model in parallel to obtain a cavity-filtered triangular mesh and network geometry information of the triangular mesh, wherein the network geometry information includes vertices and triangles; Each GPU thread solves the network topology structure of the triangular mesh in parallel; wherein the network topology structure includes: the topological composition of the edge, the adjacent triangles of the edge, and the adjacent edges of the vertex; Each GPU thread performs parallel solution of the frame field based on the network geometry information and network topology information of the triangular mesh; The frame field is combined to generate a quadrilateral grid structure of a geometric model.

2. The method according to claim 1, characterized in that Each GPU thread obtains the parametric coordinates of the patch in the UV coordinate system in parallel, and performs spatial transformation on the parametric coordinates of the patch by establishing the bounding box of the patch to obtain the patch stitching result, including: Each GPU thread obtains the parameterized coordinates of all points in the patch in the UV coordinate system in parallel , Indicates the parameterized coordinate value of a point in the patch in the horizontal direction in the UV coordinate system. Indicates the parametric coordinate value of a point in the patch in the vertical direction in the UV coordinate system; For each face of the model, a bounding box is used to enclose the face; wherein the bounding box is the minimum rectangular frame that contains the parameterized plane of the face, and the extreme value coordinates of the bounding box include: horizontal extreme value U_min, horizontal extreme value U_max, vertical extreme value V_min and vertical extreme value Vmax; Each GPU thread performs vertical offset preprocessing based on the vertical extreme value Vmax and the patch processing order index in parallel on a patch basis to obtain the vertical reference value of the patch. Each GPU thread obtains the position of the corresponding patch in the patch queue according to the patch processing order index, and obtains the total offset of the preceding patches of the patch based on the vertical reference of all patches before the patch. Each GPU thread is based on the total offset of the previous patch of the patch, the parameterized coordinate value Perform vertical spatial transformation on all points in the patch and the vertical extreme value V_min to obtain the transformed parameterized coordinate values ​​of the corresponding points ; Get the parameterized coordinate values ​​of all points in the patch after spatial transformation .

3. The method according to claim 1, characterized in that The method for generating an original triangular mesh model based on a patch splicing result includes: GPU threads executing a Delaunay decomposition algorithm in parallel.

4. The method according to claim 1, wherein Utilize GPU threads to process each triangle unit in the original triangle mesh model in parallel to obtain a cavity-filtered triangle mesh and network geometry information of the triangle mesh, including: Each GPU thread calculates the centroid coordinates of different triangle units in parallel and constructs a centroid-triangle unit index mapping table; Construct a first ray along the extended diagonal direction of the bounding box, and select an external point as the end point of the first ray; Generate a second ray starting from the centroid of the triangular element and ending at the external point; Establishing a boundary edge quadtree space partitioning structure, and combining the boundary edge quadtree space partitioning structure with calculating the intersection detection between the second ray and the edge of the triangle unit to obtain an intersection count of the triangle unit; Perform Boolean judgment on the intersection count of each triangle unit to obtain the reserved flag bitmap of the triangle unit; Compressing the original triangular mesh model based on the retention flag bitmap of the triangular unit, and after searching the corresponding triangle index using the centroid-triangle index mapping table, deleting the triangles with the mark bit of 0 in the retention flag bitmap to obtain a triangular mesh after cavity filtering; Extract the network geometry information of the triangular mesh ;in, represents a vertex, Represents a triangle.

5. The method according to claim 1, wherein Each GPU thread solves the network topology structure of the triangular mesh in parallel, including: For each vertex in the triangle mesh , traverse the triangle in parallel through GPU threads To generate a P_Tris array, the P_Tris array is used to store the adjacent triangle index set of the vertex; Determine vertices based on the P_Tris array Is there an adjacent edge at this vertex? If there are adjacent edges, perform parallel prefix sum calculations based on the P_Tris array to generate the P_Edges_Num_Prefix_Sum array, which is used to store the prefix sums of the adjacent edges of the vertex; Allocate storage space for the P_Edges array based on the prefix sum of the adjacent edges, and write the vertices through multi-threaded writing. The adjacent edge index of the vertex is written into the P_Edges array, and the P_Edges array is used to store the adjacent edge index set of the vertex; Remove duplicates from the adjacent edges in the P_Edges array and write the deduplicated adjacent edge index set into the Point_Start_Edge_Num array; Perform parallel prefix sum calculation based on the Point_Start_Edge_Num array to generate the P_Start_Edge_Num_Prefix_Sum array, where the P_Start_Edge_Num_Prefix_Sum array is used to store the prefix sums of the endpoints of the adjacent edges; Based on the prefix sum of the adjacent edge endpoints, the adjacent edge endpoints are written into the E_Fir_Point array and the E_Sec_Point array respectively; wherein, for any adjacent edge, the index value of the adjacent edge endpoint in the E_Fir_Point array is less than the index value of the adjacent edge endpoint in the E_Sec_Point array; Perform parallel edge traversal on the triangle mesh and perform a parallel edge traversal on each edge. , retrieve the vertex The adjacent triangle set P_Tris[ ];in, Represents an edge Another vertex of Based on the adjacent triangle set P_Tris[ ] Perform triangle containment detection, and based on the detection results, store the number of adjacent triangles of the edge and the index of the adjacent triangles of the edge in the E_Tris_Num array and the E_Tris array respectively; The topological composition of the edge, the adjacent triangles of the edge, and the adjacent edges of the vertex are obtained through the E_Fir_Point array and E_Sec_Point array, the E_Tris_Num_Prefix_Sum array and E_Tris array, and the P_Edges_Num_Prefix_Sum array and P_Edges array respectively.

6. The method according to claim 5, characterized in that The vertex is realized by following the steps below and vertices The network topology query operation between: According to the vertex index value, the vertex and vertices Perform sequence normalization; Get the vertex through the P_Start_Edge_Num_Prefix_Sum array As the storage starting position of the edge head node, start; According to Point_Start_Edge_Num[ ] Calculate the end position end = start + out-degree, lock the candidate edge set E_Sec_Point[start : end]; Starting from the starting position start, access the vertices one by one in memory order Associated edge, read the secondary endpoint index current_ = E_Sec_Point[idx], if current_ == : Immediately returns the current global index idx and ends the search when traversal reaches position end; if no match is found, returns -1.

7. The method according to claim 5, characterized in that Each GPU thread performs parallel solution of the frame field based on the network geometry and topology information of the triangular mesh, including: Construct the optimization problem: ; Among them, the first component , the second component , represents the frame field, represents the surface, Indicates adjustable parameters; Using the heat diffusion system to derive the optimization problem, we obtain the equation: ;in, represents the number of iterations, Represents a constant related to the diagonal length and the shortest side length of the surface bounding box. is a diagonal matrix, the elements in the diagonal matrix , 、 Represents edges The area of ​​an adjacent triangle of Represents the stiffness matrix, condition matrix is A single-column matrix, Indicates the The first component in the iteration The Quantity, Indicates the The first component in the iteration The Components, initial condition matrix The construction follows the principle of boundary distinction; Calculate the diagonal matrix based on the mesh topology information and mesh geometry information and the stiffness matrix The value of each element in; The diagonal matrix and the stiffness matrix Substitute the values ​​of each element in into the equation to solve the conditional frame field.

8. The method according to claim 7, characterized in that Calculate the diagonal matrix based on the mesh topology information and mesh geometry information and the stiffness matrix The values ​​of each element in include: Based on the adjacent triangles of the edge, parallel threads are synchronized to obtain the diagonal matrix The value of each element in; Traverse each edge and find the adjacent triangle of the current edge based on the E_Tris array and adjacent triangles ; In adjacent triangles and adjacent triangles Find the four adjacent edges of the current edge and calculate the angle between the four adjacent edges and the current edge , angle , angle and angle Then, based on the angle , angle , angle and angle Calculate the weight of the current edge; Calculate adjacent triangles and adjacent triangles The local coordinate system of the current edge is calculated based on the local coordinate system. ; Based on the angle of the current edge and weights, calculate the stiffness matrix The value of the corresponding element in .

9. The method according to claim 7, characterized in that The diagonal matrix and the stiffness matrix Substitute the values ​​of each element in into the equation to solve the conditional frame field, including: Step 1: Calculate the matrix ; Step 2: Based on the diagonal length of the bounding box and the shortest side length of the grid , calculate the initial angle and the final angle XX ; Step 3: Set the total number of loops , current number of cycles , and execute Step 4-Step 8, each cycle ; Step 4: Based on the initial angle , final angle , Current cycle number and the total number of cycles , calculate the current angle ; Step 5: Based on the current angle , diagonal matrix and the stiffness matrix Parallel matrix updates , and let the number of iterations ; Step 6: Based on the current angle and the conditional matrix , calculate the matrix ; Step 7: Solve using the parallel preconditioned conjugate gradient method ; Step 8: Parallel Update , ; Step 9: If , then go to Step 7; where c represents a set threshold; Step 10: If , then calculate the edge frame field.

10. A GPU-based parallel geometric model quadrilateral mesh generation system, characterized in that: The system comprises: The UV stitching module is used to obtain the parametric coordinates of the facets in the UV coordinate system in parallel for models composed of multiple faces. It then performs spatial transformation on the parametric coordinates of the facets by establishing a bounding box for the facets to obtain the facet stitching result. A triangular mesh generation module is used to generate an original triangular mesh model based on the facet splicing result, and use each GPU thread to process each triangular unit in the original triangular mesh model in parallel to obtain a triangular mesh after cavity filtering and network geometry information of the triangular mesh, wherein the network geometry information includes vertices and triangles; A topology structure generation module is used for each GPU thread to solve the network topology structure of the triangular mesh in parallel; wherein the network topology structure includes: the topological composition of the edge, the adjacent triangles of the edge, and the adjacent edges of the vertex; The frame field generation module is used by each GPU thread to perform parallel solution of the frame field based on the network geometry information and network topology information of the triangular mesh; The quadrilateral mesh generation module is used to generate the quadrilateral mesh structure of the geometric model in combination with the frame field.

Images