GPU parallel processing method for reverse solution of product NURBS surface mapping

Abstract

The present invention discloses a GPU parallel processing method for inverse solution of product NURBS surface mapping. A given NURBS surface, UV plane, and three-dimensional point group of the product are input; the intermediate matrix of each segment of the NURBS surface is calculated in parallel using the GPU and stored in the GPU; a subdivision tree of the NURBS surface is constructed in parallel using the GPU and stored in the GPU; and each three-dimensional point in the three-dimensional point group to be mapped and the subdivision tree are processed in parallel using the GPU to obtain the two-dimensional feature point corresponding to each three-dimensional point on the UV plane. The present invention utilizes two methods, surface subdivision method and Newton iterative numerical method, to solve the problem of feature point mapping of the product NURBS surface, optimize the representation of the NURBS surface, and improve the convergence speed of the product surface algorithm while ensuring the acquisition of the global optimal solution.

Description

Technical Field

[0001] The present invention relates to the processing of a product NURBS surface in the CAD field, and in particular to a GPU parallel processing method for reversely solving mapping feature points of a product NURBS surface. Background Art

[0002] In the field of computer-aided geometric design, NURBS surfaces have become the most commonly used way to represent three-dimensional objects in modern CAD systems due to their powerful representation capabilities and high flexibility.

[0003] A common problem in product NURBS surface research is the inverse solution of NURBS surfaces, that is, calculating the 2D (u, v) parameters corresponding to the 3D points (x, y, z) on the NURBS surface. This problem plays a vital role in solving problems such as NURBS surface fitting of point clouds, collision detection, physical simulation, and interference avoidance in CAD / CAM systems. When processing large and complex 3D models, the inverse solution operation is often time-consuming, and parallelizing the algorithm is a common method to increase computational speed.

[0004] Furthermore, the commonly used representation of NURBS surfaces is the Cox-deBoor basis function recursive representation. This requires multiple recursions to calculate points and partial derivatives on the surface. The resulting recursive function calls consume a large amount of registers and cache. Due to the limited register and cache capacity of each block in the GPU, deep recursion cannot be supported. Furthermore, cache allocation and release operations within a block cannot be parallelized and must be processed serially by the cache manager. Therefore, the basis function recursive representation is inherently unsuitable for parallel processing on the GPU.

[0005] Existing parallel methods for inverse problem solving simply deploy the code to the GPU, without considering the characteristics of GPU programming. Data storage and use lack commonality, and the parallelism of the code is not guaranteed, thus only leveraging a limited portion of the parallel performance. Summary of the Invention

[0006] In response to the shortcomings of the background technology, the present invention proposes a parallel algorithm for inverse solution of NURBS surfaces that combines the surface subdivision method and the Gauss-Newton iterative method. The matrix representation of the NURBS surface is used to calculate the three-dimensional points and partial derivatives on the surface. The intermediate matrix required for the matrix representation is calculated in advance and stored as a constant memory, which improves the public nature of the data and avoids the large amount of register occupation and frequent request and release of the cache in the recursive method.

[0007] The present invention first calculates the intermediate matrix corresponding to each node interval according to the NURBS surface degree, control points and node vectors input by the user, then constructs the surface subdivision tree in parallel on the GPU, and then performs hierarchical traversal from the root node of the subdivision tree for each point in the array in parallel. During the hierarchical traversal process, bounding box detection and Newton iteration are used alternately to find the target two-dimensional point.

[0008] To achieve the above-mentioned purpose, the technical solution adopted by the present invention is specifically the following steps:

[0009] 1) Input the NURBS surface G, the uv plane to be mapped, and the 3D point set P to be mapped. The NURBS surface G is described by a bivariate piecewise rational function S(u,v), where u and v represent two orthogonal directions on the uv plane. The NURBS surface G is divided into multiple segments by segmentation points, and adjacent segments are separated by segmentation points.

[0010] 2) Use the GPU to parallelly calculate the intermediate matrix of each segment of the NURBS surface G and save it in the constant memory area of ​​the GPU;

[0011] 3) Use the GPU to construct the subdivision tree T of the NURBS surface G in parallel and save it in the constant memory area of ​​the GPU;

[0012] 4) Use GPU to parallelize each 3D point P in the 3D point group P to be mapped i Processing is performed between the subdivision tree T to obtain each three-dimensional point P i The corresponding two-dimensional feature points on the uv plane.

[0013] The step 2) is specifically as follows:

[0014] For the NURBS surface, solve the matrix form of the following partial derivatives to obtain the coefficient matrix N of the NURBS surface in the u and v directions respectively: u and N v :

[0015] S p,q (u,v)=UN u H(N v ) T V T

[0016]

[0017]

[0018]

[0019]

[0020] Among them, S p,q (u,v) represents the three-dimensional point on the NURBS surface corresponding to the two-dimensional parameter point (u,v), and Represents the first-order partial derivative of the NURBS surface corresponding to the two-dimensional parameter point (u, v) in the u and v directions, p and q represent the order of the NURBS surface G in the u and v directions respectively, n and m represent the number of control points of the NURBS surface G in the u and v directions respectively, U and V are the segment point vectors of the NURBS surface in the u and v directions respectively, N u and N v They are the coefficient matrices of the NURBS surface in the u and v directions, H is the matrix composed of all control points on the NURBS surface, u p+1 、v q+1 Represent the p+1th and q+1th segment point values ​​in the segment point vectors U and V respectively;

[0021] Then calculate the intermediate matrix N u H(N v ) T And saved in the constant memory area of ​​the GPU.

[0022] In step 2), the coefficient matrix N u and N v The value of each element in is obtained by iterating N times according to the following recursive function:

[0023]

[0024]

[0025] In the above formula, r and c are the row index and column index of the coefficient matrix, p and q are the orders of the NURBS surface G in the u and v directions, i and j are the segment point interval indexes of the current guess parameter point on U and V during the Newton iteration process; N u (r,c) and N v (r,c) represents the coefficient matrix N u and N v The recursive function value of the element in row r and column c, a, b, k, l represent the recursive parameters, and their initial values ​​are p, q, i, j, respectively. Represents the difference between the i+1th segmentation point and the i-th segmentation point on the segmentation point vector U, that is, u i+1 -u i , Represents the difference u between the k+pa segmentation point and the kth segmentation point on the segmentation point vector U k+p-a -u k ,the remaining The same can be said.

[0026] The guessed parameter points are intermediate two-dimensional parameter points generated during the Newton iteration process.

[0027] The intermediate matrix in the NURBS surface matrix representation is frequently accessed in subsequent Newton iteration and surface tessellation calculations. The data is read frequently and does not need to be modified during the operation process. The present invention stores this part of the pre-calculated intermediate matrix data in the constant memory area of ​​the GPU, which can be broadcast to each thread in the thread warp within a single memory access cycle. The reading speed of the constant memory can almost reach the speed of the first-level cache.

[0028] In step 3, the steps of constructing the subdivision tree T of the NURBS surface G in parallel using the GPU are as follows:

[0029] 3.1) Create a blank current layer node list List and next layer node list NextList, which are used to store the tree nodes of the current layer and the next layer respectively;

[0030] 3.2) Each tree node in the subdivision tree T is a 1*1 rectangular area located on the uv plane. The 1*1 rectangular area located at the origin of the uv plane is stored as the root node R in the current layer node list List;

[0031] 3.3) Parallel processing of each tree node in the current layer node list List, performing subdivision operations on each tree node, and storing the new tree nodes generated by subdivision as the tree nodes of the next layer in the next layer node list NextList;

[0032] 3.4) If the next layer node list NextList is not empty, the current layer node list List and the next layer node list NextList are swapped. After the swap, the next layer node list NextList is cleared, and the process returns to step 3.3) to continue the iterative process.

[0033] If the next layer node list NextList is empty, the iterative process ends, thereby obtaining the tree nodes of each level of the subdivision tree T and saving them in the constant memory area of ​​the GPU.

[0034] The subdivision operation is specifically to calculate the change in the 1*1 rectangular area of ​​the tree node in the u and v directions according to the following formula:

[0035]

[0036]

[0037] Where, ΔS u and ΔS vRespectively represent the changes in the NURBS surface in the u and v directions, u0 and v0 represent the coordinate values ​​of the lower left corner of the rectangular area in the u and v directions, u1 and v1 represent the coordinate values ​​of the upper right corner of the rectangular area in the u and v directions, respectively. and Represents the first-order partial derivative of a two-dimensional parameter point on the NURBS surface in the u and v directions. The first-order partial derivative is calculated using the aforementioned NURBS surface matrix representation. e represents the number of samples on the NURBS surface, and i represents the i-th sample.

[0038] Then compare the changes in the 1*1 rectangular area of ​​the tree node in the u and v directions with the preset tolerance δ:

[0039] If the change of the 1*1 rectangular area of ​​the tree node in the u / v direction exceeds the tolerance δ, the 1*1 rectangular area of ​​the tree node is evenly divided into two in this direction, and the resulting rectangular area is used as a new tree node;

[0040] If the 1*1 rectangular area of ​​the tree node does not exceed the tolerance δ in the u / v direction, the 1*1 rectangular area of ​​the tree node is not evenly bisected in this direction, and no new tree node is generated.

[0041] In this way, if the change in a certain direction of the NURBS surface segment corresponding to the node in the tree is less than the tolerance δ, then no further subdivision is required in that direction. The subdivision tree determines in which direction of the UV plane to subdivide during the construction process.

[0042] During the construction of the tessellation tree, the present invention stores the dynamic list of the current layer tree nodes in the shared memory of each SM, avoiding direct reading of the global memory, reducing the memory read and write time, and improving the memory access bandwidth.

[0043] When the list length is insufficient, it is expanded, which avoids pre-allocation and occupies a large amount of memory. For the tessellation tree itself, the tree structure is stored in the form of an array to improve data locality, thereby improving cache hit rate and avoiding random access to memory space.

[0044] The step 4) is specifically as follows:

[0045] For each 3D point P in the 3D point set P to be mapped i , using GPU to perform level traversal in parallel starting from the root node of the subdivision tree T:

[0046] 4.1) First determine the current 3D point P i The level in the subdivision tree T, and then use the GPU in parallel for the three-dimensional point Pi The surface segment corresponding to the 1*1 rectangular area of ​​each tree node in the tree level establishes a bounding box and determines whether the bounding box contains the current 3D point P i , and make judgments:

[0047] If the bounding box contains the current 3D point P i , then the tree node corresponding to the bounding box is retained and used as a valid tree node;

[0048] If the bounding box does not contain the current 3D point P i , then discard the tree node corresponding to the bounding box;

[0049] 4.2) From P i Starting from the center point of the rectangular plane of each valid tree node in the current layer, it expands to the surrounding 1*1 rectangular area on the uv plane, and performs a limited number of Newton iterations in parallel to find a two-dimensional feature point on the uv plane. The two-dimensional feature point corresponds to a point on the NURBS surface and the current three-dimensional point P i The distance between them is less than the tolerance δ, and the three-dimensional points and first-order partial derivatives of the NURBS surface used in the Newton iteration process are calculated using the aforementioned matrix representation;

[0050] If a GPU thread finds the corresponding two-dimensional feature point to a point on the NURBS surface and the current three-dimensional point P i If the distance between them is less than the tolerance δ, the two-dimensional feature point is saved to the two-dimensional feature point group D;

[0051] If multiple threads of the GPU obtain multiple two-dimensional feature points corresponding to a point on the NURBS surface and the current three-dimensional point P i If the distance between them is less than the tolerance δ, then take the two-dimensional feature point with the smallest distance to the two-dimensional feature point group D;

[0052] 4.3) If the current 3D point P in step 4.2) i If all tree nodes under the layer of the current subdivision tree T cannot obtain a two-dimensional feature point less than the tolerance δ after a limited number of Newton iterations, then the current three-dimensional point P i Enter the next level of tree T and make the following judgments:

[0053] If the current 3D point P i If the tree node in the level of the subdivision tree T is not empty, return to step 4.2) for processing;

[0054] If the current 3D point P i The tree node in the level of the subdivision tree T is empty, then the current 3D point P i Not on the NURBS surface, for the current 3D point P iThe processing is completed, for the next three-dimensional point P i+1 Return to step 4.1) for processing;

[0055] 4.4) After the above steps are iterated, all three-dimensional points P are obtained. i Corresponding to the two-dimensional feature points on the uv plane.

[0056] In step 4.1), the bounding box formed by the control points of the NURBS surface will be used as the bounding box of the NURBS surface. The convex hull property of the NURBS surface indicates that the NURBS surface is contained in the convex hull formed by its control points. In the above subdivision process, the cuboid formed by the NURBS control points is simply used as the bounding box of the NURBS surface. Each time the algorithm enters the next layer, the algorithm compares each NURBS surface segment in the current layer of the subdivision tree with the target 3D point P. i Perform bounding box detection in parallel and eliminate the bounding box that does not contain P i The tree node.

[0057] The present invention first uses the GPU to parallelly calculate the coefficient matrix of the surface G in the u and v directions according to the NURBS surface G input by the user; then calculates the intermediate matrix and saves it in the constant memory area of ​​the GPU; then uses the GPU to parallelly construct the subdivision tree T of the NURBS surface G, which is also saved in the constant memory area of ​​the GPU; for each three-dimensional point P in the input three-dimensional point group P, i , and start the hierarchical traversal from the root node of the subdivision tree T in parallel; during the traversal, each surface segment is first connected to the three-dimensional point P in parallel. i The bounding box detection is performed, and then a certain number of Newton iterations are performed to finally obtain the target point within a certain error range. In this way, the present invention can quickly and robustly generate the inverse solution array D corresponding to the input point group P based on the input NURBS surface and the input point group P.

[0058] The method of the present invention is to project the NURBS surface of the product onto a plane to obtain accurately corresponding two-dimensional feature points, and then the two-dimensional feature points can be used for modeling, simulation and simulation analysis.

[0059] The product described in the present invention is a product such as a geometric kernel of a CAD / CAM software system, such as an aircraft, daily necessities, etc.

[0060] The beneficial effects of the present invention are: giving full play to the hardware architecture characteristics and programming characteristics of the GPU, while ensuring the robustness of the NURBS surface inverse solution, greatly optimizing the storage performance of the NURBS surface inverse solution, and greatly improving the speed of the NURBS surface inverse solution.

[0061] The present invention utilizes two methods, namely the surface subdivision method and the Newton iteration numerical method, to deal with the problem of mapping feature points of product NURBS surfaces. The Newton iteration utilizes the surface subdivision method to select initial values, and the surface subdivision method is supplemented by the Newton iteration to improve the accuracy of the results. The matrix representation is used to optimize the representation of the NURBS surface, thereby improving the convergence speed of the product surface algorithm while ensuring the acquisition of the global optimal solution.

[0062] The present invention solves the problem of achieving efficient reverse solution of NURBS surfaces using a GPU. As an important basic operation of the geometric kernel of modern CAD software systems, NURBS surface reverse solution can play an important role in solving problems such as NURBS surface fitting of point cloud models, collision detection, physical simulation, and interference avoidance in CAD / CAM systems. The algorithm used in the present invention is implemented on a GPU, making full use of the hardware advantages of GPU multi-threading. Combined with the programming characteristics of the GPU, it ensures the robustness of the NURBS surface reverse solution, greatly optimizes the storage performance of the NURBS surface reverse solution, and greatly improves the speed of the NURBS surface reverse solution. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 It is a flow chart of the present invention.

[0064] Figure 2 Schematic diagram of the surface input of the embodiment.

[0065] Figure 3 Schematic diagram of implementing tessellation on a parameter domain.

[0066] Figure 4 Schematic diagram of a bounding box of an embodiment.

[0067] Figure 5 Schematic diagram of Newton's iteration method in an embodiment. DETAILED DESCRIPTION

[0068] The present invention will be further described below with reference to the accompanying drawings and examples.

[0069] like Figure 1 As shown, the embodiments of the complete method according to the content of the present invention are as follows:

[0070] The specific product implemented is Figure 2 The aircraft model products shown in Model D. Other examples Figure 2 As shown in ABC.

[0071] 1) Input the NURBS surface G, the uv plane to be mapped, and the 3D point set P to be mapped. The NURBS surface G is described by a bivariate piecewise rational function S(u,v), where u and v represent two orthogonal directions on the uv plane. The NURBS surface G is divided into multiple segments by preset segmentation points. Adjacent segments are separated by segmentation points, and each segment constitutes a segmentation point interval.

[0072] The input product NURBS surface G includes the degree, control points and their weights, and segmentation points. The NURBS surface G is divided into multiple segments by the segmentation points, and adjacent segments are divided by segmentation points.

[0073] The orders of NURBS surface G in u and v directions are p and q respectively. The number of control points of NURBS surface G in u and v directions are n and m respectively. The control point P f,g Represented as the f-th control point in the u direction and the g-th control point in the v direction of the NURBS surface G, w f,g represents the control point P f,g The weight of the NURBS surface G in the u and v directions is represented as U and V respectively.

[0074] Each segment of the NURBS surface G in the u direction is correspondingly provided with n p-order B-spline basis functions. The s-th p-order B-spline basis function in the u direction is expressed as N s,p (u, v); Each segment of the NURBS surface G in the v direction is correspondingly set with m q-order B-spline basis functions, and the t-th q-order B-spline basis function in the v direction is expressed as N t,q (u,v). The B-spline basis function is used to calculate the weight coefficient of each control point.

[0075] The three-dimensional point group P to be mapped is represented by Pi, where Pi represents the i-th point in the three-dimensional point group P.

[0076] D represents the two-dimensional feature point group D output by the method, where Di represents the i-th two-dimensional feature point.

[0077] In this embodiment, the parameters of the NURBS surface are shown in the following table:

[0078] Testing model parameters

[0079]

[0080] The generation rule of the input point group P is as follows: first, a pair of valid parameter values ​​​​(u, v) are generated through a random number program, and then the three-dimensional point (x, y, z) corresponding to the parameter value on the NURBS surface is calculated. Finally, a random noise ε is added to (x, y, z). The result will be used as input data to test the algorithm. Each model will be tested with 65536 test cases.

[0081] 2) Use the GPU to parallelly calculate the intermediate matrix of each segment of the NURBS surface G and save it in the constant memory area of ​​the GPU;

[0082] For the NURBS surface, solve the matrix form of the following partial derivatives to obtain the coefficient matrix N of the NURBS surface in the u and v directions respectively: u and N v :

[0083] S p,q (u,v)=UN u H(N v ) T V T

[0084]

[0085]

[0086]

[0087]

[0088] Then calculate the intermediate matrix N u H(N v ) T And saved in the constant memory area of ​​the GPU.

[0089] Coefficient matrix N u and N v The value of each element in is obtained by iterating N times according to the following recursive function:

[0090]

[0091]

[0092] The computational efficiency of the NURBS surface matrix representation and the traditional recursive representation is compared. The corresponding three-dimensional coordinate points and two first-order partial derivatives are calculated respectively, and the time overhead of the two is recorded. The results are shown in the following table:

[0093] Quantitative comparison of different NURBS representations

[0094]

[0095] The above results show that matrix representation can significantly reduce the time overhead of calculating three-dimensional coordinate points and partial derivatives.

[0096] 3) Use the GPU to construct the subdivision tree T of the NURBS surface G in parallel and save it in the constant memory area of ​​the GPU. The specific steps are as follows:

[0097] 3.1) Create a blank current layer node list List and next layer node list NextList, which are used to store the tree nodes of the current layer and the next layer respectively;

[0098] 3.2) Each tree node in the subdivision tree T is a 1*1 rectangular area located on the uv plane. The 1*1 rectangular area located at the origin of the uv plane is stored as the root node R in the current layer node list List;

[0099] 3.3) Parallel processing of each tree node in the current layer node list List, performing subdivision operations on each tree node, and storing the new tree nodes generated by subdivision as the tree nodes of the next layer in the next layer node list NextList;

[0100] 3.4) If the next layer node list NextList is not empty, the current layer node list List and the next layer node list NextList are swapped. After the swap, the next layer node list NextList is cleared, and the process returns to step 3.3) to continue the iterative process.

[0101] If the next layer node list NextList is empty, the iterative process ends, thereby obtaining the tree nodes of each level of the subdivision tree T and saving them in the constant memory area of ​​the GPU.

[0102] The subdivision operation is to calculate the change in the 1*1 rectangular area of ​​the tree node in the u and v directions according to the following formula:

[0103] Specifically, two first-order partial derivatives of the NURBS surface are used and To estimate the change ΔS of the NURBS surface in the u and v directions u and ΔS v , the calculation formula is as follows:

[0104]

[0105]

[0106] Then compare the changes in the 1*1 rectangular area of ​​the tree node in the u and v directions with the preset tolerance δ:

[0107] If the change of the 1*1 rectangular area of ​​the tree node in the u / v direction exceeds the tolerance δ, the 1*1 rectangular area of ​​the tree node is evenly divided into two in this direction, and the resulting rectangular area is used as a new tree node;

[0108] If the 1*1 rectangular area of ​​the tree node does not exceed the tolerance δ in the u / v direction, the 1*1 rectangular area of ​​the tree node is not evenly bisected in this direction, and no new tree node is generated.

[0109] The schematic diagram of the subdivision process and the subdivision tree T of the embodiment is as follows Figure 3 As shown, the subdivision tree of the embodiment has 4 layers in total.

[0110] 4) Use GPU to parallelize each 3D point P in the 3D point group P to be mapped i Processing is performed between the subdivision tree T to obtain each three-dimensional point P i The corresponding two-dimensional feature points on the uv plane.

[0111] For each 3D point P in the 3D point set P to be mapped i , using GPU to perform level traversal in parallel starting from the root node of the subdivision tree T:

[0112] 4.1) First determine the current 3D point P i The level in the subdivision tree T, and then use the GPU in parallel for the three-dimensional point P i The bounding box is created for the surface segment corresponding to the 1*1 rectangular area of ​​each tree node in the tree level. Figure 4 Schematic diagram of the bounding box of the embodiment, judging whether the bounding box contains the current three-dimensional point P i , and make judgments:

[0113] If the bounding box contains the current 3D point P i , then the tree node corresponding to the bounding box is retained and used as a valid tree node;

[0114] If the bounding box does not contain the current 3D point P i , then discard the tree node corresponding to the bounding box;

[0115] 4.2) The bounding boxes created by the surface segments corresponding to the rectangular planes of different tree nodes may overlap. i Starting from the center point of the rectangular plane of each valid tree node in the current layer, it expands to the surrounding 1*1 rectangular area on the uv plane, and performs a limited number of Newton iterations in parallel to find a two-dimensional feature point on the uv plane. The two-dimensional feature point corresponds to a point on the NURBS surface and the current three-dimensional point P iThe distance between them is less than the tolerance δ, and the three-dimensional points and first-order partial derivatives of the NURBS surface used in the Newton iteration process are calculated using the aforementioned matrix representation;

[0116] The limited number of times is specifically set to no more than 10 steps.

[0117] If a GPU thread finds the corresponding two-dimensional feature point to a point on the NURBS surface and the current three-dimensional point P i If the distance between them is less than the tolerance δ, the two-dimensional feature point is saved to the two-dimensional feature point group D;

[0118] If multiple threads of the GPU obtain multiple two-dimensional feature points corresponding to a point on the NURBS surface and the current three-dimensional point P i If the distance between them is less than the tolerance δ, then take the two-dimensional feature point with the smallest distance to the two-dimensional feature point group D.

[0119] The Newton iteration schematic diagram of the embodiment is Figure 5 , generally three Newton iterations are enough to converge.

[0120] 4.3) If the current 3D point P in step 4.2) i If all tree nodes under the layer of the current subdivision tree T cannot obtain a two-dimensional feature point less than the tolerance δ after a limited number of Newton iterations, then the current three-dimensional point P i Enter the next level of tree T and make the following judgments:

[0121] If the current 3D point P i If the tree node in the level of the subdivision tree T is not empty, return to step 4.2) for processing;

[0122] If the current 3D point P i The tree node in the level of the subdivision tree T is empty, then the current 3D point P i Not on the NURBS surface, for the current 3D point P i The processing is completed, for the next three-dimensional point P i+1 Return to step 4.1) for processing, and the method ends;

[0123] 4.4) After the above steps are iterated, all three-dimensional points P are obtained. i Corresponding to the two-dimensional feature points on the uv plane.

[0124] In the embodiment, the comparison of the final acceleration effect with other commonly used algorithms is shown in the following table:

[0125] Table 3 Inverse solution calculation speed (unit: thousand points per second)

[0126]

[0127]

[0128] Among them, the abc algorithm represents:

[0129] a. Inverse solution algorithm based on NURBS surface basis function representation (single thread);

[0130] b. Use OpenMP to accelerate (16 threads) the inverse solution algorithm based on the NURBS surface basis function representation;

[0131] c. The present invention uses a GPU-accelerated inverse solution algorithm based on NURBS matrix representation.

[0132] To address the time-consuming inverse solution problem of large, complex 3D models in CAD / CAM software, this paper proposes a general parallel algorithm based on surface tessellation, Gauss-Newton iteration, and NURBS matrix representation to solve the inverse problem of NURBS surfaces. Through GPU parallel acceleration, this algorithm achieves computational speed improvements of up to two orders of magnitude compared to traditional recursive algorithms and demonstrates good robustness on test models.

[0133] The above specific embodiments are used to illustrate the present invention rather than to limit the present invention. Any modifications and changes made to the present invention within the spirit of the present invention and the protection scope of the claims shall fall within the protection scope of the present invention.

Claims

1. A GPU parallel processing method for inverse solution of product NURBS surface mapping, characterized in that the method comprises the following steps: 1) Input the NURBS surface G, the uv plane to be mapped, and the 3D point set P to be mapped. The NURBS surface G is described by a bivariate piecewise rational function S(u,v), where u and v represent two orthogonal directions on the uv plane. The NURBS surface G is divided into multiple segments by segmentation points, and adjacent segments are separated by segmentation points. 2) Use the GPU to parallelly calculate the intermediate matrix of each segment of the NURBS surface G and save it in the GPU; 3) Use the GPU to construct the subdivision tree T of the NURBS surface G in parallel and save it in the GPU; 4) Use GPU to parallelize each 3D point P in the 3D point group P to be mapped i Processing is performed between the subdivision tree T to obtain each three-dimensional point P i The corresponding two-dimensional feature points on the uv plane; The step 2) is specifically as follows: For the NURBS surface, solve the matrix form of the following partial derivatives to obtain the coefficient matrix N of the NURBS surface in the u and v directions respectively: u and N v : Among them, S p,q (u,v) represents the three-dimensional point on the NURBS surface corresponding to the two-dimensional parameter point (u,v), and Represents the first-order partial derivative of the NURBS surface corresponding to the two-dimensional parameter point (u, v) in the u and v directions, p and q represent the order of the NURBS surface G in the u and v directions respectively, n and m represent the number of control points of the NURBS surface G in the u and v directions respectively, U and V are the segment point vectors of the NURBS surface in the u and v directions respectively, N u and N v They are the coefficient matrices of the NURBS surface in the u and v directions, H is the matrix composed of all control points on the NURBS surface, u p+1 、v q+1 Represent the p+1th and q+1th segment point values ​​in the segment point vectors U and V respectively; Then calculate the intermediate matrix N u H(N v ) T And saved in the constant memory area of ​​the GPU; In step 3, the steps of constructing the subdivision tree T of the NURBS surface G in parallel using the GPU are as follows: 3.1) Create a blank current layer node list List and next layer node list NextList, which are used to store the tree nodes of the current layer and the next layer respectively; 3.2) Each tree node in the subdivision tree T is a 1*1 rectangular area located on the uv plane. The 1*1 rectangular area located at the origin of the uv plane is stored as the root node R in the current layer node list List; 3.3) Parallel processing of each tree node in the current layer node list List, performing subdivision operations on each tree node, and storing the new tree nodes generated by subdivision as the tree nodes of the next layer in the next layer node list NextList; 3.4) If the next layer node list NextList is not empty, the current layer node list List and the next layer node list NextList are swapped. After the swap, the next layer node list NextList is cleared, and the process returns to step 3.3) to continue the iterative process. If the next layer node list NextList is empty, the iterative process ends, thereby obtaining the tree nodes of each level of the subdivision tree T and saving them in the constant memory area of ​​the GPU.

2. The GPU parallel processing method for inverse solution of product NURBS surface mapping according to claim 1, characterized in that: In step 2), the coefficient matrix N u and N v The value of each element in is obtained by iterating N times according to the following recursive function: In the above formula, r and c are the row index and column index of the coefficient matrix, p and q represent the order of the NURBS surface G in the u and v directions, i and j represent the segment point interval index of the current guess parameter point on U and V during the Newton iteration process; N u (r,c) and N v (r,c) represents the coefficient matrix N u and N v The recursive function value of the element in row r and column c, a, b, k, l represent the recursive parameters, and their initial values ​​are p, q, i, j respectively. Represents the difference between the i+1th segmentation point and the i-th segmentation point on the segmentation point vector U, that is, u i+1 -u i , Represents the difference u between the k+pa segmentation point and the kth segmentation point on the segmentation point vector U k+p-a -u k , Represents the difference between the i+1th segmentation point and the i-th segmentation point on the segmentation point vector V, that is, u i+1 -u i , Represents the difference u between the k+pa segmentation point and the kth segmentation point on the segmentation point vector V k+p-a -u k .

3. The GPU parallel processing method for inverse solution of product NURBS surface mapping according to claim 1, characterized in that: The subdivision operation is specifically to calculate the change in the 1*1 rectangular area of ​​the tree node in the u and v directions according to the following formula: Where, ΔS u and ΔS v Respectively represent the changes in the NURBS surface in the u and v directions, u0 and v0 represent the coordinate values ​​of the lower left corner of the rectangular area in the u and v directions, u1 and v1 represent the coordinate values ​​of the upper right corner of the rectangular area in the u and v directions, respectively. and Represents the first-order partial derivative of a two-dimensional parameter point on the NURBS surface in the u and v directions. The first-order partial derivative is calculated using the aforementioned NURBS surface matrix representation. e represents the number of samples on the NURBS surface, and i represents the i-th sample. Then compare the changes in the 1*1 rectangular area of ​​the tree node in the u and v directions with the preset tolerance δ: If the change of the 1*1 rectangular area of ​​the tree node in the u / v direction exceeds the tolerance δ, the 1*1 rectangular area of ​​the tree node is evenly divided into two in this direction, and the resulting rectangular area is used as a new tree node; If the 1*1 rectangular area of ​​the tree node does not exceed the tolerance δ in the u / v direction, the 1*1 rectangular area of ​​the tree node is not evenly bisected in this direction, and no new tree node is generated.

4. The GPU parallel processing method for inverse solution of product NURBS surface mapping according to claim 1, characterized in that: The step 4) is specifically as follows: For each 3D point P in the 3D point set P to be mapped i , using GPU to perform level traversal in parallel starting from the root node of the subdivision tree T: 4.1) First determine the current 3D point P i The level in the subdivision tree T, and then use the GPU in parallel for the three-dimensional point P i The surface segment corresponding to the 1*1 rectangular area of ​​each tree node in the tree level establishes a bounding box and determines whether the bounding box contains the current 3D point P i , and make judgments: If the bounding box contains the current 3D point P i , then the tree node corresponding to the bounding box is retained and used as a valid tree node; If the bounding box does not contain the current 3D point P i , then discard the tree node corresponding to the bounding box; 4.2) From P i Starting from the center point of the rectangular plane of each valid tree node in the current layer, it expands to the surrounding 1*1 rectangular area on the uv plane, and performs a limited number of Newton iterations in parallel to find a two-dimensional feature point on the uv plane. The two-dimensional feature point corresponds to a point on the NURBS surface and the current three-dimensional point P i The distance between them is less than the tolerance δ; If a GPU thread finds the corresponding two-dimensional feature point to a point on the NURBS surface and the current three-dimensional point P i If the distance between them is less than the tolerance δ, the two-dimensional feature point is saved to the two-dimensional feature point group D; If multiple threads of the GPU obtain multiple two-dimensional feature points corresponding to a point on the NURBS surface and the current three-dimensional point P i If the distance between them is less than the tolerance δ, then take the two-dimensional feature point with the smallest distance to the two-dimensional feature point group D; 4.3) If the current 3D point P in step 4.2) i If all tree nodes under the layer of the current subdivision tree T cannot obtain a two-dimensional feature point less than the tolerance δ after a limited number of Newton iterations, then the current three-dimensional point P i Enter the next level of tree T and make the following judgments: If the current 3D point P i If the tree node in the level of the subdivision tree T is not empty, return to step 4.2) for processing; If the current 3D point P i The tree node in the level of the subdivision tree T is empty, then the current 3D point P i Not on the NURBS surface, for the current 3D point P i The processing is completed, for the next three-dimensional point P i+1 Return to step 4.1) for processing; 4.4) After the above steps are iterated, all three-dimensional points P are obtained. i Corresponding to the two-dimensional feature points on the uv plane.

Images