A 3D biological image morphing method based on joint local and global constraints

Abstract

The present invention relates to a three-dimensional biological image deformation method based on local and global joint constraints, comprising: obtaining a global matching point set; calculating local constraints based on label index information ind; calculating global constraints; calculating a downsampled displacement field DF based on the local and global constraints; calculating a cubic B-spline basis function matrix; obtaining a final generated displacement field DF_whole; and performing three-dimensional biological image deformation to obtain a deformed three-dimensional biological image. The present invention utilizes a biological image anatomical atlas to calculate the downsampled displacement field with local and global joint constraints, and interpolates and smoothes the downsampled displacement field blocks by sliding the influence area of ​​the cubic B-spline basis function matrix to obtain the final displacement field and perform three-dimensional biological image deformation. The present invention generates a smooth displacement field that links the local and global components, which can effectively improve the accuracy of local deformation during biological image deformation and registration, while also achieving a smooth overall deformation effect.

Description

Technical Field

[0001] The present invention relates to the technical field of image deformation and image registration, in particular to a three-dimensional biological image deformation method based on local and global joint constraints. Background Art

[0002] In the field of image deformation and image registration, a typical task is to deform the morphological positions of one image into the morphological positions of another image, or to deform and register one image with a reference image. Thin Plate Spline Interpolation (TPS) is a commonly used method. The basic principle of TPS deformation is to achieve overall image deformation by defining a set of control points on the image and deforming these control points. During TPS deformation, the deformation displacement of each pixel relative to the control points needs to be calculated. Then, based on these deformation displacements, the image can be interpolated to generate the final deformed image.

[0003] Based on TPS deformation, the TPS deformation effect can be made smoother by adding a constraint scheme, namely STPS. STPS deformation can increase the accuracy of deformation by adding more control points, but when it uses the overall control points to add smoothness constraints to STPS deformation, the control points in the local area will be affected by the smoothness constraints of the control points in other areas, resulting in errors in the calculated displacement field. Moreover, larger 3D images require more control point sets to calculate the displacement field, which will result in longer deformation time. Downsampling the deformation displacement field is a typical solution, which also leads to the interpolation problem of the downsampling displacement field. To improve the performance and efficiency of 3D image deformation, Qu et al. designed an ultrafast image warping tool (LittleQuickWarp). Their approach involves generating a set of matching points between the biological image to be registered and the target biological image. The downsampled displacement field is then calculated using the TPS algorithm. Trilinear or B-spline interpolation is then performed on the downsampled displacement field blocks. To achieve overall smoothness, the interpolation process uses locally missing displacement field blocks to interpolate a large number of displacement field values ​​from other blocks and then smoothes them. This approach achieves both improved deformation speed and overall smoothness. However, the smoothing constraints between the control point sets do not effectively utilize the local anatomical information of the biomedical image. Using the global control point set to calculate the constraints can lead to errors in the displacement field of local biological tissue in the biomedical image. Furthermore, the downsampled displacement field interpolation scheme employed by this method introduces a large number of displacement field values ​​from other blocks into the block displacement field, which can also introduce deformation errors. Summary of the Invention

[0004] In order to solve the problem of low accuracy of local deformation effects during three-dimensional biological image deformation, the purpose of the present invention is to provide a three-dimensional biological image deformation method based on local and global joint constraints, which can effectively improve the accuracy of local deformation during biological image deformation and alignment, while achieving smooth overall deformation effects.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: a three-dimensional biological image deformation method based on local and global joint constraints, the method comprising the following steps in sequence:

[0006] (1) Obtaining the overall matching point set: obtaining the control point set tar_markers based on the anatomical atlas of the target biological image, copying the control point set tar_markers and mapping it to the biological image to be deformed, obtaining the point set sub_markers of the image to be deformed, and the control point set tar_markers and the point set sub_markers of the image to be deformed form the overall matching point set;

[0007] (2) Calculate local constraints based on the label index information ind: extract the local matching point set corresponding to the label index information ind, and calculate the affine item parameter matrix D of the STPS deformation of the area corresponding to the label index information ind ind And the non-affine parameter matrix C ind ;

[0008] (3) Calculate global constraints: Use the overall matching point set to calculate the affine term parameter matrix D of the STPS deformation all And the non-affine parameter matrix C all ;

[0009] (4) Calculate the downsampled displacement field DF based on local constraints and global constraints: According to the anatomical atlas of the target biological image, use the corresponding local constraints to calculate the displacement field at the local voxel position; for the background area, use the global constraints to calculate the displacement field at the background voxel position;

[0010] (5) Calculate the cubic B-spline basis function matrix;

[0011] (6) Traverse the downsampled displacement field DF, and after preliminary interpolation of the displacement field block, perform smoothing calculation by moving the influence area of ​​the cubic B-spline basis function matrix and the displacement field block after preliminary interpolation to obtain the displacement field of the target image size after final interpolation, that is, the final generated displacement field DF_whole;

[0012] (7) The three-dimensional biological image is deformed according to the finally generated displacement field DF_whole to obtain the deformed three-dimensional biological image.

[0013] The step (1) specifically includes the following steps:

[0014] (1a) At the boundary of each local area of ​​the anatomical atlas of the target biological image, a set of point sets is generated using the PCL point cloud library, and the label index information ind of the corresponding local area is added. This set of point sets is recorded as the control point set tar_markers;

[0015] (1b) Copy the control point set tar_markers and map it to the biological image to be deformed. Adjust the mapped point set so that the mapped point set in the biological image to be deformed and the control point set in the target biological image can be accurately matched at corresponding positions. The accurately matched point set is recorded as the point set sub_markers of the image to be deformed.

[0016] (1c) The control point set tar_markers and the image point set sub_markers to be deformed form an overall matching point set; the size of the control point set tar_markers and the image point set sub_markers to be deformed are both nCpt×3, where the number of rows nCpt is the number of matching point pairs, and the three columns store the x, y, and z coordinates of each point respectively.

[0017] The step (2) specifically includes the following steps:

[0018] (2a) According to the label index information ind, the matching point set corresponding to the label index information ind area is extracted from the control point set tar_markers and the deformed image point set sub_markers and stored in the local control point set cpt_target respectively. ind and the local image point set cpt_subject to be deformed ind , local control point set cpt_target ind and the local image point set cpt_subject to be deformed ind The dimensions are all nCpt ind ×3, where the number of rows is nCpt ind The number of matching point pairs in the corresponding region of the label index information ind, and the three columns store the x, y, and z coordinates of each point respectively;

[0019] (2b) Define the cpt_target for storing the local control point set ind The matrix of the Euclidean distance between each two points is xnxn_K ind ,xnxn_K ind The size is nCpt ind ×nCpt ind ; Calculate the local control point set cpt_target ind The Euclidean distance d between each two points local_ij , and store the inverted Euclidean distance result into the matrix xnxn_Kind , calculated as follows:

[0020] d local_ij_x =cpt_target ind [i][0]-cpt_target ind [j][0];

[0021] d local_ij_y =cpt_target ind [i][1]-cpt_target ind [j][1];

[0022] d local_ij_z =cpt_target ind [i][2]-cpt_target ind [j][2];

[0023]

[0024] Among them, cpt_target ind [i][0] is the local control point set cpt_target ind The x coordinate of the i-th point, cpt_target ind [i][1] is the local control point set cpt_target ind The y coordinate of the i-th point, cpt_target ind [i][2] is the local control point set cpt_target ind The z coordinate of the i-th point in cpt_target ind [j][0] is the local control point set cpt_target ind The x-coordinate of the j-th point in cpt_target ind [j][1] is the local control point set cpt_target ind The y coordinate of the jth point in cpt_target ind [j][2] is the local control point set cpt_target ind The z coordinate of the jth point in d local_ij_x cpt_target ind [i][0] and cpt_target ind The distance between [j][0], d local_ij_y cpt_target ind [i][1] and cpt_target ind The distance between [j][1], d local_ij_z cpt_target ind[i][2] and cpt_target ind The distance between [j][2];

[0025] (2c) Allocate two nCpts ind ×4 matrix X ind With Y ind , respectively used to store the local control point set cpt_target ind and the local image point set cpt_subject to be deformed ind Position information, where the number of rows is nCpt ind The number of matching point pairs in the corresponding area of ​​the label index information ind, X ind The first column of the four columns is set to 1, and the remaining three columns store the local control point set cpt_target ind The x, y, and z coordinates of each point; ind The first column of the four columns is set to 1, and the remaining three columns store the local image point set cpt_subject to be deformed ind The x, y, and z coordinates of each point;

[0026] (2d) Allocate a size of nCpt ind ×nCpt ind The matrix Q ind , matrix Q ind The first 4 columns and matrix X ind The same assignment; use QR decomposition to transform the matrix Q ind Further decomposed into an orthogonal matrix and an upper triangular matrix R ind , and reassign the obtained orthogonal matrix to the matrix Q ind , use the Gram-Schmidt orthogonalization process to obtain the matrix Q ind Expanded to a matrix Q with orthogonal columns ind ;

[0027] (2e) Allocation matrix q1 ind ,q2 ind 、r ind , the obtained Q ind Columns 1 to 4 are stored in matrix q1 ind , Q ind 5th to nth Cpt ind Columns are stored in matrix q2 ind , the obtained matrix R ind The 1st to 4th rows and 1st to 4th columns of ind ;

[0028] (2f) Define the STPS deformation non-affine parameter matrix of the region corresponding to the label index information ind as C ind, define the STPS deformation affine term parameter matrix of the area corresponding to the label index information ind as D ind , define C ind The size is nCpt ind ×4, D ind The size is 4×4, C ind With D ind The calculation process is as follows:

[0029] A local =q2 ind T ·xnxn_K ind q2 ind +I 0.2

[0030] C ind =q2 ind ·(A local -1 q2 ind T ·Y ind )

[0031] D ind =r ind -1 ·q1 ind T ·(Y ind -xnxn_K ind ·C ind )

[0032] Where I is (nCpt ind -4)-order identity matrix; A local is the kernel matrix of local TPS deformation;

[0033] (2g) Repeat steps (2a) to (2f) using the label index information ind until the local constraints of all partitions are solved.

[0034] The step (3) specifically includes the following steps:

[0035] (3a) Define the matrix xnxn_K for storing the Euclidean distance between each two points in the control point set tar_markers, and the size of xnxn_K is nCpt×nCpt; calculate the Euclidean distance d between each two points in the control point set tar_markers global_ij , and store the inverted Euclidean distance result into xnxn_K. The calculation method is as follows:

[0036] d global_ij_x =tar_markers[i][0]-tar_markers[j][0];

[0037] d global_ij_y =tar_markers[i][1]-tar_markers[j][1];

[0038] d global_ij_z =tar_markers[i][2]-tar_markers[j][2];

[0039]

[0040] Among them, tar_markers[i][0] is the x coordinate of the i-th point in the control point set tar_markers, tar_markers[i][1] is the y coordinate of the i-th point in the control point set tar_markers, tar_markers[i][2] is the z coordinate of the i-th point in the control point set tar_markers, tar_markers[j][0] is the x coordinate of the j-th point in the control point set tar_markers, tar_markers[j][1] is the y coordinate of the j-th point in the control point set tar_markers, tar_markers[j][2] is the z coordinate of the j-th point in the control point set tar_markers, d global_ij_x is the distance between tar_markers[0] and tar_markers[j][0], d global_ij_y is the distance between tar_markers[i][1] and tar_markers[j][1], d global_ij_z is the distance between tar_markers[i][2] and tar_markers[j][2];

[0041] (3b) Define two matrices X and Y of size nCpt × 4, which are used to store the position information of the control point set tar_markers and the image point set sub_markers to be deformed, respectively. The number of rows nCpt is the number of matching point pairs. The first column of the four columns in X is set to 1, and the remaining three columns store the x, y, and z coordinates of each point in tar_markers respectively; the first column of the four columns in Y is set to 1, and the remaining three columns store the x, y, and z coordinates of each point in sub_markers respectively.

[0042] (3c) Define a matrix Q of size nCpt × nCpt, with the first four columns of Q assigned the same values ​​as X. Use QR decomposition to further decompose the matrix Q into an orthogonal matrix and an upper triangular matrix R, and reassign the obtained orthogonal matrix to Q. Use the Gram-Schmidt orthogonalization process to expand the obtained matrix Q into a matrix Q with orthogonal columns.

[0043] (3d) Define matrices q1, q2, and r. Store the first to fourth columns of the obtained Q into matrix q1, store the fifth to nCptth columns of Q into matrix q2, and store the first to fourth rows and first to fourth columns of the obtained matrix R into r.

[0044] (3e) Define the parameter matrix of the non-affine term of the overall STPS deformation as C all , the affine term parameter matrix is ​​D all , define C all The size is nCpt×4, D all The size is 4×4; C all With D all The calculation process is as follows:

[0045] A global =q2 T ·xnxn_K·q2+I a 0.2

[0046] C all =q2·(A global -1 q2 T ·Y)

[0047] D all =r -1 ·q1 T ·(Y-xnxn_K·C all )

[0048] Among them, I a is the identity matrix of order (nCpt-4), A global is the kernel matrix of the global TPS deformation.

[0049] The step (4) specifically includes the following steps:

[0050] (4a) Allocate a downsampled displacement field DF, whose size is determined by the downsampled size of 4 times the target image size, and each displacement field element contains the displacement field information in the x, y, and z directions;

[0051] (4b) Traverse the spatial position of the downsampled displacement field DF, and record the traversal index as (x dfi ,y dfi ,z dfi ), multiply the traversal index of the downsampled displacement field by the downsampling multiple to get the actual current voxel position (x ori ,y ori ,z ori ):

[0052] x ori =x dfi ×4

[0053] y ori =y dfi ×4

[0054] z ori =z dfi ×4

[0055] Allocate a matrix X to store the current voxel position information ori , whose size is 1×4, where the first column is set to 1, and the second to fourth columns store x respectively. ori 、y ori 、z ori ;

[0056] (4c) Select the current voxel position for calculating the current control point set cpt_target_ref and the corresponding constraints, which include the affine term parameter D_ref and the non-affine term parameter C_ref: Determine the area to which the current voxel position belongs on the target image anatomical atlas. If the current voxel position belongs to the area corresponding to the label index information ind, then select the local control point set cpt_target ind As cpt_target_ref, select D ind As D_ref, select C ind As C_ref, to calculate the corresponding position displacement field of the area; if the current voxel position belongs to the background area, then select the control point set tar_markers as cpt_target_ref, select the affine item parameter matrix D all As D_ref, select the non-affine parameter matrix C all As C_ref, to calculate the corresponding position displacement field of the background area;

[0057] (4d) Allocate the matrix xmxn_K used to calculate the Euclidean distance between the current voxel position and the cpt_target_ref point set, with a size of 1×nCpt_ref, where nCpt_ref is the number of points in the current control point set cpt_target_ref. The Euclidean distance calculation process between the current voxel position and the i-th point is as follows:

[0058] d ppx_i =x ori -cpt_target_ref[i][0];

[0059] d ppy_i =y ori -cpt_target_ref[i][1];

[0060] d ppz_i =z ori-cpt_target_ref[i][2];

[0061]

[0062] Among them, cpt_target_ref[i][0] is the x coordinate of the i-th point in the current control point set cpt_target_ref, cpt_target_ref[i][1] is the y coordinate of the i-th point in the current control point set cpt_target_ref, cpt_target_ref[i][2] is the z coordinate of the i-th point in the current control point set cpt_target_ref, d ppx_i is the distance between the x coordinate of the current voxel position and cpt_target_ref[i][0], d ppy_i is the distance between the y coordinate of the current voxel position and cpt_target_ref[i][1], d ppz_i The distance between the z coordinate of the current voxel position and cpt_target_ref[i][2];

[0063] (4e) Allocate the matrix x of the stps displacement field used to calculate the downsampled displacement field DF stps , whose size is 1×4, and uses the corresponding constraints D_ref and C_ref at the current voxel position to calculate the stps displacement field matrix of the current voxel position:

[0064] x stps =X ori ·D_ref+xmxn_K·C_ref

[0065] (4f) The matrix x of the stps displacement field calculated at the current voxel position stps Arrange to the current traversal index (x) of the downsampled displacement field DF dfi ,y dfi ,z dfi )

[0066] DF[x dfi ][y dfi ][z dfi ].sx=x stps (1,2)-x ori

[0067] DF[x dfi ][y dfi ][z dfi ].sy=x stps (1,3)-y ori

[0068] DF[x dfi ][ydfi ][z dfi ].sz=x stps (1,4)-z ori

[0069] Where DF[x dif ][y dfi ][z dfi ].sx is the downsampled displacement field DF in (x dfi ,y dfi ,z dfi ) The displacement field value in the x direction at the index, and the same applies to the y and z directions;

[0070] (4g) By traversing the spatial position of the downsampled displacement field DF, repeating the above steps (4a) to (4f), the final downsampled displacement field DF based on local constraints and global constraints is obtained.

[0071] The step (5) specifically includes the following steps:

[0072] (5a) Assign a 4×4 matrix B as the coefficient matrix of the cubic B-spline basis function, whose values ​​are as follows:

[0073]

[0074] (5b) Allocate a 4×4 matrix T to store the data of the cubic B-spline basis function:

[0075]

[0076] Where t1 = 0,

[0077] (5c) Calculate the preliminary cubic B-spline basis function matrix TB of size 4×4:

[0078] TB=T·B

[0079] (5d) Calculate the Kronecker product of TB and TB, and denote the resulting matrix as BXB, whose size is 4 2 ×4 2 :

[0080]

[0081] (5e) According to the calculation principle of Kronecker product in step (5d), perform Kronecker product operation on matrix BXB and matrix TB to obtain a size of 4 3 ×4 3 The cubic B-spline basis function matrix X for three-dimensional interpolation bspline .

[0082] The step (6) specifically includes the following steps:

[0083] (6a) Traverse the downsampled displacement field DF and traverse the index position (x dfi ,y dfi ,z dfi ) corresponds to the actual voxel position (x dfi *4,y dfi *4,z dfi *4) is set as startp;

[0084] (6b) Starting from startp, find 8 voxel positions with downsampled displacement fields DF. The index values ​​of the displacement fields at these 8 voxel positions in the downsampled displacement field DF are: (x dfi ,y dfi ,z dfi )、(x dfi +1,y dfi ,z dfi )、(x dfi ,y dfi +1,z dfi )、(x dfi ,y dfi ,z dfi +1)、(x dfi +1,y dfi +1,z dfi )、(x dfi +1,y dfi ,z dfi +1)、(x dfi ,y dfi +1,z dfi +1)、(x dfi +1,y dfi +1,z dfi +1); Since the size of the downsampled displacement field DF is four times the size of the target image, each of the eight voxel positions found above has a 4×4×4 voxel space behind it that lacks the displacement field, and together they form a 8×8×8 voxel space that lacks the displacement field;

[0085] (6c) For the eight 4×4×4 voxel spaces lacking displacement fields found in step (6b), preliminary displacement field interpolation is performed on each of the 4×4×4 voxel spaces lacking displacement fields according to the position of the voxel in the target anatomical atlas: for each 4×4×4 voxel space lacking displacement fields, only the eight vertices in the space have downsampled displacement field DF values. First, the regions to which the eight vertices belong are determined according to the target anatomical atlas, and then the voxel in the 4×4×4 voxel space is determined to be in the same region as the vertex according to the target anatomical atlas. The displacement field value lacking at the corresponding voxel position is interpolated to the displacement field value at the vertex in the same region as the voxel; if the voxel position is in a different region from the eight vertices, the displacement field value at the voxel position is interpolated to the displacement field value at the vertex closest to the voxel position.

[0086] (6d) Traverse the 4×4×4 voxel space starting from startp, and move the cubic B-spline basis function matrix X bspline The action area of ​​​​the displacement field block after the preliminary interpolation in this space is smoothed. The specific process is: the current voxel position of the traversal is set as currentp, and a displacement field block with a size of 4×4×4 that has been preliminarily interpolated in step (6c) is found backward from currentp. This displacement field block is named DF_smallblock, and DF_smallblock is sorted into a matrix DF_smallblock_M with a size of 64×3 according to the displacement field value. The 64 rows represent the 64 displacement fields in the 4×4×4 displacement field block, and the 3 columns represent the displacement field components of each position in the x, y, and z directions respectively; then the cubic B-spline basis function matrix X is used bspline Perform smoothing calculation on DF_smallblock_M and store the result into matrix DF_B_block:

[0087] DF_B_block=X bspline DF_smallblock_M

[0088] The first row of the DF_B_block result is used as the final interpolated displacement field at the currentp voxel position, where the values ​​of the three columns of the first row are respectively used as the x-, y-, and z-direction components of the displacement field at the currentp voxel position;

[0089] (6e) Repeat steps (6a) to (6d) until the downsampled displacement field DF is traversed, while controlling the value not to exceed the boundary, and the final generated displacement field DF_whole will be obtained after the complete cubic B-spline basis function interpolation and local and global joint constraints.

[0090] The step (7) specifically includes the following steps:

[0091] (7a) Allocate a blank space img_warp of the target image size to store the result of image deformation;

[0092] (7b) Traverse the voxel position in the blank space img_warp, use the x, y, z coordinates of the voxel position pos_tar plus the x, y, z components of the displacement field at the corresponding position of the final generated displacement field DF_whole to obtain the position pos_sub, and assign the voxel at the position pos_sub in the image to be deformed to the position pos_tar in img_warp through trilinear interpolation;

[0093] (7c) After the operation of step (7b), the three-dimensional biological image img_warp is finally obtained after deformation based on local constraints and global constraints.

[0094] It can be seen from the above technical solution that the beneficial effects of the present invention are: first, the present invention uses the anatomical atlas of biological images to calculate the downsampled displacement field with local and global joint constraints, and interpolates and smoothes the downsampled displacement field blocks by sliding the influence area of ​​the cubic B-spline basis function matrix to obtain the final displacement field and perform three-dimensional biological image deformation; second, the present invention generates a smooth displacement field that is locally and globally connected, which can effectively improve the accuracy of local deformation in the process of biological image deformation and alignment, while achieving a smooth overall deformation effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0095] Figure 1 is a flow chart of the method of the present invention;

[0096] Figure 2 This is a schematic diagram of the process structure when smoothing calculation is performed on the displacement field block after interpolation and the influence area of ​​the moving cubic B-spline basis function matrix in the present invention;

[0097] Figure 3 This figure shows the application effect of the present invention in the cross-modal registration of the 3D monkey brain, and the comparison with the effect of the overall STPS deformation scheme. DETAILED DESCRIPTION

[0098] like Figure 1 As shown, a three-dimensional biological image deformation method based on local and global joint constraints includes the following steps in sequence:

[0099] (1) Obtaining the overall matching point set: obtaining the control point set tar_markers based on the anatomical atlas of the target biological image, copying the control point set tar_markers and mapping it to the biological image to be deformed, obtaining the point set sub_markers of the image to be deformed, and the control point set tar_markers and the point set sub_markers of the image to be deformed form the overall matching point set;

[0100] (2) Calculate local constraints based on the label index information ind: extract the local matching point set corresponding to the label index information ind, and calculate the affine item parameter matrix D of the STPS deformation of the area corresponding to the label index information ind ind And the non-affine parameter matrix C ind ;

[0101] (3) Calculate global constraints: Use the overall matching point set to calculate the affine term parameter matrix D of the STPS deformation all And the non-affine parameter matrix C all ;

[0102] (4) Calculate the downsampled displacement field DF based on local constraints and global constraints: According to the anatomical atlas of the target biological image, use the corresponding local constraints to calculate the displacement field at the local voxel position; for the background area, use the global constraints to calculate the displacement field at the background voxel position;

[0103] (5) Calculate the cubic B-spline basis function matrix;

[0104] (6) Traverse the downsampled displacement field DF, and after preliminary interpolation of the displacement field block, perform smoothing calculation by moving the influence area of ​​the cubic B-spline basis function matrix and the displacement field block after preliminary interpolation to obtain the displacement field of the target image size after final interpolation, that is, the final generated displacement field DF_whole;

[0105] (7) The three-dimensional biological image is deformed according to the finally generated displacement field DF_whole to obtain the deformed three-dimensional biological image.

[0106] The step (1) specifically includes the following steps:

[0107] (1a) At the boundary of each local area of ​​the anatomical atlas of the target biological image, a set of point sets is generated using the PCL point cloud library, and the label index information ind of the corresponding local area is added. This set of point sets is recorded as the control point set tar_markers;

[0108] (1b) Copy the control point set tar_markers and map it to the biological image to be deformed. Adjust the mapped point set so that the mapped point set in the biological image to be deformed and the control point set in the target biological image can be accurately matched at corresponding positions. The accurately matched point set is recorded as the point set sub_markers of the image to be deformed.

[0109] (1c) The control point set tar_markers and the image point set sub_markers to be deformed form an overall matching point set; the size of the control point set tar_markers and the image point set sub_markers to be deformed are both nCpt×3, where the number of rows nCpt is the number of matching point pairs, and the three columns store the x, y, and z coordinates of each point respectively.

[0110] The step (2) specifically includes the following steps:

[0111] (2a) According to the label index information ind, the matching point set corresponding to the label index information ind area is extracted from the control point set tar_markers and the deformed image point set sub_markers and stored in the local control point set cpt_target respectively. ind and the local image point set cpt_subject to be deformed ind , local control point set cpt_target ind and the local image point set cpt_subject to be deformed ind The dimensions are all nCpt ind ×3, where the number of rows is nCpt ind The number of matching point pairs in the corresponding region of the label index information ind, and the three columns store the x, y, and z coordinates of each point respectively;

[0112] (2b) Define the cpt_target for storing the local control point set ind The matrix of the Euclidean distance between each two points is xnxn_K ind ,xnxn_K ind The size is nCpt ind ×nCpt ind ; Calculate the local control point set cpt_target ind The Euclidean distance d between each two points local_ij , and store the inverted Euclidean distance result into the matrix xnxn_K ind , calculated as follows:

[0113] d local_ij_x =cpt_target ind [i][0]-cpt_target ind [j][0];

[0114] d local_ij_y =cpt_target ind [i][1]-cpt_target ind [j][1];

[0115] d local_ij_z =cpt_targetind [i][2]-cpt_target ind [j][2];

[0116]

[0117] Among them, cpt_target ind [i][0] is the local control point set cpt_target ind The x coordinate of the i-th point, cpt_target ind [i][1] is the local control point set cpt_target ind The y coordinate of the i-th point, cpt_target ind [i][2] is the local control point set cpt_target ind The z coordinate of the i-th point in cpt_target ind [j][0] is the local control point set cpt_target ind The x-coordinate of the j-th point in cpt_target ind [j][1] is the local control point set cpt_target ind The y coordinate of the jth point in cpt_target ind [j]2] is the local control point set cpt_target ind The z coordinate of the jth point in d local_ij_x cpt_target ind [i][0] and cpt_target ind The distance between [j][0], q local_ij_y cpt_target ind [i][1] and cpt_target ind The distance between [j][1], dlocal _ij_z cpt_target ind [i][2] and cpt_target ind The distance between [j][2];

[0118] (2c) Allocate two nCpts ind ×4 matrix X ind With Y ind , respectively used to store the local control point set cpt_target ind and the local image point set cpt_subject to be deformed ind Position information, where the number of rows is nCpt ind The number of matching point pairs in the corresponding area of ​​the label index information ind, X indThe first column of the four columns is set to 1, and the remaining three columns store the local control point set cpt_target ind The x, y, and z coordinates of each point; ind The first column of the four columns is set to 1, and the remaining three columns store the local image point set cpt_subject to be deformed ind The x, y, and z coordinates of each point;

[0119] (2d) Allocate a size of nCpt ind ×nCpt ind The matrix Q ind , matrix Q ind The first 4 columns and matrix X ind The same assignment; use QR decomposition to transform the matrix Q ind Further decomposed into an orthogonal matrix and an upper triangular matrix R ind , and reassign the obtained orthogonal matrix to the matrix Q ind , use the Gram-Schmidt orthogonalization process to obtain the matrix Q ind Expanded to a matrix Q with orthogonal columns ind ;

[0120] (2e) Allocation matrix q1 ind ,q2 ind 、r ind , the obtained Q ind Columns 1 to 4 are stored in matrix q1 ind , Q ind 5th to nth Cpt ind Columns are stored in matrix q2 ind , the obtained matrix R ind The 1st to 4th rows and 1st to 4th columns of ind ;

[0121] (2f) Define the STPS deformation non-affine parameter matrix of the region corresponding to the label index information ind as C ind , define the STPS deformation affine term parameter matrix of the area corresponding to the label index information ind as D ind , define C ind The size is nCpt ind ×4, D ind The size is 4×4, C ind With D ind The calculation process is as follows:

[0122] A local =q2 ind T ·xnxn_K ind q2 ind +I 0.2

[0123] C ind =q2 ind ·(A local -1 q2 ind T ·Y ind )

[0124] D ind =r ind -1 ·q1 ind T ·(Y ind -xnxn_K ind ·C ind )

[0125] Where I is (nCpt ind -4)-order identity matrix; A local is the kernel matrix of local TPS deformation;

[0126] (2g) Repeat steps (2a) to (2f) using the label index information ind until the local constraints of all partitions are solved.

[0127] The step (3) specifically includes the following steps:

[0128] (3a) Define the matrix xnxn_K for storing the Euclidean distance between each two points in the control point set tar_markers, and the size of xnxn_K is nCpt×nCpt; calculate the Euclidean distance d between each two points in the control point set tar_markers global_ij , and store the inverted Euclidean distance result into xnxn_K. The calculation method is as follows:

[0129] d global_ij_x =tar_markers[i][0]-tar_markers[j][0];

[0130] d global_ij_y =tar_markers[i][1]-tar_markers[j][1];

[0131] d global_ij_z =tar_markers[i][2]-tar_markers[j][2];

[0132]

[0133] Among them, tar_markers[i][0] is the x coordinate of the i-th point in the control point set tar_markers, tar_markers[i][1] is the y coordinate of the i-th point in the control point set tar_markers, tar_markers[i][2] is the z coordinate of the i-th point in the control point set tar_markers, tar_markers[j][0] is the x coordinate of the j-th point in the control point set tar_markers, tar_markers[j][1] is the y coordinate of the j-th point in the control point set tar_markers, tar_markers[j][2] is the z coordinate of the j-th point in the control point set tar_markers, d global_ij_x is the distance between tar_markers[i][0] and tar_markers[j][0], d global_ij_y is the distance between tar_markers[i][1] and tar_markers[j][1], d global_ij_z is the distance between tar_markers[i][2] and tar_markers[j][2];

[0134] (3b) Define two matrices X and Y of size nCpt × 4, which are used to store the position information of the control point set tar_markers and the image point set sub_markers to be deformed, respectively. The number of rows nCpt is the number of matching point pairs. The first column of the four columns in X is set to 1, and the remaining three columns store the x, y, and z coordinates of each point in tar_markers respectively; the first column of the four columns in Y is set to 1, and the remaining three columns store the x, y, and z coordinates of each point in sub_markers respectively.

[0135] (3c) Define a matrix Q of size nCpt × nCpt, with the first four columns of Q assigned the same values ​​as X. Use QR decomposition to further decompose the matrix Q into an orthogonal matrix and an upper triangular matrix R, and reassign the obtained orthogonal matrix to Q. Use the Gram-Schmidt orthogonalization process to expand the obtained matrix Q into a matrix Q with orthogonal columns.

[0136] (3d) Define matrices q1, q2, and r. Store the first to fourth columns of the obtained Q into matrix q1, store the fifth to nCptth columns of Q into matrix q2, and store the first to fourth rows and first to fourth columns of the obtained matrix R into r.

[0137] (3e) Define the parameter matrix of the non-affine term of the overall STPS deformation as C all , the affine term parameter matrix is ​​D all , define C all The size is nCpt×4, Dall The size is 4×4; C all With D all The calculation process is as follows:

[0138] A global =q2 T ·xnxn_K·q2+I a 0.2

[0139] C all =q2·(A global -1 q2 T ·Y)

[0140] D all =r -1 ·q1 T ·(Y-xnxn_K·C all )

[0141] Among them, I a is the identity matrix of order (nCpt-4), A global is the kernel matrix of the global TPS deformation.

[0142] The step (4) specifically includes the following steps:

[0143] (4a) Allocate a downsampled displacement field DF, whose size is determined by the downsampled size of 4 times the target image size, and each displacement field element contains the displacement field information in the x, y, and z directions;

[0144] (4b) Traverse the spatial position of the downsampled displacement field DF, and record the traversal index as (x dfi ,y dfi ,z dfi ), multiply the traversal index of the downsampled displacement field by the downsampling multiple to get the actual current voxel position (x ori ,y ori ,z ori ):

[0145] x ori =x dfi ×4

[0146] y ori =y dfi ×4

[0147] z ori =z dfi ×4

[0148] Allocate a matrix X to store the current voxel position information ori , whose size is 1×4, where the first column is set to 1, and the second to fourth columns store x respectively. ori 、yori 、z ori ;

[0149] (4c) Select the current voxel position for calculating the current control point set cpt_target_ref and the corresponding constraints, which include the affine term parameter D_ref and the non-affine term parameter C_ref: Determine the area to which the current voxel position belongs on the target image anatomical atlas. If the current voxel position belongs to the area corresponding to the label index information ind, then select the local control point set cpt_target ind As cpt_target_ref, select D ind As D_ref, select C ind As C_ref, to calculate the corresponding position displacement field of the area; if the current voxel position belongs to the background area, then select the control point set tar_markers as cpt_target_ref, select the affine item parameter matrix D all As D_ref, select the non-affine parameter matrix C all As C_ref, to calculate the corresponding position displacement field of the background area;

[0150] (4d) Allocate the matrix xmaxn_K used to calculate the Euclidean distance between the current voxel position and the cpt_target_ref point set, with a size of 1×nCpt_ref, where nCpt_ref is the number of points in the current control point set cpt_target_ref. The Euclidean distance calculation process between the current voxel position and the i-th point is as follows:

[0151] d ppx_i =x ori -cpt_target_ref[i][0];

[0152] d ppy_i =y ori -cpt_target_ref[i][1];

[0153] d ppz_i =z ori -cpt_target_ref[i][2];

[0154]

[0155] Among them, cpt_target_ref[i][0] is the x coordinate of the i-th point in the current control point set cpt_target_ref, cpt_target_ref[i][1] is the y coordinate of the i-th point in the current control point set cpt_target_ref, cpt_target_ref[i][2] is the z coordinate of the i-th point in the current control point set cpt_target_ref, d ppx_i is the distance between the x coordinate of the current voxel position and cpt_target_ref[i][0], d ppy_i is the distance between the y coordinate of the current voxel position and cpt_target_ref[i][1], d ppz_i The distance between the z coordinate of the current voxel position and cpt_target_ref[i][2];

[0156] (4e) Allocate the matrix x of the stps displacement field used to calculate the downsampled displacement field DF stps , whose size is 1×4, and uses the corresponding constraints D_ref and C_ref at the current voxel position to calculate the stps displacement field matrix of the current voxel position:

[0157] x stps =X ori ·D_ref+xmxn_K·C_ref

[0158] (4f) The matrix x of the stps displacement field calculated at the current voxel position stps Arrange to the current traversal index (x) of the downsampled displacement field DF dfi ,y dfi ,z dfi )

[0159] DF[x dfi ][y dfi ][z dfi ].sx=x stps (1,2)-x ori

[0160] DF[x dfi ][y dfi ][z dfi ].sy=x stps (1,3)-y ori

[0161] DF[x dfi ][y dfi ][z dfi ].sz=x stps (1,4)-z ori

[0162] Where DF[x dfi ][y dfi ][z dfi ].sx is the downsampled displacement field DF in (x dfi ,y dfi ,z dfi ) The displacement field value in the x direction at the index, and the same applies to the y and z directions;

[0163] (4g) By traversing the spatial position of the downsampled displacement field DF, repeating the above steps (4a) to (4f), the final downsampled displacement field DF based on local constraints and global constraints is obtained.

[0164] The step (5) specifically includes the following steps:

[0165] (5a) Assign a 4×4 matrix B as the coefficient matrix of the cubic B-spline basis function, whose values ​​are as follows:

[0166]

[0167] (5b) Allocate a 4×4 matrix T to store the data of the cubic B-spline basis function:

[0168]

[0169] Where t1 = 0,

[0170] (5c) Calculate the preliminary cubic B-spline basis function matrix TB of size 4×4:

[0171] TB=T·B

[0172] (5d) Calculate the Kronecker product of TB and TB, and denote the resulting matrix as BXB, whose size is 4 2 ×4 2 :

[0173]

[0174] (5e) According to the calculation principle of Kronecker product in step (5d), perform Kronecker product operation on matrix BXB and matrix TB to obtain a size of 4 3 ×4 3 The cubic B-spline basis function matrix X for three-dimensional interpolation bspline .

[0175] The step (6) specifically includes the following steps:

[0176] (6a) Traverse the downsampled displacement field DF and traverse the index position (x dfi ,y dfi ,zdfi ) corresponds to the actual voxel position (x dfi *4,y dfi *4,z dfi *4) is set as startp;

[0177] (6b) Starting from startp, find 8 voxel positions with downsampled displacement fields DF. The index values ​​of the displacement fields at these 8 voxel positions in the downsampled displacement field DF are: (x dfi ,y dfi ,z dfi )、(x dfi +1,y dfi ,z dfi )、(x dfi ,y dfi +1,z dfi )、(x dfi ,y dfi ,z dfi +1)、(x dfi +1,y dfi +1,z dfi )、(x dfi +1,y dfi ,z dfi +1)、(x dfi ,z dfi +1,z dfi +1)、(x dfi +1,y dfi +1,z dfi +1); Since the size of the downsampled displacement field DF is four times the size of the target image, each of the eight voxel positions found above has a 4×4×4 voxel space behind it that lacks the displacement field, and together they form a 8×8×8 voxel space that lacks the displacement field;

[0178] (6c) For the eight 4×4×4 voxel spaces lacking displacement fields found in step (6b), preliminary displacement field interpolation is performed on each of the 4×4×4 voxel spaces lacking displacement fields according to the position of the voxel in the target anatomical atlas: for each 4×4×4 voxel space lacking displacement fields, only the eight vertices in the space have downsampled displacement field DF values. First, the regions to which the eight vertices belong are determined according to the target anatomical atlas, and then the voxel in the 4×4×4 voxel space is determined to be in the same region as the vertex according to the target anatomical atlas. The displacement field value lacking at the corresponding voxel position is interpolated to the displacement field value at the vertex in the same region as the voxel; if the voxel position is in a different region from the eight vertices, the displacement field value at the voxel position is interpolated to the displacement field value at the vertex closest to the voxel position.

[0179] (6d) Traverse the 4×4×4 voxel space starting from startp, and move the cubic B-spline basis function matrix X bspline The action area of ​​​​the displacement field block after the preliminary interpolation in this space is smoothed. The specific process is: the current voxel position of the traversal is set as currentp, and a displacement field block with a size of 4×4×4 that has been preliminarily interpolated in step (6c) is found backward from currentp. This displacement field block is named DF_smallblock, and DF_smallblock is sorted into a matrix DF_smallblock_M with a size of 64×3 according to the displacement field value. The 64 rows represent the 64 displacement fields in the 4×4×4 displacement field block, and the 3 columns represent the displacement field components of each position in the x, y, and z directions respectively; then the cubic B-spline basis function matrix X is used bspline Perform smoothing calculation on DF_smallblock_M and store the result into matrix DF_B_block:

[0180] DF_B_block=X bspline DF_smallblock_M

[0181] The first row of the DF_B_block result is used as the final interpolated displacement field at the currentp voxel position, where the values ​​of the three columns of the first row are respectively used as the x-, y-, and z-direction components of the displacement field at the currentp voxel position;

[0182] (6e) Repeat steps (6a) to (6d) until the downsampled displacement field DF is traversed, while controlling the value not to exceed the boundary, and the final generated displacement field DF_whole will be obtained after the complete cubic B-spline basis function interpolation and local and global joint constraints.

[0183] The step (7) specifically includes the following steps:

[0184] (7a) Allocate a blank space img_warp of the target image size to store the result of image deformation;

[0185] (7b) Traverse the voxel position in the blank space img_warp, use the x, y, z coordinates of the voxel position pos_tar plus the x, y, z components of the displacement field at the corresponding position of the final generated displacement field DF_whole to obtain the position pos_sub, and assign the voxel at the position pos_sub in the image to be deformed to the position pos_tar in img_warp through trilinear interpolation;

[0186] (7c) After the operation of step (7b), the three-dimensional biological image img_warp is finally obtained after deformation based on local constraints and global constraints.

[0187] like Figure 2 As shown, the present invention performs smoothing calculation by moving the influence area of ​​the cubic B-spline basis function matrix and the displacement field block after interpolation as follows: traverse the downsampled displacement field DF, and define the real voxel position corresponding to the traversal index position of the downsampled displacement field DF as startp; with startp as the starting point, find 8 voxel positions with DF downsampled displacement field, and each of the 8 voxel positions found has a 4×4×4 voxel space without displacement field behind it; for each 4×4×4 voxel space without displacement field, perform preliminary displacement field interpolation according to the position of the voxel in the target anatomical atlas; traverse the 4×4×4 voxel space with startp as the starting point, define the current voxel position traversed as currentp, and use currentp as the starting point to find a displacement field block of size 4×4×4 that has undergone preliminary interpolation, and move the action area "BsplinebasisMatrixinfluence" of the cubic B-spline basis function matrix X_bspline. area", performs smooth calculation on the 4×4×4 displacement field block after preliminary interpolation starting from currentp, and fills the initial value of the calculation result into the currentp voxel position until the displacement field in the 4×4×4 voxel space starting from startp is interpolated; according to the traversal index of the downsampled displacement field DF, move the position of startp and repeat the above interpolation steps until the displacement field is interpolated and filled.

[0188] like Figure 3 As shown, the upper left is a slice of the monkey brain image to be deformed, the upper right is the corresponding slice of the cross-modal target monkey brain image, the lower left is the effect diagram of the LittleQuickWarp deformation scheme, and the lower right is the deformation effect diagram of the present invention. By comparison, it can be seen that the deformation scheme proposed in the present invention has a significantly improved effect.

[0189] In summary, the present invention utilizes the anatomical atlas of biological images to perform downsampled displacement field calculation with local and global joint constraints, and performs interpolation and smoothing operations on the downsampled displacement field blocks by sliding the influence area of ​​the cubic B-spline basis function matrix to obtain the final displacement field and perform three-dimensional biological image deformation; the present invention generates a smooth displacement field that is locally and globally connected, which can effectively improve the accuracy of local deformation in the process of biological image deformation and alignment, while achieving a smooth overall deformation effect.

Claims

1. A three-dimensional biological image deformation method based on local and global joint constraints, characterized by: The method comprises the following steps in sequence: (1) Obtaining the overall matching point set: obtaining the control point set tar_markers based on the anatomical atlas of the target biological image, copying the control point set tar_markers and mapping it to the biological image to be deformed, obtaining the point set sub_markers of the image to be deformed, and the control point set tar_markers and the point set sub_markers of the image to be deformed form the overall matching point set; (2) Calculate local constraints based on the label index information ind: extract the local matching point set corresponding to the label index information ind, and calculate the affine item parameter matrix D of the STPS deformation of the area corresponding to the label index information ind ind And the non-affine parameter matrix C ind ; (3) Calculate global constraints: Use the overall matching point set to calculate the affine term parameter matrix D of the STPS deformation all And the non-affine parameter matrix C all ; (4) Calculate the downsampled displacement field DF based on local constraints and global constraints: According to the anatomical atlas of the target biological image, use the corresponding local constraints to calculate the displacement field at the local voxel position; for the background area, use the global constraints to calculate the displacement field at the background voxel position; (5) Calculate the cubic B-spline basis function matrix; (6) Traverse the downsampled displacement field DF, and after preliminary interpolation of the displacement field block, perform smoothing calculation by moving the influence area of ​​the cubic B-spline basis function matrix and the displacement field block after preliminary interpolation to obtain the displacement field of the target image size after final interpolation, that is, the final generated displacement field DF_whole; (7) The three-dimensional biological image is deformed according to the finally generated displacement field DF_whole to obtain the deformed three-dimensional biological image.

2. The 3D biological image deformation method based on local and global joint constraints according to claim 1, characterized in that: The step (1) specifically includes the following steps: (1a) At the boundary of each local area of ​​the anatomical atlas of the target biological image, a set of point sets is generated using the PCL point cloud library, and the label index information ind of the corresponding local area is added. This set of point sets is recorded as the control point set tar_markers; (1b) Copy the control point set tar_markers and map it to the biological image to be deformed. Adjust the mapped point set so that the mapped point set in the biological image to be deformed and the control point set in the target biological image can be accurately matched at corresponding positions. The accurately matched point set is recorded as the point set sub_markers of the image to be deformed. (1c) The control point set tar_markers and the image point set sub_markers to be deformed form an overall matching point set; the size of the control point set tar_markers and the image point set sub_markers to be deformed are both nCpt×3, where the number of rows nCpt is the number of matching point pairs, and the three columns store the x, y, and z coordinates of each point respectively.

3. The method for 3D biological image deformation based on local and global joint constraints according to claim 1, characterized in that: The step (2) specifically includes the following steps: (2a) According to the label index information ind, the matching point set corresponding to the label index information ind area is extracted from the control point set tar_markers and the image point set to be deformed sub_markers and stored in the local control point set cpt_target respectively. ind and the local image point set cpt_subject to be deformed ind , local control point set cpt_target ind and the local image point set cpt_subject to be deformed ind The dimensions are all nCpt ind ×3, where the number of rows is nCpt ind The number of matching point pairs in the corresponding region of the label index information ind, and the three columns store the x, y, and z coordinates of each point respectively; (2b) Define the cpt_target for storing the local control point set ind The matrix of the Euclidean distance between each two points is xnxn_K ind ,xnxn_K ind The size is nCpt ind ×nCpt ind ; Calculate the local control point set cpt_target ind The Euclidean distance d between each two points local_ij , and store the inverted Euclidean distance result into the matrix xnxn_K ind , calculated as follows: d local_ij_x =cpt_target ind [i][0]-cpt_target ind [j][0]; d local_ij_y =cpt_target ind [i][1]-cpttarget ind [j][1]; d local_ij_z =cpt_target ind [i][2]-cpt_target ind [j][2]; Among them, cpt_target ind [i][0] is the local control point set cpt_target ind The x coordinate of the i-th point, cpt_target ind [i][1] is the local control point set cpt_target ind The y coordinate of the i-th point in cpt_target ind [i][2] is the local control point set cpt_target ind The z coordinate of the i-th point in cpt_target ind [j][0] is the local control point set cpt_target ind The x-coordinate of the j-th point in cpt_target ind [j][1] is the local control point set cpt_target ind The y coordinate of the jth point in cpt_target ind [j][2] is the local control point set cpt_target ind The z coordinate of the jth point in d local_ij_x cpt_target ind [i][0] and cpt_target ind The distance between [j][0], d local_ij_y cpt_target ind [i][1] and cpt_target ind The distance between [j][1], d local_ij_z cpt_target ind [i][2] and cpt_target ind The distance between [j][2]; (2c) Allocate two nCpt ind ×4 matrix X ind With Y ind , respectively used to store the local control point set cpt_target ind and the local image point set cpt_subject to be deformed ind Position information, where the number of rows is nCpt ind The number of matching point pairs in the corresponding area of ​​the label index information ind, X ind The first column of the four columns is set to 1, and the remaining three columns store the local control point set cpt_target ind The x, y, and z coordinates of each point; ind The first column of the four columns is set to 1, and the remaining three columns store the local image point set cpt_subject to be deformed ind The x, y, and z coordinates of each point; (2d) Allocate a size of nCpt ind ×nCpt ind The matrix Q ind , matrix Q ind The first 4 columns and matrix X ind The same assignment; use QR decomposition to transform the matrix Q ind Further decomposed into an orthogonal matrix and an upper triangular matrix R ind , and reassign the obtained orthogonal matrix to the matrix Q ind , use the Gram-Schmidt orthogonalization process to obtain the matrix Q ind Expanded to a matrix Q with orthogonal columns ind ; (2e) Allocation matrix q1 ind ,q2 ind 、r ind , the obtained Q ind Columns 1 to 4 are stored in matrix q1 ind , Q ind 5th to nth Cpt ind Columns are stored in matrix q2 ind , the obtained matrix R ind The 1st to 4th rows and 1st to 4th columns of ind ; (2f) Define the STPS deformation non-affine parameter matrix of the region corresponding to the label index information ind as C ind , define the STPS deformation affine term parameter matrix of the area corresponding to the label index information ind as D ind , define C ind The size is nCpt ind ×4, D ind The size is 4×4, C ind With D ind The calculation process is as follows: <h2 style=";text-align:left;direction:ltr">A<h2 style=";text-align:left;direction:ltr"> local <h2 style=";text-align:left;direction:ltr"> =q2<h2 style=";text-align:left;direction:ltr"> ind <h2 style=";text-align:left;direction:ltr"> T <h2 style=";text-align:left;direction:ltr"> ·xnxn_K<h2 style=";text-align:left;direction:ltr"> ind <h2 style=";text-align:left;direction:ltr"> q2<h2 style=";text-align:left;direction:ltr"> ind <h2 style=";text-align:left;direction:ltr"> +I·0.2 C ind =q2 ind ·(A local -1 ·q2 ind T ·Y ind ) D ind =r ind -1 ·q1 ind T ·(Y ind -xnxn_K ind ·C ind ) Where I is (nCpt ind -4)-order identity matrix; A local is the kernel matrix of local TPS deformation; (2g) Repeat steps (2a) to (2f) using the label index information ind until the local constraints of all partitions are solved.

4. The 3D biological image deformation method based on local and global joint constraints according to claim 1, characterized in that: The step (3) specifically includes the following steps: (3a) Define the matrix xnxn_K for storing the Euclidean distance between each two points in the control point set tar_markers, and the size of xnxn_K is nCpt×nCpt; calculate the Euclidean distance d between each two points in the control point set tar_markers global_ij , and store the inverted Euclidean distance result into xnxn_K. The calculation method is as follows: d global_ij_x =tar_markers[i][0]-tar_markers[j][0]; d global_ij_y =tar_markers[i][1]-tar_markers[j][1]; d global_ij_z =tar_markers[i][2]-tar_markers[j][2]; Among them, tar_markers[i][0] is the x coordinate of the i-th point in the control point set tar_markers, tar_markers[i][1] is the y coordinate of the i-th point in the control point set tar_markers, tar_markers[i][2] is the z coordinate of the i-th point in the control point set tar_markers, tar_markers[j][0] is the x coordinate of the j-th point in the control point set tar_markers, tar_markers[i][1] is the y coordinate of the j-th point in the control point set tar_markers, tar_markers[j][2] is the z coordinate of the j-th point in the control point set tar_markers, d global_ij_x is the distance between tar_markers[j][0] and tar_markers[j][0], d global_ij_y is the distance between tar_markers[i][1] and tar_markers[j][1], d global_ij_z is the distance between tar_markers[i][2] and tar_markers[j][2]; (3b) Define two matrices X and Y of size nCpt × 4, which are used to store the position information of the control point set tar_markers and the image point set sub_markers to be deformed, respectively. The number of rows nCpt is the number of matching point pairs. The first column of the four columns in X is set to 1, and the remaining three columns store the x, y, and z coordinates of each point in tar_markers respectively; the first column of the four columns in Y is set to 1, and the remaining three columns store the x, y, and z coordinates of each point in sub_markers respectively. (3c) Define a matrix Q of size nCpt × nCpt, with the first four columns of Q assigned the same values ​​as X. Use QR decomposition to further decompose the matrix Q into an orthogonal matrix and an upper triangular matrix R, and reassign the obtained orthogonal matrix to Q. Use the Gram-Schmidt orthogonalization process to expand the obtained matrix Q into a matrix Q with orthogonal columns. (3d) Define matrices q1, q2, and r. Store the first to fourth columns of the obtained Q into matrix q1, store the fifth to nCptth columns of Q into matrix q2, and store the first to fourth rows and first to fourth columns of the obtained matrix R into r. (3e) Define the parameter matrix of the non-affine term of the overall STPS deformation as C all , the affine term parameter matrix is ​​D all , define C all The size is nCpt×4, D all The size is 4×4; C all With D all The calculation process is as follows: A global =q2 T ·xnxn_K·q2+Ia·0.2 C all =q2·(A global -1 ·q2 T ·Y) D all =r -1 ·q1 T ·(Y-xnxn_K·C all ) Among them, I a is the identity matrix of order (nCpt-4), A global is the kernel matrix of the global TPS deformation.

5. The 3D biological image deformation method based on local and global joint constraints according to claim 1, characterized in that: The step (4) specifically includes the following steps: (4a) Allocate a downsampled displacement field DF, whose size is determined by the downsampled size of 4 times the target image size, and each displacement field element contains the displacement field information in the x, y, and z directions; (4b) Traverse the spatial position of the downsampled displacement field DF, and record the traversal index as (x dfi ,y dfi ,z dfi ), multiply the traversal index of the downsampled displacement field by the downsampling multiple to get the actual current voxel position (x ori ,y ori ,z ori ): x ori =x dfi ×4 and ori =and dfi ×4 With ori =z dfi ×4 Allocate a matrix x to store the current voxel position information ori , whose size is 1×4, where the first column is set to 1, and the second to fourth columns store x respectively. ori 、y ori 、z ori ; (4c) Select the current voxel position for calculating the current control point set cpt_target_ref and the corresponding constraints, which include the affine term parameter D_ref and the non-affine term parameter C_ref: Determine the area to which the current voxel position belongs on the target image anatomical atlas. If the current voxel position belongs to the area corresponding to the label index information ind, then select the local control point set cpt_target ind As cpt_target_ref, select D ind As D_ref, select C ind As C_ref, to calculate the corresponding position displacement field of the area; if the current voxel position belongs to the background area, then select the control point set tar_markers as cpt_target_ref, select the affine item parameter matrix D all As D_ref, select the non-affine parameter matrix C all As C_ref, to calculate the corresponding position displacement field of the background area; (4d) Allocate the matrix xmxn_K used to calculate the Euclidean distance between the current voxel position and the cpt_target_ref point set, with a size of 1×nCpt_ref, where nCpt_ref is the number of points in the current control point set cpt_target_ref. The Euclidean distance calculation process between the current voxel position and the i-th point is as follows: d ppx_i =x ori -cpt_target_ref[i][0]; d ppy_i =y ori -cpt_target_ref[i][1]; d ppz_i =z ori -cpt_target_ref[i][2]; Among them, cpt_target_ref[i][0] is the x coordinate of the i-th point in the current control point set cpt_target_ref, cpt_target_ref[i][1] is the y coordinate of the i-th point in the current control point set cpt_target_ref, cpt_target_ref[i][2] is the z coordinate of the i-th point in the current control point set cpt_target_ref, d ppx_i is the distance between the x coordinate of the current voxel position and cpt_target_ref[i][0], d ppy_i is the distance between the y coordinate of the current voxel position and cpt_target_ref[i][1], d ppz_i The distance between the z coordinate of the current voxel position and cpt_target_ref[i][2]; (4e) Allocate the matrix x of the stps displacement field used to calculate the downsampled displacement field DF stps , whose size is 1×4, and uses the corresponding constraints D_ref and C_ref at the current voxel position to calculate the stps displacement field matrix of the current voxel position: x stps =X ori ·D_ref+xmxn_K·C_ref (4f) The matrix x of the stps displacement field calculated at the current voxel position stps Arrange to the current traversal index (x) of the downsampled displacement field DF dfi ,y dfi ,z dfi ) DF[x dfi ][y dfi ][z dfi ].sx=x stps (1,2)-x ori DF[x dfi ][y dfi ][z dfi ].sy=x stps (1,3)-y ori DF[x dfi ][y dfi ][z dfi ].sz=x stps (1,4)-z ori Where DF[x dfi ][y dfi ][z dfi ].sx is the downsampled displacement field DF in (x dfi ,y dfi ,z dfi ) The displacement field value in the x direction at the index, and the same applies to the y and z directions; (4g) By traversing the spatial position of the downsampled displacement field DF, repeating the above steps (4a) to (4f), the final downsampled displacement field DF based on local constraints and global constraints is obtained.

6. The 3D biological image deformation method based on local and global joint constraints according to claim 1, characterized in that: The step (5) specifically includes the following steps: (5a) Assign a 4×4 matrix B as the coefficient matrix of the cubic B-spline basis function, whose values ​​are as follows: (5b) Allocate a 4×4 matrix T to store the data of the cubic B-spline basis function: Where t1 = 0, (5c) Calculate the preliminary cubic B-spline basis function matrix TB of size 4×4: TB=T·B (5d) Calculate the Kronecker product of TB and TB, and denote the resulting matrix as BXB, whose size is 4 2 ×4 2 : (5e) According to the calculation principle of Kronecker product in step (5d), perform Kronecker product operation on matrix BXB and matrix TB to obtain a size of 4 3 ×4 3 The cubic B-spline basis function matrix X for three-dimensional interpolation bspline .

7. The 3D biological image deformation method based on local and global joint constraints according to claim 1, characterized in that: The step (6) specifically includes the following steps: (6a) Traverse the downsampled displacement field DF and traverse the index position (x dfi ,y dfi ,z dfi ) corresponds to the actual voxel position (x dfi *4,y dfi *4,z dfi *4) Defined as start p; (6b) Starting from start p, find 8 voxel positions with downsampled displacement fields DF. The index values ​​of the displacement fields at these 8 voxel positions in the downsampled displacement field DF are: (x dfi ,y dfi ,z dfi )、(x dfi +1,y dfi ,z dfi )、(x dfi ,y dfi +1,z dfi )、(x dfi ,y dfi ,z dfi +1)、(x dfi +1,y dfi +1,z dfi )、(x dfi +1,y dfi ,z dfi +1)、(x dfi ,y dfi +1,z dfi +1)、(x dfi +1,y dfi +1,z dfi +1); Since the size of the downsampled displacement field DF is four times the size of the target image, each of the eight voxel positions found above has a 4×4×4 voxel space behind it that lacks the displacement field, and together they form a 8×8×8 voxel space that lacks the displacement field; (6c) For the eight 4×4×4 voxel spaces lacking displacement fields found in step (6b), preliminary displacement field interpolation is performed on each of the 4×4×4 voxel spaces lacking displacement fields according to the position of the voxel in the target anatomical atlas: for each 4×4×4 voxel space lacking displacement fields, only the eight vertices in the space have downsampled displacement field DF values. First, the regions to which the eight vertices belong are determined according to the target anatomical atlas, and then the voxel in the 4×4×4 voxel space is determined to be in the same region as the vertex according to the target anatomical atlas. The displacement field value lacking at the corresponding voxel position is interpolated to the displacement field value at the vertex in the same region as the voxel; if the voxel position is in a different region from the eight vertices, the displacement field value at the voxel position is interpolated to the displacement field value at the vertex closest to the voxel position. (6d) Traverse the 4×4×4 voxel space starting from start p, and move the cubic B-spline basis function matrix X bspline The action area of ​​​​the displacement field block after the preliminary interpolation in this space is smoothed. The specific process is: the current voxel position of the traversal is set as current p, and a displacement field block of size 4×4×4 that has been preliminarily interpolated in step (6c) is found backward from current p. This displacement field block is named DF_smallblock, and DF_smallblock is sorted into a matrix DF_smallblock_M of size 64×3 according to the displacement field value, in which 64 rows represent 64 displacement fields in the 4×4×4 displacement field block, and 3 columns represent the displacement field components of each position displacement field in the x, y, and z directions respectively; then the cubic B-spline basis function matrix X is used bspline Perform smoothing calculation on DF_smallblock_M and store the result into matrix DF_B_block: DF_B_block=X bspline ·DF_smallblock_M The first row of the DF_B_block result is used as the final interpolated displacement field at the current p voxel position, where the values ​​of the three columns of the first row are respectively used as the x-, y-, and z-direction components of the displacement field at the current p voxel position; (6e) Repeat steps (6a) to (6d) until the downsampled displacement field DF is traversed, while controlling the value not to exceed the boundary, and the final generated displacement field DF_whole will be obtained after the complete cubic B-spline basis function interpolation and local and global joint constraints.

8. The 3D biological image deformation method based on local and global joint constraints according to claim 1, characterized in that: The step (7) specifically includes the following steps: (7a) Allocate a blank space img_warp of the target image size to store the result of image deformation; (7b) Traverse the voxel position in the blank space img_warp, use the x, y, z coordinates of the voxel position pos_tar plus the x, y, z components of the displacement field at the corresponding position of the final generated displacement field DF_whole to obtain the position pos_sub, and assign the voxel at the position pos_sub in the image to be deformed to the position pos_tar in img_warp through trilinear interpolation; (7c) After the operation of step (7b), the three-dimensional biological image img_warp is finally obtained after deformation based on local constraints and global constraints.

Images