Thin plate spline deformation parameter matrix calculation and linear interpolation method based on CUDA acceleration

Abstract

The application provides a CUDA acceleration-based thin-plate spline deformation parameter matrix calculation and linear interpolation method, which processes three-dimensional to-be-registered feature points and template feature points through three-dimensional affine transformation; the to-be-registered feature points and the template feature points processed by the affine transformation matrix are used to calculate a thin-plate spline deformation parameter matrix based on functions in cublas and cusolver in the CUDA architecture; the to-be-registered three-dimensional image coordinates are inversely transformed by the obtained thin-plate spline deformation parameter matrix, and three-linear interpolation is performed. Through parallelization of large matrix algebra operation and linear interpolation by using the parallel architecture of the GPU, fast calculation of thin-plate spline transformation deformation parameter matrix calculation and linear interpolation is realized, and the running efficiency of the image registration algorithm is effectively improved.

Description

Technical Field

[0001] The present invention relates to the field of image registration and parallel acceleration technology, and in particular to a thin plate spline deformation parameter matrix calculation and linear interpolation method based on CUDA acceleration. Background Art

[0002] Warping images to a standard coordinate space is crucial for many image-based computational tasks. However, for multidimensional and high-resolution images, accurate warping operations can be memory-intensive and computationally expensive. For high-throughput image analysis studies, such as brain mapping and registration projects, high-performance image transformation tools that are compatible with common image analysis pipelines are required.

[0003] Image registration can be achieved through a variety of methods, some of the most commonly used methods include Gaussian smoothing and heat diffusion, B-spline method and thin plate spline (TPS) method. Compared to TPS, the first three methods require additional parameters but usually produce less accurate warping results.

[0004] However, TPS often produces significant distortion when fitting images. Therefore, we combine affine transformations with TPS transformations to improve registration quality. However, due to the computational complexity of the TPS algorithm, when the number of control points is in the hundreds or thousands, the calculation process still takes a long time to complete. Furthermore, to obtain better registration results, hundreds or even thousands of iterations are typically performed to improve registration accuracy. This makes the entire process extremely time-consuming, and even some powerful workstations struggle to run it.

[0005] Therefore, how to achieve a deformation method that combines affine transformation and thin plate spline transformation to achieve fast deformation while ensuring accuracy has become a technical problem that needs to be solved urgently. Summary of the Invention

[0006] In view of this, the purpose of the embodiments of the present invention is to solve the defects of thin plate spline deformation and linear interpolation in the existing technology that the calculation time is long when the number of control points in the image registration task is large, and to provide a thin plate spline deformation parameter matrix calculation and linear interpolation method based on CUDA acceleration.

[0007] A first aspect of the present invention provides a thin plate spline deformation parameter matrix calculation and linear interpolation method based on CUDA acceleration, the method comprising:

[0008] S1, processing the three-dimensional feature points to be registered and the template feature points through three-dimensional affine transformation;

[0009] S2, using the feature points to be registered and the template feature points obtained in S1 and processed by the affine transformation matrix, and based on the functions in cublas and cusolver in the CUDA architecture, calculating the thin plate spline deformation parameter matrix;

[0010] S3, performing an inverse transformation of the coordinates of the three-dimensional image to be registered using the thin plate spline deformation parameter matrix obtained in S2, and performing trilinear interpolation.

[0011] Furthermore, the step S1, processing the three-dimensional feature points to be registered and the template feature points through three-dimensional affine transformation, includes:

[0012] Obtain input N groups of feature points to be registered and template feature points, and normalize the N groups of feature points to be registered and template feature points;

[0013] According to the affine transformation formula, the affine transformation matrix is ​​calculated;

[0014] Denormalize the affine transformation matrix to obtain the final affine transformation matrix, and then multiply the N groups of feature points to be registered by the final affine transformation matrix to obtain N groups of feature points to be registered after affine transformation;

[0015] Among them, the affine transformation formula is as follows:

[0016]

[0017] Furthermore, S2 uses the feature points to be registered and the template feature points obtained in S1 and processed by the affine transformation matrix, and calculates the thin plate spline deformation parameter matrix based on the functions in cublas and cusolver in the CUDA architecture, including:

[0018] S21, constructing a process matrix required for calculating the thin plate spline deformation parameter matrix based on the interpolation function of the three-dimensional thin plate spline transformation; the interpolation function formula is as follows:

[0019]

[0020] Where U is the basis function, which is defined as follows:

[0021]

[0022]

[0023] Using the feature points to be registered and the template feature points processed by the affine transformation matrix, the process matrices R, P, Y, and L are constructed:

[0024]

[0025]

[0026]

[0027]

[0028] Among them, the R matrix stores the basis function values ​​of N groups of feature points, and its size is N×N. The first column of the P matrix is ​​all 1, and the last three columns store N groups of template feature points, and its size is N×4. The L matrix is ​​composed of the R matrix, the P matrix, and the P transposed matrix. The remaining positions are 0, and its size is (N+4)×(N+4). The first N rows of the Y matrix store the feature points to be registered after being processed by the affine transformation matrix, and the last 4 rows are all 0. Its size is (N+4)×3;

[0029] S22, based on the functions in the CUDA cublas library and cusolver library, the pseudo-inverse matrix L is obtained by calculating its pseudo-inverse matrix L -1 ;

[0030] S23, according to the pseudo inverse matrix L obtained in S22 -1 , and the Y matrix, construct the thin plate spline deformation parameter matrix W, whose formula is as follows:

[0031]

[0032] Finally, the pseudo-inverse matrix L - Multiplying with the Y matrix yields the W matrix, whose size is (N+4)×3.

[0033] Further, S22, the functions in the CUDA-based cublas library and cusolver library calculate the pseudo-inverse of the L matrix to obtain the pseudo-inverse matrix L -1 ,include:

[0034] S221, based on the SVD decomposition function in the CUDA cusolver library, performs SVD decomposition on the L matrix and obtains the U matrix, S matrix and V matrix. T The matrix is ​​expressed as follows:

[0035] L=USV T

[0036] Among them, because the size of L matrix is ​​(N+4)×(N+4), the U matrix and V T The matrix is ​​an orthogonal matrix of (N+4)×(N+4), and the S matrix is ​​a diagonal matrix of (N+4)×(N+4);

[0037] S222, convert the L matrix into a one-dimensional array for storage, and pass it from the CPU to the GPU, and finally pass it to the SVD decomposition function in the cusolver library, and finally obtain the U matrix, the diagonal values ​​in the S matrix and the VT Matrix, and take the value from GPU to CPU;

[0038] S223, since the S matrix values ​​obtained by the SVD decomposition function in the cusolver library only have values ​​on its diagonal, it is padded to form a (N+4)×(N+4) diagonal matrix;

[0039] S224, since the SVD decomposition function in the cusolver library obtains the U matrix, the U matrix is ​​transposed and the U matrix is ​​obtained. T Matrix, whose size is (N+4)×(N+4);

[0040] S225, since the SVD decomposition function in the cusolver library yields V T Matrix, so for V T The matrix is ​​transposed to obtain the V matrix, whose size is (N+4)×(N+4);

[0041] S226, for the diagonal matrix S of size (N+4)×(N+4), use the function in the cusolver library in CUDA to obtain the inverse matrix S of S -1 ;

[0042] S227, calculate the pseudo-inverse of the L matrix using the SVD method, and the formula is as follows:

[0043] L -1 =VS -1 U T

[0044] Among them, V matrix, S -1 Matrix and U T The matrices have been obtained in the previous steps. Finally, the functions in the cublas library based on CUDA are used to accelerate the multiplication of the three matrices, and finally the pseudo-inverse matrix L of the L matrix is ​​obtained. -1 , whose size is (N+4)×(N+4).

[0045] Furthermore, the step S3 performs an inverse transformation of the coordinates of the three-dimensional image to be registered using the thin plate spline deformation parameter matrix obtained in step S2, and performs trilinear interpolation, including:

[0046] S31, using CUDA to convert the three-dimensional serial traversal of the voxels of the three-dimensional image to be registered into a combination of one-dimensional serial traversal and parallel traversal;

[0047] S32, the 3D template image data, deformation parameter matrix W, and template feature points are passed to the GPU. Through the CUDA kernel function constructed in the previous step, in the requested GPU thread, the voxel coordinates in the 3D image to be registered are inversely transformed by using the thin plate spline deformation parameter matrix. The specific formula is as follows:

[0048]

[0049] Multiplying the voxel coordinates in the three-dimensional image to be registered with the thin plate spline deformation parameter matrix W obtained in step 23 to obtain the corresponding voxel coordinates in the template registration image;

[0050] S33, based on the three-dimensional voxel coordinates obtained in the previous step, the corresponding three-dimensional image to be registered is interpolated using trilinear interpolation method, and the formula is as follows:

[0051]

[0052] For a point (x, y, z), x0 represents the previous point closest to x, and x1 represents the next point closest to x. Similarly, for y0, y1, z0, z1, x d ,y d , z d , represents the corresponding weight;

[0053] Interpolation in the X-axis direction is performed using the following formula:

[0054]

[0055] Where V[x0,y0,z0] represents the value of the voxel in the three-dimensional template image;

[0056] Then interpolate the Y-axis direction, the formula is as follows:

[0057]

[0058] Interpolation is performed along the Z axis, and the formula is as follows:

[0059] c=c0(1-z d )+c1z d

[0060] Finally, c is the value of a voxel in the three-dimensional image to be registered. Through the GPU parallel architecture, a combination of one-dimensional serial traversal and parallel traversal is used to calculate the value of each voxel in the three-dimensional image to be registered, thereby completing the linear interpolation of the three-dimensional image to be registered.

[0061] Furthermore, the S31, using CUDA to complete the conversion of the three-dimensional serial traversal of the voxels of the three-dimensional image to be registered into a combination of one-dimensional serial traversal and parallel traversal, includes:

[0062] By utilizing CUDA's parallel design architecture, the serial loops for the first two dimensions of the three-dimensional image to be registered are placed on the GPU and simplified to one-dimensional parallel processing:

[0063] First, assign N to the CUDA kernel function z ×N Y threads, where N z Refers to the Z-axis size of the three-dimensional image to be registered, where N Y Refers to the Y-axis size of the three-dimensional image to be registered, and then the blockDim in the kernel function is fixed to a two-dimensional block with a size of 32*32;

[0064] Each thread has a unique row and column number, and the row number ranges from 0 to N. z , whose column numbers range from 0 to N Y , corresponding to the Z and Y axis sizes of the three-dimensional image data;

[0065] Then put the CUDA kernel function into a size of 0 to N x In the cycle, N x Refers to the X-axis size of the three-dimensional image to be registered; finally, each thread processes the work of each voxel in the three-dimensional image to be registered. Assuming that this is the k-th layer loop, the thread number (m, n) processes the voxel corresponding to the X-axis coordinate k, Y-axis coordinate n, and Z-axis coordinate m in the three-dimensional voxel, that is, each thread number corresponds to a voxel in the three-dimensional image to be registered.

[0066] In addition, the second aspect of the present invention provides a storage medium storing a computer program; the program is loaded and executed by a processor to implement the above-mentioned CUDA-accelerated thin plate spline deformation parameter matrix calculation and linear interpolation method steps.

[0067] In the solution of the present invention, three-dimensional feature points to be registered and template feature points are processed through a three-dimensional affine transformation; the feature points to be registered and template feature points processed by the affine transformation matrix are used, and based on the functions in cublas and cusolver in the CUDA architecture, a thin plate spline deformation parameter matrix is ​​calculated; the obtained thin plate spline deformation parameter matrix is ​​used to perform an inverse transformation of the coordinates of the three-dimensional image to be registered, and trilinear interpolation is performed. Compared with the existing technology, the thin plate spline transformation deformation parameter matrix calculation and linear interpolation method of the present invention based on CUDA acceleration utilizes the parallel architecture of the GPU to parallelize large matrix algebraic operations and linear interpolation, achieving rapid calculation of the thin plate spline transformation deformation parameter matrix and linear interpolation, effectively improving the operating efficiency of the image registration algorithm. BRIEF DESCRIPTION OF THE DRAWINGS

[0068] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for use in the embodiments. It should be understood that the following drawings only illustrate certain embodiments of the present invention and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other relevant drawings can be obtained based on these drawings without paying any creative work.

[0069] Figure 1 This is a flow chart of a thin plate spline deformation parameter matrix calculation and linear interpolation method based on CUDA acceleration disclosed in an embodiment of the present invention;

[0070] Figure 2 It is a schematic structural diagram of an electronic device disclosed in an embodiment of the present invention. DETAILED DESCRIPTION

[0071] Example embodiments will now be described more fully with reference to the accompanying drawings. However, example embodiments can be implemented in many forms and should not be construed as limited to the examples set forth herein; rather, these embodiments are provided so that this application will be thorough and complete and will fully convey the concepts of the example embodiments to those skilled in the art.

[0072] In addition, described feature, structure or characteristic can be combined in one or more embodiments in any suitable manner.In the following description, many specific details are provided so as to provide a full understanding of the embodiments of the present application. However, it will be appreciated by those skilled in the art that the technical scheme of the present application can be put into practice without one or more of the specific details, or other methods, components, devices, steps etc. can be adopted. In other cases, known methods, devices, implementations or operations are not shown or described in detail to avoid blurring the various aspects of the application.

[0073] The block diagrams shown in the accompanying drawings are merely functional entities and do not necessarily correspond to physically separate entities. That is, these functional entities may be implemented in software, in one or more hardware modules or integrated circuits, or in different networks and / or processor devices and / or microcontroller devices.

[0074] The flowcharts shown in the accompanying drawings are for illustrative purposes only and do not necessarily include all contents and operations / steps, nor must they be executed in the order described. For example, some operations / steps may be decomposed, while others may be combined or partially combined. Therefore, the actual execution order may vary depending on the actual situation.

[0075] It should be noted that the “plurality” mentioned in this article refers to two or more.

[0076] The following is a detailed description of the implementation details of the technical solution of the embodiment of the present application:

[0077] See also Figure 1 , Figure 1 This is a flow chart of a thin plate spline deformation parameter matrix calculation and linear interpolation method based on CUDA acceleration disclosed in an embodiment of the present invention. Figure 1 As shown, a thin plate spline deformation parameter matrix calculation and linear interpolation method based on CUDA acceleration in an embodiment of the present invention includes:

[0078] S1, processing the three-dimensional feature points to be registered and the template feature points through three-dimensional affine transformation;

[0079] Furthermore, the step S1, processing the three-dimensional feature points to be registered and the template feature points through three-dimensional affine transformation, includes:

[0080] Obtain input N groups of feature points to be registered and template feature points, and normalize the N groups of feature points to be registered and template feature points;

[0081] According to the affine transformation formula, the affine transformation matrix is ​​calculated;

[0082] Denormalize the affine transformation matrix to obtain the final affine transformation matrix, and then multiply the N groups of feature points to be registered by the final affine transformation matrix to obtain N groups of feature points to be registered after affine transformation;

[0083] Among them, the affine transformation formula is as follows:

[0084]

[0085] S2, using the feature points to be registered and the template feature points obtained in S1 and processed by the affine transformation matrix, and based on the functions in cublas and cusolver in the CUDA architecture, calculating the thin plate spline deformation parameter matrix;

[0086] Specifically, in this embodiment, the thin plate spline transformation deformation parameter matrix is ​​calculated: the functions in the CUDA-based cublas library and cusolver library are used to perform SVD decomposition on the L matrix, and the U matrix, S matrix and V matrix are obtained. T Matrix, then solve the pseudo-inverse matrix of the L matrix, and finally obtain the thin plate spline transformation deformation parameter matrix W. Because SVD decomposition, large matrix multiplication, and large matrix inversion operations of large matrices are computationally time-consuming on the CPU, by leveraging the hardware advantages of the GPU itself, the SVD decomposition, matrix multiplication, and matrix inversion processes of large matrices are accelerated, saving a considerable amount of time compared to traditional CPU algorithms.

[0087] Furthermore, S2 uses the feature points to be registered and the template feature points obtained in S1 and processed by the affine transformation matrix, and calculates the thin plate spline deformation parameter matrix based on the functions in cublas and cusolver in the CUDA architecture, including:

[0088] S21, constructing a process matrix required for calculating the thin plate spline deformation parameter matrix based on the interpolation function of the three-dimensional thin plate spline transformation; the interpolation function formula is as follows:

[0089]

[0090] Where U is the basis function, which is defined as follows:

[0091]

[0092]

[0093] Using the feature points to be registered and the template feature points processed by the affine transformation matrix, the process matrices R, P, Y, and L are constructed:

[0094]

[0095]

[0096]

[0097]

[0098] Among them, the R matrix stores the basis function values ​​of N groups of feature points, and its size is N×N. The first column of the P matrix is ​​all 1, and the last three columns store N groups of template feature points, and its size is N×4. The L matrix is ​​composed of the R matrix, the P matrix, and the P transposed matrix. The remaining positions are 0, and its size is (N+4)×(N+4). The first N rows of the Y matrix store the feature points to be registered after being processed by the affine transformation matrix, and the last 4 rows are all 0. Its size is (N+4)×3;

[0099] S22, based on the functions in the CUDA cublas library and cusolver library, the pseudo-inverse matrix L is obtained by calculating its pseudo-inverse matrix L -1 ;

[0100] S23, according to the pseudo inverse matrix L obtained in S22 -1 , and the Y matrix, construct the thin plate spline deformation parameter matrix W, whose formula is as follows:

[0101]

[0102] Finally, the pseudo-inverse matrix L -Multiplying with the Y matrix yields the W matrix, whose size is (N+4)×3.

[0103] Further, S22, the functions in the CUDA-based cublas library and cusolver library calculate the pseudo-inverse of the L matrix to obtain the pseudo-inverse matrix L -1 ,include:

[0104] S221, based on the SVD decomposition function in the CUDA cusolver library, performs SVD decomposition on the L matrix and obtains the U matrix, S matrix and V matrix. T The matrix is ​​expressed as follows:

[0105] L=USV T

[0106] Among them, because the size of L matrix is ​​(N+4)×(N+4), the U matrix and V T The matrix is ​​an orthogonal matrix of (N+4)×(N+4), and the S matrix is ​​a diagonal matrix of (N+4)×(N+4);

[0107] S222, convert the L matrix into a one-dimensional array for storage, and pass it from the CPU to the GPU, and finally pass it to the SVD decomposition function in the cusolver library, and finally obtain the U matrix, the diagonal values ​​in the S matrix and the V T Matrix, and take the value from GPU to CPU;

[0108] S223, since the S matrix values ​​obtained by the SVD decomposition function in the cusolver library only have values ​​on its diagonal, it is padded to form a (N+4)×(N+4) diagonal matrix;

[0109] S224, since the SVD decomposition function in the cusolver library obtains the U matrix, the U matrix is ​​transposed and the U matrix is ​​obtained. T Matrix, whose size is (N+4)×(N+4);

[0110] S225, since the SVD decomposition function in the cusolver library yields V T Matrix, so for V T The matrix is ​​transposed to obtain the V matrix, whose size is (N+4)×(N+4);

[0111] S226, for the diagonal matrix S of size (N+4)×(N+4), use the function in the cusolver library in CUDA to obtain the inverse matrix S of S -1 ;

[0112] S227, calculate the pseudo-inverse of the L matrix using the SVD method, and the formula is as follows:

[0113] L -1 =VS -1 U T

[0114] Among them, V matrix, S -1 Matrix and U T The matrices have been obtained in the previous steps. Finally, the functions in the cublas library based on CUDA are used to accelerate the multiplication of the three matrices, and finally the pseudo-inverse matrix L of the L matrix is ​​obtained. -1 , whose size is (N+4)×(N+4).

[0115] S3, performing an inverse transformation of the coordinates of the three-dimensional image to be registered using the thin plate spline deformation parameter matrix obtained in S2, and performing trilinear interpolation.

[0116] Specifically, this embodiment uses the thin plate spline deformation parameter matrix W obtained by S2 to obtain the three-dimensional template image coordinates corresponding to the three-dimensional image coordinates to be registered, and uses trilinear interpolation to interpolate the three-dimensional image to be registered. Since voxel traversal of a three-dimensional image is time-consuming in the CPU, the GPU parallel architecture is used to convert the three-dimensional serial traversal into a method that combines one-dimensional serial traversal with parallel traversal, which can significantly save computing time. This embodiment places the inverse transformation of the three-dimensional registered image coordinates and the linear interpolation of the image to be registered into the designed CUDA kernel function, and uses a method that combines one-dimensional serial traversal with parallel traversal to save computing time and is easy to introduce.

[0117] Furthermore, the step S3 performs an inverse transformation of the coordinates of the three-dimensional image to be registered using the thin plate spline deformation parameter matrix obtained in step S2, and performs trilinear interpolation, including:

[0118] S31, using CUDA to convert the three-dimensional serial traversal of the voxels of the three-dimensional image to be registered into a combination of one-dimensional serial traversal and parallel traversal;

[0119] S32, the 3D template image data, deformation parameter matrix W, and template feature points are passed to the GPU. Through the CUDA kernel function constructed in the previous step, in the requested GPU thread, the voxel coordinates in the 3D image to be registered are inversely transformed by using the thin plate spline deformation parameter matrix. The specific formula is as follows:

[0120]

[0121] Multiplying the voxel coordinates in the three-dimensional image to be registered with the thin plate spline deformation parameter matrix W obtained in step 23 to obtain the corresponding voxel coordinates in the template registration image;

[0122] S33, based on the three-dimensional voxel coordinates obtained in the previous step, the corresponding three-dimensional image to be registered is interpolated using trilinear interpolation method, and the formula is as follows:

[0123]

[0124] For a point (x, y, z), x0 represents the previous point closest to x, and x1 represents the next point closest to x. Similarly, for y0, y1, z0, z1, x d ,y d , z d , represents the corresponding weight;

[0125] Interpolation in the X-axis direction is performed using the following formula:

[0126]

[0127] Where V[x0,y0,z0] represents the value of the voxel in the three-dimensional template image;

[0128] Then interpolate the Y-axis direction, the formula is as follows:

[0129]

[0130] Interpolation is performed along the Z axis, and the formula is as follows:

[0131] c=c0(1-z d )+c1z d

[0132] Finally, c is the value of a voxel in the three-dimensional image to be registered. Through the GPU parallel architecture, a combination of one-dimensional serial traversal and parallel traversal is used to calculate the value of each voxel in the three-dimensional image to be registered, thereby completing the linear interpolation of the three-dimensional image to be registered.

[0133] Furthermore, the S31, using CUDA to complete the conversion of the three-dimensional serial traversal of the voxels of the three-dimensional image to be registered into a combination of one-dimensional serial traversal and parallel traversal, includes:

[0134] By utilizing CUDA's parallel design architecture, the serial loops for the first two dimensions of the three-dimensional image to be registered are placed on the GPU and simplified to one-dimensional parallel processing:

[0135] First, assign N to the CUDA kernel function z ×N Y threads, where N z Refers to the Z-axis size of the three-dimensional image to be registered, where N Y Refers to the Y-axis size of the three-dimensional image to be registered, and then the blockDim in the kernel function is fixed to a two-dimensional block with a size of 32*32;

[0136] Each thread has a unique row number and column number, and its row number ranges from 0 to N z , whose column numbers range from 0 to N Y , corresponding to the Z and Y axis sizes of the three-dimensional image data;

[0137] Then put the CUDA kernel function into a size of 0 to N x In the cycle, N x Refers to the X-axis size of the three-dimensional image to be registered; finally, each thread processes the work of each voxel in the three-dimensional image to be registered. Assuming that this is the k-th layer loop, the thread number (m, n) processes the voxel corresponding to the X-axis coordinate k, Y-axis coordinate n, and Z-axis coordinate m in the three-dimensional voxel, that is, each thread number corresponds to a voxel in the three-dimensional image to be registered.

[0138] Furthermore, in this embodiment, the running time results of the thin plate spline deformation parameter matrix calculation and linear interpolation method based on CUDA acceleration of this embodiment are tested under 1742 groups of control points, as shown in Table 1.

[0139] Table 1

[0140] Points: 1742 groups First run time Second run time Third run time Average running time Not accelerated 2176.11s 2189.26s 2148.43s 2171.27s After acceleration 103.80s 101.34s 100.90s 102.01s

[0141] As can be seen from Table 1, based on the platform with a CPU of i5-10400F and a GPU of RTX-3060, the experimental results were compared with the unaccelerated experimental results when the number of control points was 1742, and the acceleration ratio can reach about 21 times.

[0142] In addition, the present application also discloses an electronic device, such as Figure 2 As shown, the electronic device includes: one or more processors, a memory, the memory is used to store one or more computer programs; the computer program is configured to be executed by the one or more processors, and the program includes steps for executing the thin plate spline deformation parameter matrix calculation and linear interpolation method based on CUDA acceleration as described above.

[0143] In addition, an embodiment of the present application further provides a storage medium storing a computer program; the program is loaded and executed by a processor to implement the above-mentioned CUDA-accelerated thin plate spline deformation parameter matrix calculation and linear interpolation method steps.

[0144] Those skilled in the art will appreciate that the units and algorithm steps of each example described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, computer software, or a combination of the two. In order to clearly illustrate the interchangeability of hardware and software, the above description has generally described the composition and steps of each example according to function. Whether these functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Professional and technical personnel can use different methods to implement the described functions for each specific application, but such implementation should not be considered to be beyond the scope of the present invention.

[0145] The specific implementation methods described above further illustrate the objectives, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above description is only a specific implementation method of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A thin plate spline deformation parameter matrix calculation and linear interpolation method based on CUDA acceleration, characterized in that: The method comprises: S1, processing the three-dimensional feature points to be registered and the template feature points through three-dimensional affine transformation; S2, using the feature points to be registered and the template feature points obtained in S1 and processed by the affine transformation matrix, and based on the functions in cublas and cusolver in the CUDA architecture, calculating the thin plate spline deformation parameter matrix; The S2 uses the feature points to be registered and the template feature points obtained in the S1 and processed by the affine transformation matrix, and calculates the thin plate spline deformation parameter matrix based on the functions in cublas and cusolver in the CUDA architecture, including: S21, constructing a process matrix required for calculating the thin plate spline deformation parameter matrix based on the interpolation function of the three-dimensional thin plate spline transformation; the interpolation function formula is as follows: Where U is the basis function, which is defined as follows: Using the feature points to be registered and the template feature points processed by the affine transformation matrix, the process matrices R, P, Y, and L are constructed: Among them, the R matrix stores the basis function values ​​of N groups of feature points, and its size is N×N. The first column of the P matrix is ​​all 1, and the last three columns store N groups of template feature points, and its size is N×4. The L matrix is ​​composed of the R matrix, the P matrix, and the P transposed matrix. The remaining positions are 0, and its size is (N+4)×(N+4). The first N rows of the Y matrix store the feature points to be registered after being processed by the affine transformation matrix, and the last 4 rows are all 0, and its size is (N+4)×3; S22, based on the functions in the CUDA cublas library and cusolver library, the pseudo-inverse matrix L is obtained by calculating its pseudo-inverse matrix L -1 ; S23, according to the pseudo inverse matrix L obtained in S22 -1 , and the Y matrix, construct the thin plate spline deformation parameter matrix W, whose formula is as follows: Finally, the pseudo-inverse matrix L - Multiplying with the Y matrix yields the W matrix, whose size is (N+4)×3; S3, performing an inverse transformation of the coordinates of the three-dimensional image to be registered using the thin plate spline deformation parameter matrix obtained in S2, and performing trilinear interpolation.

2. The thin plate spline deformation parameter matrix calculation and linear interpolation method based on CUDA acceleration according to claim 1 is characterized in that: The step S1, processing the three-dimensional feature points to be registered and the template feature points through three-dimensional affine transformation, includes: Obtain input N groups of feature points to be registered and template feature points, and normalize the N groups of feature points to be registered and template feature points; According to the affine transformation formula, the affine transformation matrix is ​​calculated; Denormalize the affine transformation matrix to obtain the final affine transformation matrix, and then multiply the N groups of feature points to be registered by the final affine transformation matrix to obtain N groups of feature points to be registered after affine transformation; Among them, the affine transformation formula is as follows:

3. The thin plate spline deformation parameter matrix calculation and linear interpolation method based on CUDA acceleration according to claim 2, characterized in that: S22, the functions in the CUDA-based cublas library and cusolver library calculate the pseudo-inverse of the L matrix to obtain the pseudo-inverse matrix L -1 ,include: S221, based on the SVD decomposition function in the CUDA cusolver library, performs SVD decomposition on the L matrix and obtains the U matrix, S matrix and V matrix. T The matrix is ​​expressed as follows: L=USV T Among them, because the size of L matrix is ​​(N+4)×(N+4), the U matrix and V T The matrix is ​​an orthogonal matrix of (N+4)×(N+4), and the S matrix is ​​a diagonal matrix of (N+4)×(N+4); S222, convert the L matrix into a one-dimensional array for storage, and pass it from the CPU to the GPU, and finally pass it to the SVD decomposition function in the cusolver library, and finally obtain the U matrix, the diagonal values ​​in the S matrix and the V T Matrix, and take the value from GPU to CPU; S223, since the S matrix values ​​obtained by the SVD decomposition function in the cusolver library only have values ​​on its diagonal, it is padded to form a (N+4)×(N+4) diagonal matrix; S224, since the SVD decomposition function in the cusolver library obtains the U matrix, the U matrix is ​​transposed and the U matrix is ​​obtained. T Matrix, whose size is (N+4)×(N+4); S225, since the SVD decomposition function in the cusolver library yields V T Matrix, so for V T The matrix is ​​transposed to obtain the V matrix, whose size is (N+4)×(N+4); S226, for the diagonal matrix S of size (N+4)×(N+4), use the function in the cusolver library in CUDA to obtain the inverse matrix S of S -1 ; S227, calculate the pseudo-inverse of the L matrix using the SVD method, and the formula is as follows: L -1 =VS -1 U T Among them, V matrix, S -1 Matrix and U T The matrices have been obtained in the previous steps. Finally, the functions in the cublas library based on CUDA are used to accelerate the multiplication of the three matrices, and finally the pseudo-inverse matrix L of the L matrix is ​​obtained. -1 , whose size is (N+4)×(N+4).

4. The thin plate spline deformation parameter matrix calculation and linear interpolation method based on CUDA acceleration according to claim 1 or 3, characterized in that: The step S3 performs an inverse transformation of the coordinates of the three-dimensional image to be registered using the thin plate spline deformation parameter matrix obtained in step S2, and performs trilinear interpolation, including: S31, using CUDA to convert the three-dimensional serial traversal of the voxels of the three-dimensional image to be registered into a combination of one-dimensional serial traversal and parallel traversal; S32, the 3D template image data, deformation parameter matrix W, and template feature points are passed to the GPU. Through the CUDA kernel function constructed in the previous step, in the requested GPU thread, the voxel coordinates in the 3D image to be registered are inversely transformed by using the thin plate spline deformation parameter matrix. The specific formula is as follows: Multiplying the voxel coordinates in the three-dimensional image to be registered with the thin plate spline deformation parameter matrix W obtained in step 23 to obtain the corresponding voxel coordinates in the template registration image; S33, based on the three-dimensional voxel coordinates obtained in the previous step, the corresponding three-dimensional image to be registered is interpolated using trilinear interpolation method, and the formula is as follows: For a point (x, y, z), x0 represents the previous point closest to x, and x1 represents the next point closest to x. Similarly, for y0, y1, z0, z1, x d ,y d , z d , represents the corresponding weight; Interpolation is performed in the X-axis direction, and the formula is as follows: Where V[x0,y0,z0] represents the value of the voxel in the three-dimensional template image; Then interpolate the Y-axis direction, the formula is as follows: Interpolation is performed along the Z axis, and the formula is as follows: c=c0(1-z d )+c1z d Finally, c is the value of a voxel in the three-dimensional image to be registered. Through the GPU parallel architecture, a combination of one-dimensional serial traversal and parallel traversal is used to calculate the value of each voxel in the three-dimensional image to be registered, thereby completing the linear interpolation of the three-dimensional image to be registered.

5. The thin plate spline deformation parameter matrix calculation and linear interpolation method based on CUDA acceleration according to claim 4 is characterized in that: The step S31, using CUDA to convert a three-dimensional serial traversal of voxels of the three-dimensional image to be registered into a combination of a one-dimensional serial traversal and a parallel traversal, includes: By utilizing CUDA's parallel design architecture, the serial loops for the first two dimensions of the three-dimensional image to be registered are placed on the GPU and simplified to one-dimensional parallel processing: First, assign N to the CUDA kernel function z ×N Y threads, where N z Refers to the Z-axis size of the three-dimensional image to be registered, where N Y Refers to the Y-axis size of the three-dimensional image to be registered, and then the blockDim in the kernel function is fixed to a two-dimensional block with a size of 32*32; Each thread has a unique row number and column number, and its row number ranges from 0 to N z , whose column numbers range from 0 to N Y , corresponding to the Z and Y axis sizes of the three-dimensional image data; Then put the CUDA kernel function into a size of 0 to N x In the cycle, N x Refers to the X-axis size of the three-dimensional image to be registered; finally, each thread processes the work of each voxel in the three-dimensional image to be registered. Assuming that this is the k-th layer loop, the thread number (m, n) processes the voxel corresponding to the X-axis coordinate k, Y-axis coordinate n, and Z-axis coordinate m in the three-dimensional voxel, that is, each thread number corresponds to a voxel in the three-dimensional image to be registered.

6. An electronic device, comprising: One or more processors, a memory, the memory being used to store one or more computer programs; characterized in that the computer program is configured to be executed by the one or more processors, and the program includes steps for executing the CUDA-accelerated thin plate spline deformation parameter matrix calculation and linear interpolation method as described in any one of claims 1 to 5.

7. A storage medium storing a computer program; characterized in that: The program is loaded and executed by a processor to implement the steps of the thin plate spline deformation parameter matrix calculation and linear interpolation method based on CUDA acceleration as described in any one of claims 1 to 5.

Images