A method of non-rigid registration of medical images

Abstract

The application provides a medical image non-rigid registration method, and the specific steps include: (1) input floating and fixed images, and send the images into a pyramid feature encoder to generate feature maps with different sampling rates; (2) starting from the highest sampling rate feature, generate an adaptive displacement space, a displacement field and estimate an uncertainty in stages; (3) according to the uncertainty of the last stage and the floating feature of the current stage, form a displacement base and a displacement fine adjustment value respectively, and then add them to obtain an adaptive displacement space; (4) according to the deformation feature, the target feature and the displacement space of the current stage, obtain a displacement field through a multi-head correlation, cost aggregation, normalization and displacement fusion step, and estimate the uncertainty of the field at the same time; (5) repeat steps (2) to (4) until a displacement field with the same size as the original image is obtained, and then obtain the final registration result after deforming the image.

Description

Technical Field

[0001] This invention belongs to the field of medical image processing technology, and relates to a medical image registration method, particularly a non-rigid, unsupervised medical image registration method based on deep learning. Background Technology

[0002] Medical image registration involves finding a spatial transformation to align or match all anatomical points or medically diagnostic reference points in two images. Specifically, registration methods can be broadly categorized into rigid transformations and flexible transformations based on the transformation relationship. Rigid transformations are essentially global linear transformations, with typical methods including rotation, scaling, and affine transformations. However, their drawbacks include the inability to simulate local geometric differences, and the lack of robustness to deformation due to the simplicity of the model parameters. Non-rigid transformations overcome these shortcomings by using a dense displacement field (also known as a deformation field or transformation field) as the model parameters. In this case, each pixel or voxel has its own displacement vector, thus exhibiting local deformation capabilities and being widely used in the registration of complex and variable medical images.

[0003] For medical image registration, researchers have designed many traditional registration algorithms to automate this process using computers. However, these traditional methods still have limitations, such as slow computation speed, complex registration procedures, and inability to utilize massive amounts of data. Compared to traditional registration methods, deep learning-based methods can shift the long iteration time to model training time, resulting in significantly shorter actual inference time compared to traditional methods. Furthermore, end-to-end deep regression network design replaces the complex procedures of traditional registration methods, greatly reducing the difficulty of the registration process. Additionally, deep learning networks can obtain more generalized and robust registration results by fitting massive amounts of data.

[0004] Due to the scarcity of gold standard annotations and the high cost of manual annotation, the development of supervised medical image registration is constrained. Unsupervised medical image registration, however, has lower or no annotation requirements and holds greater promise. Unsupervised deep learning image registration originated from Jaderberg's application of Spatial Transformer (STN) to registration networks. Based on the deformation field predicted by the network, the source image is interpolated and sampled to obtain the deformed image. This allows for the calculation of image similarity loss during training, forming a loss function that drives the backpropagation of the deep network, achieving the training goal. Balakrishnan's VoxelMorph registration framework utilizes a 3D U-Net to output a global 3D deformation field without generating intermediate displacement parameters. This method's excellent performance has garnered significant attention for unsupervised deep learning registration methods. However, this framework does not perform well with images exhibiting large deformations because the network's ability to capture long-distance deformations is insufficient. Therefore, many methods based on a "coarse-to-fine" approach have been proposed. Zhao employs a cascaded approach, performing progressive registration on the images one after another. Each stage of the cascade is a VoxelMorph. While this method is effective, it has a large number of model parameters, which limits the registration performance and speed. Hu proposed a two-stream pyramid registration network called Dual-PRNet, which utilizes image pyramids to directly generate multiple registration fields of different scales from the feature maps generated by the decoder. Summary of the Invention

[0005] The technical problem this invention aims to solve is that current methods are ineffective at registering images with significant deformation, and the commonly used "coarse-to-fine" approach suffers from error accumulation. Therefore, a non-rigid registration method for medical images based on adaptive correlation and uncertainty estimation is proposed.

[0006] The technical solution of the present invention is as follows:

[0007] A non-rigid registration method for medical images, comprising the following steps:

[0008] Step (1) Pyramid Feature Encoding

[0009] The input to this step is a floating image and a fixed image, which are then fed into a pyramid feature encoding network with shared parameters to obtain floating or fixed feature maps at different sampling ratios. The pyramid feature encoding network generates feature maps at different sampling ratios, i.e., a feature pyramid. The main component of the feature encoding network is a 3D convolutional layer, stacked using the residual block structure of ResNet as the basic unit, and different sampling ratios are obtained by setting the stride of the convolutional layers. This encoding network can effectively autoencode 3D medical images and extract features suitable for registration. The following steps take each feature map at a different ratio as a stage, starting with the highest sampling ratio, i.e., the smallest size feature map, and generate larger and more precise fields stage by stage.

[0010] Step (2) Preprocessing

[0011] The period from the start of step (2) to the end of step (4) when the displacement field is generated is defined as one stage.

[0012] This step first upsamples the displacement field φ and uncertainty σ from the previous stage to obtain the same size as the feature map; if it is the first stage, the displacement field and uncertainty are initialized to zero vectors. Next, the floating feature of the current stage obtained in step (1) is deformed, which can be done using only the upsampled displacement field φ, resulting in the deformed floating feature of the current stage, referred to simply as the deformation feature:

[0013] Feats warp (c,x,y,z)=φ(x,y,z)⊙Feats mov (c,x,y,z)

[0014] Where x, y, and z represent the voxel positions, corresponding to the three coordinate axes of the 3D image; Feats mov and Feats warp These are the floating feature and the deformation feature, respectively. The displacement field φ contains three 3D matrices, each with the same size as the floating feature, representing the relationship matrix for transforming from the floating feature to the fixed feature. ⊙ represents the deformation operation, and the deformation operation performed by different channels is the same.

[0015] Step (3) Generate displacement space

[0016] The input for this step consists of two parts: one is the floating feature map of the current stage, which is fed into a multi-layer convolutional layer to generate shift fine-tuning values, specifically in the form of a four-dimensional tensor.

[0017] Disp shift (t,x,y,z)=Convs[Feats mov (c,x,y,z)]

[0018] Among them, Disp shift represents the displacement fine-tuning value, and t represents the number of points in the displacement space; c represents the number of channels in the feature map; Convs represents a series of self-network structures, consisting of multiple layers of 3D convolution, instance normalization, and non-linear activation functions.

[0019] Another part of the input for this step is the uncertainty from the previous stage, which is used to generate the displacement basis value through a linear transformation.

[0020] Disp space (t,x,y,z)=[α*σ(x,y,z)+β]*space(t)

[0021] Among them, Disp space σ represents the displacement basis value; α and β represent the parameters of the linear transformation; space represents a cubic sampling space in three-dimensional space with t points; before performing the second multiplication, the tensor needs to be extended to the same dimension.

[0022] Summing the displacement base value and the displacement fine-tuning value yields a coordinate vector, which is the final discrete displacement space:

[0023] Disp(t, x, y, z) = Disp space (t, x, y, z) + Disp shift (t, x, y, z)

[0024] Where Disp represents the final displacement space, specifically a vector of coordinate values.

[0025] The displacement space is dynamically transformed based on the prior structural information of the floating image, thus enabling registration to bypass noise and avoid discontinuities in the displacement field; at the same time, it can accurately expand the sampling range, effectively cope with image registration with large deformations, and thus improve the accuracy of the initial matching.

[0026] Step (4) Estimation of multiple positions

[0027] The input for this step consists of two parts. The first part is the deformation features of the current stage obtained through step (2) and the fixed feature map obtained through step (1). The other part is the adaptive coordinates (displacement space Disp) generated in step (3). After multi-head correlation, cost aggregation, normalization and displacement fusion, the displacement field is obtained, and the uncertainty of the field is estimated at the same time. Specifically, the deformation features Feats are first... warp Resampling is performed based on adaptive coordinates, note that this is done in t rounds, therefore the generated features will have an additional t dimension:

[0028] Featsdisp (c,t,x,y,z)=Disp(t,x,y,z)⊙Feats warp (c, x, y, z)

[0029] Feats disp This represents the multi-displacement features obtained from sampling.

[0030] Feats with multi-displacement and fixed features fix Dividing the channel into k headers is essentially grouping the channel, i.e., the K_divide operation:

[0031]

[0032]

[0033] Then Perform relevant operations along the dimension to obtain the k-head similarity cost body Sim:

[0034]

[0035] in This indicates the inner product operation.

[0036] Then, cost aggregation is performed using a k-head convolutional network to finally obtain the cost volume (Cost).

[0037] Cost(k,t,x,y,z)=K_Convs[Sim(k,t,x,y,z)]

[0038] Where K_Convs represents a series of self-network structures, consisting of multiple layers of 3D convolutions, instance normalization, and non-linear activation functions, and the number of groups in the convolutional layers is k.

[0039] Then, normalize the displacements generated by all different heads to obtain the probability vector Pro corresponding to the displacement space:

[0040] Pro(k,t,x,y,z)=Softmax[Cost(k,t,x,y,z)]

[0041] Here, Softmax represents the normalization function, which is performed on both k and t dimensions.

[0042] Finally, the displacement space Disp is multiplied by the probability vector Pro to obtain the fused displacement field φ:

[0043]

[0044] in This represents the inner product operation, performed in both k and t dimensions, with the tensor expanded beforehand. Since the displacement field φ is essentially a first-order moment result, the square root of the second-order central moment yields the uncertainty σ. The specific calculation process is as follows:

[0045]

[0046] Once a displacement field and uncertainty of the same scale as the floating feature map are generated, the stage is considered complete and the next stage will begin. Jump to step (2); if all stages are completed, proceed to step (5).

[0047] Step (5) Construct the unsupervised loss function

[0048] After obtaining the displacement field φ at the same scale as the original image, the floating image I is then processed. Mov The deformation is performed, and the resulting deformed image is the final registration result I. φ(Mov) :

[0049] I φ(Mov) =φ(x,t,z)⊙I Mov

[0050] Since training a neural network model requires constructing a loss function, this invention uses an unsupervised loss function, which consists of two parts: one is a similarity loss function. Second, the smoothness loss function of the displacement field. The sum of the two yields the final loss function, where λ is the balancing parameter, also known as the regularization coefficient.

[0051]

[0052] The beneficial effects of this invention are as follows: This method uses adaptive correlation, i.e., generating an adaptive displacement space to accurately control the sampling range of features, which can effectively improve the registration effect of images with large deformations. It can also cope with noise interference, weak texture and displacement field discontinuity, and improve the similarity after registration. In addition, uncertainty estimation is used to alleviate the problem of error accumulation, and a multi-head mechanism is used to deepen the registration. Attached Figure Description

[0053] Figure 1 This is a basic schematic diagram of the method of the present invention.

[0054] Figure 2 The overall flowchart provided for the method of the present invention.

[0055] Figure 3 A detailed schematic diagram of step (3) of the method of the present invention.

[0056] Figure 4 A detailed schematic diagram of step (4) of the method of the present invention.

[0057] Figure 5 An example diagram of the displacement field generated by the registration method of the present invention.

[0058] Figures 6(a) and 6(b) are examples of image comparison of the registration results provided by the method of the present invention, wherein Figure 6(a) is the fixed image (target) and Figure 6(b) is the deformed image (result). Detailed Implementation

[0059] The specific embodiments of the present invention will be further described below with reference to the accompanying drawings and technical methods.

[0060] The basic process of a non-rigid registration method for medical images according to the present invention is as follows: Figure 1 and Figure 2 As shown, the basic steps include: (1) inputting floating and fixed images and sending them to the pyramid feature encoder to generate feature maps with different sampling ratios; (2) starting from the feature with the highest sampling ratio (i.e. the smallest size), generating adaptive displacement space and displacement field in stages, and estimating the uncertainty; (3) forming displacement base and displacement fine-tuning value respectively based on the uncertainty of the previous stage and the floating feature of the current stage, and then summing them to obtain the adaptive displacement space; (4) obtaining the displacement field after multi-head correlation, cost aggregation, normalization and displacement fusion steps based on the deformation feature, target feature and displacement space of the current stage, and estimating the uncertainty of the field at the same time; (5) repeating steps (2) to (4) until obtaining a displacement field with the same size as the original image, and obtaining the final registration result after deformation.

[0061] This embodiment uses respiratory biphasic registration of lung CT images. The specific registration process is as follows:

[0062] Step (1) Pyramid Feature Encoding

[0063] This step inputs the inspiratory and expiratory phases of a lung CT scan, as floating and fixed images, respectively (e.g., ...). Figure 2 (As shown on the left) The feature encoder that sends the shared parameters (such as...) Figure 2 (As shown by the arrows at the top and bottom of the image), four feature maps of 1 / 2, 1 / 4, 1 / 8, and 1 / 16 of the input image size are obtained. Specifically, this is achieved through 3D convolutional layers, IN normalization layers, and pooling layers. Each pooling operation yields a feature map of the next size. Before each pooling operation, there are 2, 3, 4, and 6 layers of 3D convolutions, and each 3D convolutional layer has a jumper structure followed by an IN normalization layer. The four feature maps generated in the end have 12, 24, 24, and 24 channel dimensions, respectively.

[0064] Step (2) Preprocessing

[0065] The input for this step first upsamples the displacement field φ and uncertainty σ from the previous stage using trilinear interpolation to obtain the same size as the feature map; if it is the first stage, the displacement field and uncertainty are initialized to zero vectors. Then, the floating features of the current stage obtained in step (1) are deformed (e.g., Figure 2 The module shown in the figure can be used to sample the upsampled displacement field φ as coordinates to obtain the floating characteristics after deformation in the current stage, which is simply referred to as deformation characteristics.

[0066] Step (3) Generate displacement space

[0067] like Figure 2 The specific steps for using the CE module are as follows: Figure 3 As shown. Specifically, 27 points are selected in the displacement space, corresponding to values ​​{-1, 0, 1} on each axis in three-dimensional space, forming a total of 3... 3 =27-point sampling space. The displacement fine-tuning value is obtained by passing the floating features through a series of self-convolutional structures, which consist of three hourglass-shaped convolutional modules stacked together. Each hourglass convolutional module has four convolutional layers. The first three convolutional layers have strides of 1, 2, and 1, respectively, while the fourth layer is a deconvolution with a stride of 2, connected by a skip structure. After passing through these four layers, the feature scale remains unchanged. The displacement cardinality is generated by linearly changing the uncertainty. Finally, the displacement fine-tuning value and the displacement cardinality are summed to obtain the displacement space.

[0068] Step (4) Estimation of multiple positions

[0069] like Figure 2 The DE module in the middle, the specific steps are as follows: Figure 4 As shown. First, the deformation features are resampled according to the displacement space. For example, when the number of spatial points is 27 and the number of channels is 24, the 24*H*W*D vector becomes 27*24*H*W*D. Then, the two types of features are grouped. For example, the 24-channel dimension is divided into 4 groups, each with 6 dimensions. After obtaining similar cost volumes, cost aggregation is performed through a K-head convolutional network. The K-head convolutional network is defined as a three-layer hourglass convolutional module (similar to...). Figure 3 The network structure is as follows: (the upper part of the network structure), and the last hourglass convolution layer contains a concatenated connection. Additionally, the number of groups in the convolutional layer needs to be defined as K for each K-head. Then, the features generated by all different heads are normalized to obtain the probability vector corresponding to the displacement space. Finally, the displacement space and the probability vector are multiplied by the inner product to obtain the fused displacement field. Finally, the second-order central moment analysis of the displacement field is performed to obtain the uncertainty.

[0070] Once a displacement field and uncertainty of the same scale as the floating feature map are generated, the stage is considered complete and the next stage will begin. Jump to step (2); if all stages are completed, proceed to step (5).

[0071] Step (5) Construct the supervised loss function

[0072] like Figure 2 As shown, after obtaining the displacement field with the same scale as the original image, the floating image is deformed and a similarity loss function is constructed with the fixed image. At the same time, a smoothness loss function is constructed, and the two are summed to obtain the final loss function.

[0073] Step (6) Training and Validation Methods

[0074] The network optimizer uses Adam with an exponentially decaying learning rate curve. The initial learning rate is set to 0.0005, and it is reduced to 0.95 every 10,000 iterations for a total of 100,000 iterations. Moving averages are used to update the model parameters. If the dataset is insufficient, K-fold cross-validation is used.

[0075] Step (7) Evaluation Indicators

[0076] Since the lung parenchyma often occupies a large number of voxels in CT images, TRE between landmarks is used to assess the registration quality within the lung, that is, the degree of similarity after registration.

[0077]

[0078] Where x f ,x m These represent the sets of landmarks corresponding to the floating and fixed images, respectively. These are typically manually marked or machine-generated reference points with anatomical value. In addition, there is a smoothness index, which primarily examines the folded voxels of the deformed image. The location of folded voxels in the deformation field φ can be determined using the Jacobian determinant. Define a voxel with coordinates (i,j,k), and the Jacobian determinant value is:

[0079]

[0080] The smoothness index is generally defined as the number of voxels whose Jacobian determinant is negative.

[0081]

[0082] The displacement field obtained in this embodiment is as follows: Figure 5As shown in Figures 6(a) and 6(b), the registration results are compared. Evaluation using an evaluation index was conducted. Under the condition of ensuring smoothness (i.e., the number of negative values ​​in the Jacobian determinant is 0), the deformed image showed a high degree of similarity to the fixed image, indicating a good registration result and demonstrating the effectiveness of the method of this invention.

Claims

1. A non-rigid registration method for medical images, characterized in that, The steps are as follows: Step (1) Pyramid Feature Encoding The input to this step is a floating image and a fixed image, which are then fed into a pyramid feature encoding network with shared parameters to obtain floating or fixed feature maps with different sampling ratios. The pyramid feature encoding network generates feature maps with different sampling ratios, i.e., a feature pyramid. The main component of the feature encoding network is a 3D convolutional layer, which is stacked with the residual block structure of ResNet as the basic unit, and different sampling ratios are obtained by setting the stride of the convolutional layer. The next step takes the feature map with each ratio as a stage, and starts with the highest sampling ratio, i.e. the smallest size feature map, to generate larger and more precise fields stage by stage. Step (2) Preprocessing This step first inputs the displacement field from the previous stage. Uncertainty Upsampling is performed to obtain the same size as the feature map; if it is the first stage, the displacement field and uncertainty are initialized to zero vectors; then the floating features of the current stage obtained in step (1) are deformed, which only requires the upsampled displacement field. This process can then be completed, yielding the floating characteristics after deformation at the current stage, which is simply referred to as deformation. feature: ； Where x, y, and z represent the voxel positions, corresponding to the three coordinate axes of the 3D image; and These are floating features and deformation features, respectively; Displacement field It contains three 3D matrices, each with the same size as the floating feature, representing the relationship matrix for transforming the floating feature into the fixed feature. This is a deformation operation; the deformation operation performed on different channels is the same. Step (3) Generate displacement space The input for this step consists of two parts: one is the floating feature map of the current stage, which is fed into a multi-layer convolutional layer to generate shift fine-tuning values, specifically in the form of a four-dimensional tensor. ； in represents the displacement fine-tuning value, and t represents the number of points in the displacement space; c represents the number of channels in the feature map; Convs represents a series of self-network structures, consisting of multiple layers of 3D convolution, instance normalization, and non-linear activation functions; Another part of the input for this step is the uncertainty from the previous stage, which is transformed into a displacement basis value through a linear transformation; ； in Indicates the base value of the displacement; Indicates uncertainty, and The parameter represents the linear transformation, and space represents a cubic sampling space in three-dimensional space with t points; before performing the second multiplication, the tensor needs to be expanded to the same dimension. Summing the displacement base value and the displacement fine-tuning value yields a coordinate vector, which is the final discrete displacement space: = + ； Where Disp represents the final displacement space, specifically a coordinate vector; Step (4) Estimation of multiple positions The input for this step consists of two parts. The first part is the deformation features of the current stage obtained through step (2) and the fixed feature map obtained through step (1). The other part is the adaptive coordinates, i.e., the displacement space, generated in step (3). The displacement field is obtained through multi-head correlation, cost aggregation, normalization, and displacement fusion, and the uncertainty of the field is estimated simultaneously. Specifically, the deformation characteristics are first... Sampling is performed based on adaptive coordinates, and this involves t rounds of sampling. Therefore, the generated features will have an additional dimension of t. ； in This represents the multi-displacement features obtained from sampling; Multi-displacement features, fixed features Dividing the channel into k heads is essentially grouping the channels. operate: ； ； Then Perform relevant calculations along the dimension to obtain the k-head similarity cost body. : ； in Indicates inner product operation; Then, cost aggregation is performed using a k-head convolutional network to finally obtain the cost volume. : Cost ] Where K_Convs represents a series of self-network structures, consisting of multiple layers of 3D convolutions, instance normalization, and non-linear activation functions, and the number of groups in the convolutional layers is k; Then, normalize the displacements generated by all different heads to obtain the probability vector Pro corresponding to the displacement space: Pro ] ; in Softmax represents the normalization function, which is performed in two dimensions, k and t. Finally, the displacement space With probability vector Perform the inner product to obtain the fused displacement field. : ； in This indicates an inner product operation performed along two dimensions, k and t, with the tensor expanded beforehand. Due to displacement field Essentially, it's the result of the first moment, while the root result of the second central moment yields the uncertainty. The specific calculation process is as follows: ； Once a displacement field and uncertainty of the same scale as the floating feature map are generated, the stage is considered complete and the next stage will begin. Jump to step (2); if all stages are completed, proceed to step (5). Step (5) Construct the unsupervised loss function Obtain the displacement field at the same scale as the original image. Then, for the floating image The deformation process is performed, and the resulting deformed image is the final registration result. : ； Training a neural network model requires constructing a loss function. An unsupervised loss function is used, and its construction consists of two parts: one is the similarity loss function. Secondly, the smoothness loss function of the displacement field. The two are summed to obtain the final loss function. This is the balance parameter between the two, also known as the regularization coefficient: 。

Images