A deformable registration method for 3D multimodal images for rigid targets

Abstract

The present invention discloses a deformable registration method for three-dimensional multimodal images for rigid targets. By constructing a point cloud from a multimodal three-dimensional image, using normal and curvature information to construct a local feature coordinate system, the iterative nearest neighbor algorithm is corrected and the error points are eliminated to improve the robustness of the registration. Based on the structural information of the point cloud itself, partition registration is performed to meet the differences between regions on the basis of local rigid registration, to achieve deformable rigid registration, and to improve the overall registration accuracy. The existing iterative nearest neighbor registration method does not make sufficient use of the local features of the point cloud and has large errors. Therefore, the present invention makes full use of the multiple features of the point cloud to construct a local feature reference coordinate system to improve the registration process, and combines the structural characteristics of the point cloud to partition the point cloud to achieve a better registration effect between the registration point cloud and the reference point cloud.

Description

Technical Field

[0001] The invention relates to a three-dimensional registration method for multi-modal images, in particular to a deformable registration method for three-dimensional multi-modal images for rigid targets. Background Art

[0002] In modern medical diagnosis, it has become a trend to use a variety of medical imaging technologies to obtain three-dimensional images for comprehensive diagnosis. Multimodal three-dimensional registration can combine the advantages of different imaging technologies and improve the accuracy of diagnosis. It has a high practical value in clinical medicine. The three-dimensional images obtained by medical imaging technology are three-dimensional images combined with two-dimensional sequence tomography images. The main multimodal three-dimensional registration methods include registration methods based on feature points, optimization registration methods based on image similarity measurement, methods based on rigid and non-rigid registration, and methods for image registration using neural networks. (1) Registration methods based on feature points are divided into manually marked feature points and feature points obtained using feature extraction algorithms. Manually marked feature points are based on anatomical structures, which are highly subjective and difficult to automate. Feature points obtained using feature extraction algorithms are affected by image grayscale differences and noise, and it is difficult to find a sufficient number of corresponding feature points. (2) For optimization methods based on image similarity measurement, the registration accuracy depends on the design of similarity indicators. In multimodal images, since different imaging methods have the same anatomical structure corresponding to different imaging structures, it is a challenge to design a similarity measurement indicator with sufficient robustness. (3) Based on rigid and non-rigid registration methods, the registration accuracy depends on the registration algorithm. Improving the registration algorithm will improve the registration accuracy. (4) Methods of image registration using neural networks. Manual image registration is very cumbersome and time-consuming. The lack of sufficient labeled samples makes it difficult to meet the training requirements of the neural network. Summary of the invention

[0003] Purpose of the invention: The technical problem to be solved by the present invention is to provide a deformable registration method for three-dimensional multimodal images of rigid targets in view of the deficiencies in the prior art.

[0004] The technical problem to be solved by the present invention is the three-dimensional registration problem of multimodal imaging. Based on an improved iterative nearest neighbor registration algorithm, the normal and curvature information of the point cloud is used to construct a local feature reference coordinate system, the registration process is corrected and error points are eliminated, and a regional registration method based on structural information is used to achieve deformable rigid registration, thereby achieving more accurate registration of multimodal imaging.

[0005] In order to solve the above technical problems, the present invention discloses a deformable registration method for three-dimensional multimodal images for rigid targets, the steps of which are as follows:

[0006] Step 1, obtaining a three-dimensional multimodal image, inputting three-dimensional images of different modalities into segmentation network models adapted to the images of their respective modalities, and extracting the objects to be registered in the three-dimensional images of different modalities;

[0007] Step 2: unify the x, y, and z resolutions of the three-dimensional images of different modalities to the actual physical size of the target to be registered; construct the three-dimensional images of different modalities into a three-dimensional dense point cloud; downsample the three-dimensional dense point cloud to obtain a reference point cloud and a point cloud to be registered respectively;

[0008] Step 3, calculate the local features of the point cloud, estimate the normal vector of the reference point cloud and the point cloud to be registered using a method based on local surface fitting; calculate the surface curvature of each point of the reference point cloud and the point cloud to be registered;

[0009] Step 4, using the calculated local features of the point cloud to establish a local feature reference coordinate system, the local feature reference coordinate system is designed based on the spherical coordinate system;

[0010] Step 5: According to the structural characteristics of the point cloud to be registered, the point cloud to be registered is divided into regions in the x, y, and z directions, and the divided sub-point clouds to be registered are registered in sequence, that is, on the basis of local rigid registration, the differences between the regions are met to achieve deformable rigid registration;

[0011] Step 6, based on the improved iterative nearest neighbor algorithm of the local feature reference coordinate system established in step 4 and the deformable rigid registration method in step 5, the reference point cloud and the point cloud to be registered are registered;

[0012] Step 7: restore the registered point cloud to a three-dimensional image to complete the registration of the multimodal three-dimensional image.

[0013] In the present invention, step 1 comprises:

[0014] Step 2.1, respectively label the objects to be registered in the multi-modal images, and prepare the training sets of the objects to be registered in each modality image;

[0015] Step 2.2, using the segmentation network to train the segmentation network model under the respective modality images on different training sets of the target to be registered;

[0016] In step 2.3, the target to be registered is extracted from the multimodal image using the segmentation network model of each modality.

[0017] In the present invention, as described in step 2, unifying the x, y, and z resolutions of three-dimensional images of different modalities to the actual physical size of the target to be registered includes: adjusting the x, y resolutions of two-dimensional images of different modalities according to the actual physical size of the target to be registered on the x and y axes, wherein the adjustment is two-dimensional scaling, and the plane formed by the x and y axes becomes the scaling plane; and interpolating the z directions of images of different modalities according to the actual physical size of the target to be registered in the z direction.

[0018] In the present invention, as described in step 3, the method of estimating the normal vector based on the local surface fitting includes: for each point p in the point cloud, searching for the K nearest neighboring points, and calculating the best fitting local plane P for the K neighboring points, wherein the best fitting is to minimize the error in the least squares calculation, and the local fitting plane is expressed as:

[0019]

[0020] In the formula, is the normal vector of plane P, d is the distance from point p to the origin of the coordinate system, i=1, 2, 3, ..., K, that is, the sequence number of point p, and the normal vector of the plane fitted by the K nearest neighboring points is regarded as the normal vector of the current point.

[0021] In the present invention, in step 3, the normal vector of the plane P is obtained by principal component analysis, and the steps include: fitting the plane P through the centroid P of K neighborhood points 0 , find the center of mass P 0 The covariance matrix Q of the K neighboring points is decomposed into eigenvalues ​​to obtain the eigenvalue λ of Q 1 ,λ 2 ,λ 3 , take λ 1 ,λ 2 ,λ 3 The minimum eigenvalue λ in min The corresponding eigenvector is the normal vector of plane P. The formula for calculating the covariance matrix Q is as follows:

[0022] Where T is the matrix transpose symbol.

[0023] In the present invention, in step 3, the calculation of the surface curvature of each point includes solving the eigenvalue λ of the covariance matrix Q 1 ,λ 2 ,λ 3 , the minimum eigenvalue λ min , then the surface curvature at point p is:

[0024]

[0025] In the present invention, step 4 includes: the unit normal vector obtained for each point in step 3 is (x, y, z), the surface curvature is δ, and the normal vector and the surface curvature are converted to the established local feature reference coordinate system In the coordinate transformation method, the coordinate transformation method is:

[0026] r = δ

[0027]

[0028]

[0029] In the formula, r represents the surface curvature of point p. The smaller the value, the flatter the neighborhood of the point. θ represents the angle between the normal vector of point p and the positive direction of the z-axis. is the angle between the normal vector and the x-axis as seen from the positive z-axis.

[0030] In the present invention, step 5 comprises:

[0031] Step 5.1: Count the maximum and minimum coordinates of the point cloud M to be registered on the x, y, and z axes. The maximum value of the point cloud M to be registered in the x direction is x max , the minimum value is x min ; The maximum value in the y direction is y max , the minimum value is y min ; The maximum value in the z direction is z max , the minimum value is z min ;

[0032] Step 5.2: Observe the structural features of the point cloud to be registered and determine the number of divisions of the point cloud in each direction, such as n in the x, y, and z directions respectively. x , n y , n z sub-point clouds, then the point cloud M to be registered is divided into N sub-point clouds, N = n x *n y *n z , N sub-point clouds are represented as M 1 、M 2 、M 3 ......M N ,have

[0033] In the present invention, step 6 comprises:

[0034] Step 6.1: In a three-dimensional rectangular coordinate system, use the Euclidean distance to calculate each sub-point cloud M of the point cloud M to be registered i (i=1,2,3……N) The nearest neighbor point set C in the reference point cloud F, the sub-point cloud M to be registered iThere is a one-to-one correspondence between the points in the nearest neighbor point set C; the specific implementation is to traverse M i , calculate the Euclidean distance between the point and all the points in the reference point cloud F, select the point with the smallest distance as the corresponding nearest neighbor point of the point, and the obtained series of points constitute the nearest neighbor point set C;

[0035] Step 6.2: Modify the nearest neighbor point set C in step 6.1 using the local feature reference coordinate system constructed in step 4, and define the similarity S between two points in the local feature reference coordinate system. sph ,

[0036]

[0037] Where point p in the local feature reference coordinate system 1 、p 2 The descriptions are For each point p in the nearest neighbor point set C obtained in step 6.1 c Find its neighborhood point set K in the reference point cloud F p0 , and point p c The corresponding sub-point cloud M i The point in Ask separately With the neighborhood point set K p0 The similarity S of all points in sph , select similarity S sph The smallest point Replace point p c , perform the above correction on each point in the nearest neighbor point set C to obtain the corrected nearest neighbor point set C′;

[0038] Step 6.3: Use the local feature reference coordinate system constructed in step 4 to remove error points from the corrected nearest neighbor point set C′; i And the corrected nearest neighbor point set C′ is counted in local coordinates r,θ, value, and find both r, θ, The difference in values ​​Δr, Δθ, Setting the threshold:

[0039]

[0040]

[0041]

[0042] Among them, r threshold Represents the sub-point cloud M to be registered iThe average value of the maximum and minimum values ​​of the difference between the corrected nearest neighbor point set C′ in the local coordinates, θ threshold Represents the sub-point cloud M to be registered i The average of the maximum and minimum values ​​of the difference between the values ​​of θ and the corrected nearest neighbor point set C′ in the local coordinates, Represents the sub-point cloud M to be registered i And the corrected nearest neighbor point set C′ in local coordinates The average of the maximum and minimum values ​​of the difference between the values ​​retains Δr, Δθ, are all less than the corresponding threshold r threshold ,θ threshold , The sub-point cloud M to be registered i The remaining points are removed as error points, and the corrected nearest neighbor point set C′ is obtained. 0 ;

[0043] Step 6.4: According to the nearest neighbor point set C′ obtained in step 6.3 after removing the error and correcting it 0 and the sub-point cloud M to be registered i , the singular value decomposition method is used to calculate the rotation matrix and translation matrix of the point cloud, and the improved iterative closest point algorithm is used to iteratively optimize the rotation matrix and translation matrix. The loss function of the iterative process is defined as,

[0044]

[0045] in, Represents the sub-point cloud M to be registered i The total number of midpoints, d k Represents the point set C′ 0 and point set M i The Euclidean distance between corresponding points in k , Δθ k , Respectively represent the point set C′ in the local feature reference coordinate system 0 and point set M i The corresponding point r, θ, When the Loss value is less than the set threshold, the iteration is stopped and the point set M is obtained. i and point set C′ 0 The final rotation matrix R of the registration i and the translation matrix T i ;

[0046] Step 6.5, for each sub-point cloud M described in step 5.2 1 、M 2 、M 3 ......M N, respectively perform the operations of step 6.1 to step 6.4 to obtain the rotation matrix R and translation matrix T of the point cloud M to be registered, and then obtain the result M after the registration of the point cloud M to be registered F :

[0047]

[0048]

[0049] In the present invention, step 7 includes,

[0050] Step 7.1: Compare the reference point cloud F and the registered point cloud M F , according to the number of interpolation in the z direction in step 2, the reference point cloud F and the registered point cloud M in the z direction F To stratify;

[0051] Step 7.2, treating the points in each layer of the point cloud in step 7.1 as a sphere with a certain radius, projecting each layer of the point cloud onto the xy plane, and restoring it to a two-dimensional image sequence, that is, restoring it to a registered three-dimensional multimodal image;

[0052] In step 7.3, the similarity measure between sets is used to measure the registration effect of 3D images of different modalities.

[0053] Beneficial effects: The present invention provides a deformable registration method for three-dimensional multimodal images for rigid targets. Point clouds are constructed through three-dimensional images, and local feature reference coordinate systems are constructed using the normal and curvature information of the point clouds. The registration process is corrected and error points are eliminated, and the multiple features of the point clouds are fully utilized to improve the robustness of the registration. The regional registration method based on structural information is used to meet the differences between regions on the basis of local rigid registration, that is, non-rigid changes are added to the rigid registration, deformable rigid registration is achieved, and the registration accuracy is improved. The existing iterative nearest neighbor registration method lacks the elimination of error points, resulting in errors in the rotation matrix and translation matrix obtained by iteration. Therefore, the present invention makes full use of the multiple features of the point cloud to construct a local feature reference coordinate system to improve the registration process, and combines the point cloud structure information to partition the point cloud to achieve a better registration effect between the registration point cloud and the reference point cloud. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] The present invention will be further described in detail below in conjunction with the accompanying drawings and specific embodiments, and the above and / or other advantages of the present invention will become more clear.

[0055] Figure 1 It is a schematic diagram of the process of the present invention.

[0056] Figure 2It is a schematic diagram of an improved method of the iterative closest point algorithm (ICP) in the present invention. DETAILED DESCRIPTION

[0057] Example

[0058] like Figure 1 As shown, the present invention discloses a deformable registration method for three-dimensional multimodal images for rigid targets, comprising the following steps:

[0059] Step 1: In this embodiment, pelvic region images under different modalities are registered, namely, computed tomography (CT) images and magnetic resonance imaging (MRI) images, and the registration target is the pelvic region in the image. First, the acquired CT and MRI images are respectively input into two pre-trained pelvic segmentation network models to extract the pelvic region in the CT and MRI images respectively;

[0060] Step 2: Unify the resolution in each direction and construct a point cloud. The two directions that constitute the cross section of the CT and MRI images are called the x direction and the y direction, and the direction perpendicular to the cross section is called the z direction. The resolution in the three directions is unified to the actual physical size of the human pelvis; the three-dimensional CT and MRI images are respectively constructed as three-dimensional dense point clouds; the three-dimensional dense point clouds are downsampled to obtain a reference point cloud and a point cloud to be registered, respectively. The reference point cloud refers to the pelvic CT point cloud, and the point cloud to be registered refers to the pelvic MRI point cloud;

[0061] Step 3, calculate the local features of the point cloud, estimate the normal vector of the pelvic CT point cloud and the pelvic MRI point cloud using a method based on local surface fitting; and calculate the surface curvature of each point of the two point clouds;

[0062] Step 4, using the calculated local features of the point cloud to establish a local feature reference coordinate system, the local feature reference coordinate system is designed based on the spherical coordinate system;

[0063] Step 5: According to the structural characteristics of the pelvic MRI point cloud, including roughly left-right symmetry and rigidity in the z direction, the pelvic MRI point cloud is divided into regions along the x and y directions, and the divided pelvic MRI sub-point clouds are registered in turn, that is, the differences between the regions are satisfied on the basis of local rigid registration, so as to achieve deformable rigid registration;

[0064] Step 6: Based on the improved iterative nearest neighbor algorithm of the local feature reference coordinate system established in step 4 (such as Figure 2 As shown) and the deformable rigid registration method in step 5, the pelvic CT point cloud and the pelvic MRI point cloud are registered;

[0065] Step 7, restore the registered pelvic MRI point cloud to a three-dimensional image, and realize the registration of the pelvis in the CT and MRI three-dimensional images.

[0066] In this embodiment, in step 1, the LabelMe annotation tool is first used to outline the pelvic area in the CT and MRI images to produce pelvic segmentation data sets of CT and MRI, and two segmentation network models are used for training respectively. The segmentation network model is a U-Net segmentation network. Refer to (Ronneberger O, Fischer P, Brox TU-Net: Convolutional networks for biomedical image segmentation.arXiv[J].Lecture Notes in Computer Science, 2015, 2015.). During the U-Net network training process, the input image is first encoded and then decoded to output the segmentation result. According to the difference between this result and the actual segmentation, the segmentation network is trained by back propagation.

[0067] In this embodiment, in step 2, the resolution of the CT and MRI images in the x, y, and z directions is unified to the actual physical size of the pelvis, including scaling and interpolation. In this embodiment, the CT image obtained is 512*512 pixels, and there are about 100 CT images of the pelvic area of ​​a case; the MRI image obtained is 784*784 pixels, and there are about 40 MRI images of the pelvic area of ​​a case. There is a difference between the size of the pelvic area extracted from the CT and MRI images in step 1 and the actual pelvic size. The scaling refers to adjusting the x, y resolution of the CT and MRI images according to the actual physical size of the pelvis in the x, y directions, and the adjustment is two-dimensional scaling; the interpolation is to interpolate the z direction of the CT and MRI images according to the actual physical size of the pelvis in the z direction, and the interpolation method is bicubic interpolation, and the CT and MRI images are interpolated by 2 times and 5 times, respectively.

[0068] In this embodiment, in step 3, the method based on local surface fitting for normal vector estimation and surface curvature calculation includes:

[0069] Step 3.1, for each point p in the pelvic CT point cloud and MRI point cloud, search for its nearest K neighboring points, and calculate the best fitting local plane P for the K neighboring points, where the best fitting is to minimize the error in the least squares calculation. This local fitting plane can be expressed as:

[0070]

[0071] In the formula, is the normal vector of plane P, d is the distance from point p to the origin of the coordinate system, i=1,2,3,……,K, which is the serial number of point p. The normal vector of the plane fitted by the K nearest neighboring points can be regarded as the normal vector of the current point.

[0072] The normal vector of the plane P can be obtained by principal component analysis, and the steps include: fitting the plane P through the centroid P of K neighborhood points 0 , find the center of mass P 0 The covariance matrix Q of the K neighboring points is decomposed into eigenvalues ​​to obtain the eigenvalue λ of Q 1 ,λ 2 ,λ 3 . Take λ 1 ,λ 2 ,λ 3 The minimum eigenvalue λ in min The corresponding eigenvector is the normal vector of plane P. The formula for calculating the covariance matrix Q is as follows:

[0073]

[0074] Step 3.2, the calculation of the surface curvature of each point includes solving the eigenvalue λ of the covariance matrix Q 1 ,λ 2 ,λ 3 , the minimum eigenvalue λ min , then the surface curvature at point p is:

[0075]

[0076] In this embodiment, in step 4, the establishment of a local feature reference coordinate system using the calculated point cloud local features includes converting the unit normal vector (x, y, z) obtained in step 3.1 and the surface curvature δ obtained in step 3.2 to the local feature reference coordinate system In the coordinate transformation method, the coordinate transformation method is:

[0077] r = δ

[0078]

[0079]

[0080] In the formula, r represents the surface curvature of point p. The smaller the value, the flatter the neighborhood of the point. θ represents the angle between the normal vector of point p and the positive direction of the z-axis. is the angle between the normal vector and the x-axis as seen from the positive z-axis.

[0081] In this embodiment, in step 5, the step of dividing the pelvic MRI point cloud into regions in each axis direction according to the structural characteristics of the pelvic MRI point cloud includes:

[0082] Step 5.1: Count the maximum and minimum coordinates of the pelvic MRI point cloud M on the x, y, and z axes, and record the maximum value in the x direction as x max , the minimum value is x min ; The maximum value in the y direction is y max , the minimum value is y min ; The maximum value in the z direction is z max , the minimum value is z min ;

[0083] Step 5.2, observe the structural features of the point cloud to be registered, including the roughly symmetric left and right and the rigid features in the z direction, so as to determine the area division of the pelvic MRI point cloud along the x and y directions and the number of divisions in each direction, that is, divide the point cloud into n areas in the x, y, and z directions respectively. x =2,n y =2,n z = 1 sub-point cloud, then the pelvic MRI point cloud M is divided into N sub-point clouds, N = n x *n y *n z =4. The 4 sub-point clouds are represented as M 1 、M 2 、M 3 、M 4 ,have

[0084] In this embodiment, in step 6, an improved iterative nearest neighbor algorithm (such as Figure 2 As shown) and the deformable rigid registration method in step 5 include:

[0085] Step 6.1, in a three-dimensional rectangular coordinate system, use the Euclidean distance to calculate each sub-point cloud M of the pelvic MRI point cloud M i (i=1, 2, 3, 4) The nearest neighbor point set C in the pelvic CT point cloud F, each sub-point cloud M i There is a one-to-one correspondence between the points in the nearest neighbor point set C, which is implemented by traversing M i Calculate the Euclidean distance between the point and all points in the pelvic CT point cloud F, and select the point with the smallest distance as the corresponding nearest neighbor point of the point. Such a series of points constitute the nearest neighbor point set C;

[0086] Step 6.2: Modify the nearest neighbor point set C in step 6.1 using the local feature reference coordinate system constructed in step 4, and define the similarity S between two points in the local feature reference coordinate system. sph ,

[0087]

[0088] Where point p in the local feature reference coordinate system 1 、p 2 The descriptions are For each point p in the nearest neighbor point set C obtained in step 6.1 c Find its neighborhood point set K in the reference point cloud F p0 , and point p c The corresponding sub-point cloud M i The point in Ask separately With the neighborhood point set K p0 The similarity S of all points in sph , select similarity S sph The smallest point Replace point p c , perform the above correction on each point in the nearest neighbor point set C to obtain the corrected nearest neighbor point set C′;

[0089] Step 6.3, using the local feature reference coordinate system constructed in step 4, remove the error points from the corrected nearest neighbor point set C′; for the pelvic MRI sub-point cloud M to be registered in step 6.1 i And the corrected nearest neighbor point set C′ is counted in local coordinates r,θ, value, and find both r, θ, The difference in values ​​Δr, Δθ, Setting the threshold:

[0090]

[0091]

[0092]

[0093] Among them, r threshold Represents the sub-point cloud M to be registered i The average value of the maximum and minimum values ​​of the difference between the corrected nearest neighbor point set C′ in the local coordinates, θ threshold Represents the sub-point cloud M to be registered i The average of the maximum and minimum values ​​of the difference between the values ​​of θ and the corrected nearest neighbor point set C′ in the local coordinates, Represents the sub-point cloud M to be registered i And the corrected nearest neighbor point set C′ in local coordinates The average of the maximum and minimum values ​​of the difference between the values ​​retains Δr, Δθ, are all less than the corresponding threshold r threshold ,θ threshold , The sub-point cloud M to be registered iThe remaining points are removed as error points, and the corrected nearest neighbor point set C′ is obtained. 0 ;

[0094] Step 6.4: According to the nearest neighbor point set C′ obtained in step 6.3 after removing the error and correcting it 0 and the MRI sub-point cloud M to be registered i , the singular value decomposition method is used to calculate the rotation matrix and translation matrix of the point cloud, and the improved iterative closest point algorithm is used to iteratively optimize the rotation matrix and translation matrix. The loss function of the iterative process is defined as,

[0095]

[0096] in, Represents the sub-point cloud M to be registered i The total number of midpoints, d k Represents the point set C′ 0 and point set M i The Euclidean distance between corresponding points in k , Δθ k , Respectively represent the point set C′ in the local feature reference coordinate system 0 and point set M i The corresponding point r, θ, When the Loss value is less than the set threshold, the iteration is stopped and the point set M is obtained. i and point set C′ 0 The final rotation matrix R of the registration i and the translation matrix T i ;

[0097] Step 6.5, for each sub-point cloud M described in step 5.2 1 、M 2 、M 3 、M 4 Perform the operations from step 6.1 to step 6.4 respectively to obtain the rotation matrix R and translation matrix T of the pelvic MRI point cloud M, and then obtain the result M after the registration of the pelvic MRI point cloud M. F .

[0098]

[0099]

[0100]

[0101] In this embodiment, step 7 includes:

[0102] Step 7.1: Compare the pelvic CT point cloud F and the registered pelvic MRI point cloud MF According to the number of layers before interpolation in the z direction in step 2, the pelvic CT point cloud F and the registered pelvic MRI point cloud M in the z direction are interpolated. F To stratify;

[0103] Step 7.2, treating the points in each layer of the point cloud in step 7.1 as a sphere with a certain radius, projecting each layer of the point cloud onto the xy plane, and restoring it to a three-dimensional MRI image after registration;

[0104] Step 7.3, using the similarity measure between sets to measure the registration effect of CT image and MRI image. The above steps realize the registration of three-dimensional multimodal imaging and achieve a relatively ideal registration effect.

[0105] The present invention provides a method and idea for a deformable registration method of three-dimensional multimodal images for rigid targets. There are many methods and ways to implement the technical solution. The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention. All components not specified in this embodiment can be implemented by existing technologies.

Claims

1. A deformable registration method for three-dimensional multimodal images for rigid targets, characterized in that: The steps include: Step 1, obtaining a three-dimensional multimodal image, inputting three-dimensional images of different modalities into segmentation network models adapted to the images of their respective modalities, and extracting the objects to be registered in the three-dimensional images of different modalities; Step 2: unify the x, y, and z resolutions of the three-dimensional images of different modalities to the actual physical size of the target to be registered; construct the three-dimensional images of different modalities into a three-dimensional dense point cloud; Downsampling the three-dimensional dense point cloud to obtain a reference point cloud and a point cloud to be registered; Step 3, calculate the local features of the point cloud, estimate the normal vector of the reference point cloud and the point cloud to be registered using a method based on local surface fitting; calculate the surface curvature of each point of the reference point cloud and the point cloud to be registered; Step 4, using the calculated local features of the point cloud to establish a local feature reference coordinate system, the local feature reference coordinate system is designed based on the spherical coordinate system; Step 5: According to the structural characteristics of the point cloud to be registered, the point cloud to be registered is divided into regions in the x, y, and z directions, and the divided sub-point clouds to be registered are registered in sequence, that is, on the basis of local rigid registration, the differences between the regions are satisfied to achieve deformable rigid registration; Step 6, based on the improved iterative nearest neighbor algorithm of the local feature reference coordinate system established in step 4 and the deformable rigid registration method in step 5, the reference point cloud and the point cloud to be registered are registered; Step 7, restoring the registered point cloud to a three-dimensional image, completing the registration of the multi-modal three-dimensional image; Wherein, step 6 comprises: Step 6.1: In a three-dimensional rectangular coordinate system, use the Euclidean distance to calculate each sub-point cloud M of the point cloud M to be registered i (i=1,2,3……N) The nearest neighbor point set C in the reference point cloud F, the sub-point cloud M to be registered i There is a one-to-one correspondence between the points in the nearest neighbor point set C; the specific implementation is to traverse M i , calculate the Euclidean distance between the point and all the points in the reference point cloud F, select the point with the smallest distance as the corresponding nearest neighbor point of the point, and the obtained series of points constitute the nearest neighbor point set C; Step 6.2: Modify the nearest neighbor point set C in step 6.1 using the local feature reference coordinate system constructed in step 4, and define the similarity S between two points in the local feature reference coordinate system. sph , The descriptions of points p1 and p2 in the local feature reference coordinate system are For each point p in the nearest neighbor point set C obtained in step 6.1 c Find its neighborhood point set K in the reference point cloud F p0 , and point p c The corresponding sub-point cloud M i The point in Ask separately With the neighborhood point set K p0 The similarity S of all points in sph , select similarity S sph The smallest point Replace point p c , perform the above correction on each point in the nearest neighbor point set C to obtain the corrected nearest neighbor point set C′; Step 6.3: Use the local feature reference coordinate system constructed in step 4 to remove error points from the corrected nearest neighbor point set C′; i And the corrected nearest neighbor point set C′ is counted in the local coordinates r,θ, value, and find both r, θ, The difference in values ​​Δr, Δθ, Setting the threshold: Among them, r threshold Represents the sub-point cloud M to be registered i The average value of the maximum and minimum values ​​of the difference between the corrected nearest neighbor point set C′ in the local coordinates, θ threshold Represents the sub-point cloud M to be registered i The average of the maximum and minimum values ​​of the difference between the values ​​of θ and the corrected nearest neighbor point set C′ in the local coordinates, Represents the sub-point cloud M to be registered i And the corrected nearest neighbor point set C′ in local coordinates The average of the maximum and minimum values ​​of the difference between the values ​​retains Δr, Δθ, are all less than the corresponding threshold r threshold ,θ threshold , The sub-point cloud M to be registered i The points in the corrected nearest neighbor point set C′ are taken as the error points and the remaining points are removed to obtain the corrected nearest neighbor point set C′0 after the error is removed; Step 6.4: Based on the nearest neighbor point set C′0 obtained in step 6.3 after removing the error and correcting it, and the sub-point cloud M to be registered i , the singular value decomposition method is used to calculate the rotation matrix and translation matrix of the point cloud, and the improved iterative closest point algorithm is used to iteratively optimize the rotation matrix and translation matrix. The loss function of the iterative process is defined as, in, Represents the sub-point cloud M to be registered i The total number of midpoints, d k Denote the point set C′0 and the point set M i The Euclidean distance between corresponding points in k , Δθ k , They represent the point set C′0 and the point set M in the local feature reference coordinate system respectively. i The corresponding point r, θ, When the Loss value is less than the set threshold, the iteration is stopped and the point set M is obtained. i The final rotation matrix R registered with the point set C′0 i and the translation matrix T i ; Step 6.5, for each sub-point cloud M1, M2, M3, ..., M in step 5 N , respectively perform the operations of step 6.1 to step 6.4 to obtain the rotation matrix R and translation matrix T of the point cloud M to be registered, and then obtain the result M after the registration of the point cloud M to be registered F :

2. The deformable registration method for three-dimensional multimodal images for rigid targets according to claim 1, characterized in that: Step 1 includes: Step 2.1, respectively label the objects to be registered in the multi-modal images, and prepare the training sets of the objects to be registered in each modality image; Step 2.2, using the segmentation network to train the segmentation network model under the respective modality images on different training sets of the target to be registered; In step 2.3, the target to be registered is extracted from the multimodal image using the segmentation network model of each modality.

3. The deformable registration method for three-dimensional multimodal images for rigid targets according to claim 2, characterized in that: As described in step 2, unifying the x, y, and z resolutions of three-dimensional images of different modalities to the actual physical size of the target to be registered includes: adjusting the x, y resolutions of two-dimensional images of different modalities according to the actual physical size of the target to be registered on the x and y axes, wherein the adjustment is two-dimensional scaling, and the plane formed by the x and y axes becomes the scaling plane; and interpolating the z directions of images of different modalities according to the actual physical size of the target to be registered in the z direction.

4. The deformable registration method for three-dimensional multimodal images for rigid targets according to claim 3, characterized in that: As described in step 3, the normal vector estimation based on the local surface fitting method includes: for each point p in the point cloud, searching for its nearest K neighboring points, calculating the best fitting local plane P for the K neighboring points, the best fitting is to minimize the error in the least squares calculation, and this local fitting plane is expressed as: In the formula, is the normal vector of plane P, d is the distance from point p to the origin of the coordinate system, i=1, 2, 3, ..., K, that is, the sequence number of point p, and the normal vector of the plane fitted by the K nearest neighboring points is regarded as the normal vector of the current point.

5. The deformable registration method for three-dimensional multimodal images for rigid targets according to claim 4, characterized in that: In step 3, the normal vector of plane P is obtained by principal component analysis, and the steps include: fitting plane P through the centroid P0 of K neighboring points, obtaining the covariance matrix Q of the centroid P0 and the K neighboring points, performing eigenvalue decomposition on Q, obtaining the eigenvalues ​​λ1, λ2, λ3 of Q, and taking the minimum eigenvalue λ among λ1, λ2, λ3 min The corresponding eigenvector is the normal vector of plane P. The formula for calculating the covariance matrix Q is as follows: Where T is the matrix transpose symbol.

6. The deformable registration method for three-dimensional multimodal images for rigid targets according to claim 5, characterized in that: In step 3, the calculation of the surface curvature of each point includes solving the eigenvalues ​​λ1, λ2, λ3 of the covariance matrix Q, the minimum eigenvalue λ min , then the surface curvature at point p is:

7. The deformable registration method for three-dimensional multimodal images for rigid targets according to claim 6, characterized in that: In step 4 Including, the unit normal vector (x, y, z) for each point obtained in step 3, the surface curvature is δ, and the normal vector and surface curvature are converted to the established local feature reference coordinate system In the coordinate transformation method, the coordinate transformation method is: r = δ In the formula, r represents the surface curvature of point p. The smaller the value, the flatter the neighborhood of the point. θ represents the angle between the normal vector of point p and the positive direction of the z-axis. is the angle between the normal vector and the x-axis as seen from the positive z-axis.

8. The deformable registration method for three-dimensional multimodal images for rigid targets according to claim 7, characterized in that: Step 5 includes: Step 5.1: Count the maximum and minimum coordinates of the point cloud M to be registered on the x, y, and z axes. The maximum value of the point cloud M to be registered in the x direction is x max , the minimum value is x min ; The maximum value in the y direction is y max , the minimum value is y min ; The maximum value in the z direction is z max , the minimum value is z min ; Step 5.2: Observe the structural features of the point cloud to be registered and determine the number of divisions of the point cloud in each direction, such as n in the x, y, and z directions respectively. x , n y , n z sub-point clouds, then the point cloud M to be registered is divided into N sub-point clouds, N = n x *n y *n z , the N sub-point clouds are represented as M1, M2, M3...M N ,have 9. The deformable registration method for three-dimensional multimodal images for rigid targets according to claim 8, characterized in that: Step 7 includes, Step 7.1: Compare the reference point cloud F and the registered point cloud M F , according to the number of interpolation in the z direction in step 2, the reference point cloud F and the registered point cloud M in the z direction F To stratify; Step 7.2, treating the points in each layer of the point cloud in step 7.1 as a sphere with a certain radius, projecting each layer of the point cloud onto the xy plane, and restoring it to a two-dimensional image sequence, that is, restoring it to a registered three-dimensional multimodal image; In step 7.3, the similarity measure between sets is used to measure the registration effect of 3D images of different modalities.

Images