CUDA-accelerated method for calculating deformation parameters of smooth thin plate splines

Abstract

The present invention relates to a method for calculating deformation parameters of a smooth thin plate spline based on CUDA acceleration, which solves the defect of long calculation time of smooth thin plate spline deformation when the number of control points is large in the image registration task compared with the prior art. The present invention comprises the following steps: processing of three-dimensional feature points to be registered and template feature points; performing QR decomposition on the matrix of the points to be registered; and calculating deformation parameters of the smooth thin plate spline. The present invention utilizes the parallel architecture of the GPU to parallelize large matrix algebraic operations, realizes fast calculation of deformation parameters of the smooth thin plate spline, and effectively improves the operating efficiency of the image registration algorithm.

Description

Technical Field

[0001] The invention relates to the technical fields of image registration and parallel acceleration, and in particular to a method for calculating deformation parameters of a smooth thin plate spline based on CUDA acceleration. Background Art

[0002] Warping images to a standard coordinate space is essential for many image computing-related tasks. However, for multi-dimensional and high-resolution images, accurate deformation operations are expensive in terms of computer memory consumption and computation time. For high-throughput image analysis studies, such as brain mapping projects, high-performance image warping 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 with TPS, the first three methods require additional parameters, but usually produce less accurate warping results. The goal of TPS deformation is to solve a function f such that f(Pi) = Pi' (1≤i≤n) and the bending energy function is minimized, while other points on the image can also be well corrected by interpolation. The deformation function can be imagined as bending a thin steel plate so that the steel plate passes through a given n points. The deformation formula can be expressed as follows:

[0004]

[0005] where u j is the deformation method, C 0 is the template point, A is the affine component, w k is the non-affine component, U(r)=r 2 log(r) is the TPS kernel function, p i is the coordinate of the template point, and p is the coordinate of the point to be registered.

[0006] However, TPS often suffers from large distortion when fitting images, and unrealistic deformations may occur locally. Wahba et al. proposed the smoothed thin plate spline deformation method (STPS) in 1990 to solve these problems by adding a smoothing component. However, when the number of control points is in the hundreds or thousands (which is common in current biological microscopy image registration tasks), it still takes a long time to complete the calculation process.

[0007] Moreover, at present, in order to obtain better registration results, the registration accuracy is usually improved by iterating hundreds or even thousands of times, which makes the entire process extremely time-consuming and even difficult to run on some very powerful workstations.

[0008] Therefore, how to achieve a smooth thin plate spline deformation method that can achieve rapid deformation while ensuring accuracy has become a technical problem that needs to be solved urgently. Summary of the invention

[0009] The purpose of the present invention is to solve the defect of long calculation time when the number of control points of smooth thin plate spline deformation in image registration task is large (hundreds or even thousands) in the prior art, and to provide a smooth thin plate spline deformation parameter calculation method based on CUDA acceleration to solve the above problem.

[0010] In order to achieve the above object, the technical solution of the present invention is as follows:

[0011] A method for calculating deformation parameters of a smooth thin plate spline based on CUDA acceleration comprises the following steps:

[0012] Processing of 3D feature points to be registered and template feature points: Obtain N groups of feature points to be registered and template feature points as input, calculate the distance between each feature point to be registered and each template feature point, and put them into a matrix U of size N×N, put the 3D coordinates of the feature points to be registered into a matrix X of size N×N, and put the 3D coordinates of the template feature points into a matrix Y of size N×4;

[0013] Perform QR decomposition on the registration point matrix: Use the functions and orthogonal extension kernel functions in CUDA-based cusolver to perform QR decomposition on the matrix X to obtain the Q matrix and R matrix;

[0014] Calculation of deformation parameters of smooth thin plate spline: Using matrices U and Y and matrices Q and R obtained in the previous step, a multiplication kernel function suitable for calculation is designed in CUDA, and based on the functions in CUDA's cublas, the affine component and non-affine component of the smooth thin plate spline deformation are calculated respectively to obtain the deformation parameters of the smooth thin plate spline.

[0015] The QR decomposition of the registration point matrix comprises the following steps:

[0016] The QR decomposition function in the CUDA-based cusolver library performs QR decomposition on the matrix X to obtain the Q matrix and R matrix, which are expressed as follows:

[0017]

[0018] At this time, the Q matrix has only the first four columns with values ​​and is orthogonal, and the remaining columns are 0. The R matrix has only the first four rows with values, and its size is N×4. Q1 is the first four columns of Q, and its size is N×4. Q2 is the part that needs to be expanded, and its size is N×(N-4). R1 ​​is the first four columns and the first four rows of R, and its size is 4×4.

[0019] The extended orthogonal kernel function designed based on CUDA completes the extended orthogonalization of the matrix Q, starting from the fifth column; when filling the i-th column, the following steps are used, i ≥ 5:

[0020] Find the sum of each row of Q and then get its square root to get the row number k of the row with the smallest value;

[0021] Pass Q to the GPU and allocate N threads, corresponding to each element in the i-th column;

[0022] The calculation content of the nth thread is: the nth row of the first column to the i-1th column of Q multiplied by the first to i-1th elements of the kth row;

[0023] After all threads are executed, the initial value of the expanded i-th column is obtained, and the obtained result is taken out from the GPU and put into a matrix T of size N×1;

[0024] Add 1 to the kth element of T, normalize T, and fill it into the i-th column of Q;

[0025] The orthogonal Q matrix is ​​calculated.

[0026] The calculation of the deformation parameters of the smooth thin plate spline comprises the following steps:

[0027] The following formula, CUDA-based cublas library function, and CUDA-based matrix multiplication kernel function are used to calculate the affine component and the non-affine component. The calculation formula is as follows:

[0028]

[0029]

[0030] in B=C+λI, W of size N×4 is the non-affine component, λ is the smoothing coefficient, I of size (N-4)×(N-4) is the identity matrix, and A of size 4×4 is the affine component;

[0031] Calculate W:

[0032] Done with CUDA The multiplication of U and Q2 results in a matrix C of size (N-4)×(N-4);

[0033] The result of the previous step is added to the I matrix multiplied by the smoothing coefficient to obtain a matrix B of size (N-4)×(N-4), and the inverse matrix of B is obtained using the function in the cusolver library based on CUDA;

[0034] Use CUDA to complete Q2 and B -1 , The product of gets the matrix D of size N×N,

[0035] First, transform the input matrices Q2 and B -1 As a one-dimensional array, allocate N×(N-4) threads and set the blockDim in the kernel function to a fixed two-dimensional block size of 32*32;

[0036] Each thread has a unique row and column number, with row numbers ranging from 0 to N and column numbers ranging from 0 to N-4, which is consistent with the matrix Q2B -1 The size corresponds to

[0037] Each thread calculates the matrix Q2B -1 An element of the thread number (m, n) thread calculation matrix Q2B -1 The element at row m and column n,

[0038] The method is to multiply the mth row of Q2 by B -1 Get the element at the corresponding position in the nth column;

[0039] Allocate N×N threads again and set the blockDim in the kernel function to a fixed two-dimensional block size of 32*32;

[0040] Each thread has a unique row and column number, with row numbers ranging from 0 to N and column numbers ranging from 0 to N, which is similar to the matrix The size corresponds to

[0041] Each thread calculates the matrix An element of the thread calculation matrix of thread number (m, n) The element in the mth row and nth column, using the Q2B method -1 Multiply the mth row of The nth column of gets the element at the corresponding position, where The mth row of is the mth column of Q2. The row index is converted to Q2 column index, thereby reducing one input; finally, the matrix D is obtained;

[0042] The multiplication of D and Y is completed on the CPU to obtain the non-affine component W;

[0043] The calculation of A is completed on the CPU to obtain the affine component, and thus all the deformation parameters of the smooth thin plate spline are obtained.

[0044] The calculation of the matrix C comprises the following steps:

[0045] Convert the input matrices U and Q2 into one-dimensional arrays, allocate (N-4)×N threads in the GPU, and set the blockDim in the kernel function to a fixed two-dimensional block size of 32*32;

[0046] Set each thread to have a unique row and column number, with the row number ranging from 0 to N-4 and the column number ranging from 0 to N, which is the same as the matrix The size corresponds to

[0047] Set the thread calculation matrix of thread number (m, n) The element at row m and column n,

[0048] The method is The mth row of U is multiplied by the nth column of U, where The mth row of is the mth column of Q2. The row index is converted to Q2 column index, thus reducing one input;

[0049] Allocate (N-4)×(N-4) threads again, and set blockDim in the kernel function to a fixed two-dimensional block size of 32*32;

[0050] Set each thread to have a unique row and column number, with the row number ranging from 0 to N-4 and the column number ranging from 0 to N-4, which is the same as the matrix The size corresponds to

[0051] Set each thread to calculate the matrix An element of the thread calculation matrix of thread number (m, n) The element at row m and column n,

[0052] The method is Multiply the mth row of U by the nth column of U to get the element at the corresponding position;

[0053] Finally, we get the matrix C.

[0054] Beneficial Effects

[0055] Compared with the prior art, the CUDA-accelerated smooth thin plate spline deformation parameter calculation method of the present invention utilizes the parallel architecture of GPU to parallelize large matrix algebraic operations, realizes fast calculation of smooth thin plate spline deformation parameters, and effectively improves the operating efficiency of the image registration algorithm.

[0056] Experiments have shown that this method achieves a speedup ratio of about 9 times on a platform with a CPU of i7-10700f and a GPU of RTX-2080, greatly shortening the deformation time. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] Figure 1It is a method sequence diagram of the present invention. DETAILED DESCRIPTION

[0058] In order to have a further understanding and recognition of the structural features and the effects achieved by the present invention, a preferred embodiment and accompanying drawings are used for detailed description as follows:

[0059] like Figure 1 As shown, the method for calculating deformation parameters of a smooth thin plate spline based on CUDA acceleration described in the present invention comprises the following steps:

[0060] The first step is to process the 3D feature points to be registered and the template feature points.

[0061] Get the input N groups of feature points to be registered and template feature points, calculate the distance between each feature point to be registered and each template feature point, and put them into a matrix U of size N×N, put the three-dimensional coordinates of the feature points to be registered into a matrix X of size N×N, and put the three-dimensional coordinates of the template feature points into a matrix Y of size N×4.

[0062] The second step is to perform QR decomposition on the matrix of feature points to be registered: using the functions and orthogonal extension kernel functions in cusolver based on CUDA, perform QR decomposition on the matrix X to obtain the Q matrix and the R matrix. By analyzing the time-consuming part of the algorithm and putting it into the GPU for calculation, the corresponding Q matrix and R matrix can be quickly calculated by taking advantage of the hardware advantages of the GPU itself, which can save a lot of time compared to the traditional CPU algorithm. The specific steps are as follows:

[0063] (1) Based on the QR decomposition function in the cusolver library of CUDA, the matrix X is decomposed into Q matrix and R matrix, which are expressed as follows:

[0064]

[0065] At this time, only the first four columns of the Q matrix have values ​​and are regular orthogonal, and the remaining columns are 0. Only the first four rows of the R matrix have values, and the size is N×4. Q1 is the first four columns of Q, and the size is N×4. Q2 is the part that needs to be expanded, and the size is N×(N-4). R1 ​​is the first four columns and the first four rows of R, and the size is 4×4.

[0066] (2) Based on the extended orthogonal kernel function designed by CUDA, the extended orthogonal matrix Q is completed, starting from the fifth column; when filling the i-th column, the following steps are adopted, i ≥ 5:

[0067] A1) Find the sum of each row of Q and then get its square root to get the row number k of the row with the smallest value;

[0068] A2) Pass Q to the GPU and allocate N threads to each element of the i-th column;

[0069] The calculation content of the nth thread is: the nth row of the first column to the i-1th column of Q multiplied by the first to i-1th elements of the kth row;

[0070] After all threads are executed, the initial value of the expanded i-th column is obtained, and the obtained result is taken out from the GPU and put into a matrix T of size N×1;

[0071] A3) Add 1 to the kth element of T, normalize T, and fill it into the i-th column of Q.

[0072] (3) Calculate and obtain the orthogonal Q matrix.

[0073] The third step is to calculate the deformation parameters of the smooth thin plate spline: using the matrices U and Y and the matrices Q and R1 obtained in the previous step, a multiplication kernel function suitable for calculation is designed in CUDA, and based on the functions in CUDA's cublas, the affine component and non-affine component of the smooth thin plate spline deformation are calculated respectively to obtain the deformation parameters of the smooth thin plate spline. Since large matrix multiplication takes a long time to calculate in the CPU, and the calculation of each element in the matrix does not interfere with each other during calculation, it can be calculated in parallel, thereby saving a lot of calculation time. The present invention designs a kernel function according to the algorithm formula to put the matrix multiplication into the kernel function, and reduces the input by modifying the data index, saving data transfer time. And the matrix inversion is completed based on the function in CUDA's cublas, which is simple and convenient. The specific steps are as follows:

[0074] (1) The calculation of affine components and non-affine components is completed using the following formula and the CUDA-based cublas library function and the matrix multiplication kernel function designed based on CUDA. The calculation formula is as follows:

[0075]

[0076]

[0077] in B=C+λI, W of size N×4 is the non-affine component, λ is the smoothing coefficient, I of size (N-4)×(N-4) is the identity matrix, and A of size 4×4 is the affine component.

[0078] (2) Calculate W:

[0079] A1) Using CUDA The multiplication of U and Q2 obtains a matrix C of size (N-4)×(N-4). To calculate the matrix C, first convert U and Q2 row by row into a one-dimensional array and pass it to the GPU. Allocate (N-4)×(N-4) threads in the GPU, and set blockDim in the kernel function to 32*32.

[0080] The calculation of matrix C includes the following steps:

[0081] A11) converts the input matrices U and Q2 into one-dimensional arrays, allocates (N-4)×N threads in the GPU, and sets the blockDim in the kernel function to a fixed two-dimensional block size of 32*32;

[0082] A12) Set each thread to have a unique row number and column number, with the row number ranging from 0 to N-4 and the column number ranging from 0 to N, which is the same as the matrix The size corresponds to

[0083] A13) Set the thread calculation matrix of thread number (m, n) The element at row m and column n,

[0084] The method is The mth row of U is multiplied by the nth column of U, where The mth row of is the mth column of Q2. The row index is converted to Q2 column index, thus reducing one input;

[0085] A14) reallocate (N-4)×(N-4) threads and set the blockDim in the kernel function to a fixed two-dimensional block size of 32*32;

[0086] A15) Set each thread to have a unique row number and column number, with the row number ranging from 0 to N-4 and the column number ranging from 0 to N-4, which is the same as the matrix The size corresponds to

[0087] A16) Set each thread to calculate the matrix An element of the thread calculation matrix of thread number (m, n) The element at row m and column n,

[0088] The method is Multiply the mth row of U by the nth column of U to get the element at the corresponding position;

[0089] A17) Finally, we get the matrix C.

[0090] A2) The result of the previous step is added to the I matrix multiplied by the smoothing coefficient to obtain a matrix B of size (N-4)×(N-4), and the inverse matrix of B is obtained by using the function in the cusolver library based on CUDA;

[0091] A3) Use CUDA to complete Q2 and B -1 , The product of gets the matrix D of size N×N,

[0092] First, transform the input matrices Q2 and B -1 As a one-dimensional array, allocate N×(N-4) threads and set the blockDim in the kernel function to a fixed two-dimensional block size of 32*32;

[0093] Each thread has a unique row and column number, with row numbers ranging from 0 to N and column numbers ranging from 0 to N-4, which is consistent with the matrix Q2B -1 The size corresponds to

[0094] Each thread calculates the matrix Q2B -1 An element of the thread number (m, n) thread calculation matrix Q2B -1 The element at row m and column n,

[0095] The method is to multiply the mth row of Q2 by B -1 Get the element at the corresponding position in the nth column;

[0096] Allocate N×N threads again and set the blockDim in the kernel function to a fixed two-dimensional block size of 32*32;

[0097] Each thread has a unique row and column number, with row numbers ranging from 0 to N and column numbers ranging from 0 to N, which is similar to the matrix The size corresponds to

[0098] Each thread calculates the matrix An element of the thread calculation matrix of thread number (m, n) The element in the mth row and nth column, where the principle is Q2B -1 Multiply the mth row of The nth column of gets the element at the corresponding position, where The mth row of is the mth column of Q2. The row index is converted to Q2 column index, thereby reducing one input; finally, the matrix D is obtained.

[0099] (3) Using the traditional method on the CPU to complete the multiplication of D and Y to obtain the non-affine component W;

[0100] (4) Using the traditional method on the CPU, the calculation of A is completed to obtain the affine component, and thus all the deformation parameters of the smooth thin plate spline are obtained.

[0101] As shown in Table 1, it is a schematic table of the running time results of the smooth thin plate spline deformation based on GPU acceleration under 1742 sets of control points.

[0102] Table 1 Schematic diagram of the running time results of the smooth thin plate spline deformation based on GPU acceleration under 1742 control points

[0103] Points: 1742 groups First deformation time Second deformation time The third deformation time Average deformation time Not accelerated 25.83s 26.26s 26.18s 26.09s After acceleration 3.01s 2.96s 3.08s 3.01s

[0104] As can be seen from Table 1, based on the platform with CPU i7-10700f and GPU RTX-2080, 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 9 times.

[0105] The above shows and describes the basic principles, main features and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The above embodiments and descriptions only describe the principles of the present invention. The present invention may be subject to various changes and improvements without departing from the spirit and scope of the present invention. These changes and improvements fall within the scope of the present invention. The scope of protection claimed by the present invention is defined by the attached claims and their equivalents.

Claims

1. A method for calculating deformation parameters of smooth thin plate splines for image registration based on CUDA acceleration, characterized in that: The following steps are involved: 11) Processing of three-dimensional feature points to be registered and template feature points: Obtain N groups of feature points to be registered and template feature points as input, calculate the distance between each feature point to be registered and each template feature point, and put them into a matrix U of size N×N, put the three-dimensional coordinates of the feature points to be registered into a matrix X of size N×N, and put the three-dimensional coordinates of the template feature points into a matrix Y of size N×4; 12) Perform QR decomposition on the matrix of the points to be registered: Use the functions and orthogonal extension kernel functions in the CUDA-based cusolver to perform QR decomposition on the matrix X to obtain the Q matrix and the R matrix; 13) Calculation of deformation parameters of smooth thin plate spline: Using matrices U and Y and matrices Q and R obtained in the previous step, a multiplication kernel function suitable for calculation is designed in CUDA, and based on the functions in CUDA's cublas, the affine component and non-affine component of the deformation of the smooth thin plate spline are calculated respectively to obtain the deformation parameters of the smooth thin plate spline; The calculation of the deformation parameters of the smooth thin plate spline comprises the following steps: 131) The calculation of affine component and non-affine component is completed by using the following formula and CUDA-based cublas library function and CUDA-based matrix multiplication kernel function. The calculation formula is as follows: in, B=C+λI, W of size N×4 is the non-affine component, λ is the smoothing coefficient, I of size (N-4)×(N-4) is the identity matrix, and A of size 4×4 is the affine component; 132) Calculate W: 1321)Completed using CUDA The multiplication of U and Q2 results in a matrix C of size (N-4)×(N-4); 1322) Add the result of the previous step calculation to the I matrix multiplied by the smoothing coefficient to obtain a matrix B with a size of (N-4)×(N-4), and use the function in the cusolver library based on CUDA to obtain the inverse matrix of B; 1323) Use CUDA to complete Q2 and B -1 , The product of gets the matrix D of size N×N, First, transform the input matrices Q2 and B -1 As a one-dimensional array, allocate N×(N-4) threads and set the blockDim in the kernel function to a fixed two-dimensional block size of 32*32; Each thread has a unique row and column number, with row numbers ranging from 0 to N and column numbers ranging from 0 to N-4, which is consistent with the matrix Q2B -1 The size corresponds to Each thread calculates the matrix Q2B -1 An element of the thread number (m, n) thread calculation matrix Q2B -1 The element at row m and column n, The method is to multiply the mth row of Q2 by B -1 Get the element at the corresponding position in the nth column; Allocate N×N threads again and set the blockDim in the kernel function to a fixed two-dimensional block size of 32*32; Each thread has a unique row and column number, with row numbers ranging from 0 to N and column numbers ranging from 0 to N, which is similar to the matrix The size corresponds to Each thread calculates the matrix An element of the thread calculation matrix of thread number (m, n) The element in the mth row and nth column, using the Q2B method -1 Multiply the mth row of The nth column of gets the element at the corresponding position, where The mth row of is the mth column of Q2. The row index is converted to Q2 column index, thereby reducing one input; finally, the matrix D is obtained; 133) Perform multiplication of D and Y on the CPU to obtain the non-affine component W; 134) The calculation of A is completed on the CPU to obtain the affine component, and thus all the deformation parameters of the smooth thin plate spline are obtained.

2. The method for calculating deformation parameters of smooth thin plate splines for image registration based on CUDA acceleration according to claim 1, characterized in that: The QR decomposition of the registration point matrix comprises the following steps: 21) Based on the QR decomposition function in the CUDA cusolver library, the matrix X is decomposed by QR to obtain the Q matrix and R matrix, which are expressed as follows: At this time, the Q matrix has only the first four columns with values ​​and is orthogonal, and the remaining columns are 0. The R matrix has only the first four rows with values, and its size is N×4. Q1 is the first four columns of Q, and its size is N×4. Q2 is the extended part, and its size is N×(N-4). R1 ​​is the first four columns and the first four rows of R, and its size is 4×4. 22) Based on the extended orthogonal kernel function designed by CUDA, the extended orthogonal matrix Q is completed, starting from the fifth column; when filling the i-th column, the following steps are adopted, i ≥ 5: 221) Find the sum of each row of Q and then get its square root to get the row number k of the row with the smallest value; 222) Pass Q into the GPU and allocate N threads corresponding to each element of the i-th column; The calculation content of the nth thread is: the nth row of the first column to the i-1th column of Q multiplied by the first to i-1th elements of the kth row; After all threads are executed, the initial value of the expanded i-th column is obtained, and the obtained result is taken out from the GPU and put into a matrix T of size N×1; 223) Add 1 to the kth element of T, normalize T, and fill it into the i-th column of Q; 23) Calculate the orthogonal Q matrix.

3. The method for calculating deformation parameters of smooth thin plate splines for image registration based on CUDA acceleration according to claim 1, characterized in that: The calculation of matrix C includes the following steps: 31) Convert the input matrices U and Q2 into one-dimensional arrays, allocate (N-4)×N threads in the GPU, and set the blockDim in the kernel function to a fixed two-dimensional block size of 32*32; 32) Set each thread to have a unique row and column number, with the row number ranging from 0 to N-4 and the column number ranging from 0 to N, which is consistent with the matrix The size corresponds to 33) Set the thread calculation matrix of thread number (m, n) The element at row m and column n, The method is The mth row of U is multiplied by the nth column of U, where The mth row of is the mth column of Q2. The row index is converted to Q2 column index, thus reducing one input; 34) Allocate (N-4)×(N-4) threads again and set the blockDim in the kernel function to a fixed two-dimensional block size of 32*32; 35) Set each thread to have a unique row and column number, with the row number ranging from 0 to N-4 and the column number ranging from 0 to N-4, which is consistent with the matrix The size corresponds to 36) Set each thread to calculate the matrix An element of the thread calculation matrix of thread number (m, n) The element at row m and column n, The method is Multiply the mth row of U by the nth column of U to get the element at the corresponding position; 37) Finally, we get the matrix C.

Images