A method for registering and fusing a low-dose CT image with an x-ray thermoacoustic image

Abstract

The application discloses a method for registration and fusion of low-dose CT images and X-ray thermoacoustic images, and belongs to the technical field of medical image processing. The method comprises the following steps: collecting low-dose CT image and X-ray thermoacoustic image data, and respectively pre-processing the low-dose CT image and the X-ray thermoacoustic image data; constructing a registration model comprising selective large-step patch sampling and an initial momentum field patch prediction network, constructing a motion image in a to-be-registered image based on the registration model, and deforming the motion image; and inputting the deformed motion image and a fixed image into an image fusion network to obtain a fusion image. The selective large-step patch sampling takes patches as input data, thereby reducing the calculation amount. The initial momentum prediction network is trained by using a traditional iterative method, thereby greatly improving the speed and accuracy of registration and fusion. The low-dose CT image and the X-ray thermoacoustic image are registered and fused, thereby effectively improving the imaging quality.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of medical image processing, and particularly relates to a method for registration and fusion of low-dose CT images and X-ray thermoacoustic images. BACKGROUND

[0002] Image fusion is a technology of combining information of multiple images to form an image by using specific fusion rules. The fused image has more superior performance and higher image quality than each sub-image, and has important significance in the field of medical images. Common medical images include low-dose CT images and X-ray thermoacoustic images. Low-dose CT is a main screening method for lung diseases, and is more widely used in clinical practice due to low radiation dose. CT imaging is an image formed by using different absorption levels of different tissue structures to X-rays. Low-dose CT images can not only achieve the diagnostic effect of conventional CT images, but also reduce the radiation dose of patients. X-ray thermoacoustic imaging is an imaging method in which ultrasonic waves are generated by thermal expansion of X-rays when they are transmitted through tissues, and the ultrasonic waves are propagated spherically in all directions from the X-ray absorber, and the X-ray thermoacoustic image is reconstructed by back projection after being detected by a sensor. Studies have shown that low-dose CT images have certain advantages in imaging depth compared to X-ray thermoacoustic images, while X-ray thermoacoustic images have advantages of high contrast and high ultrasonic resolution compared to low-dose CT images. Therefore, registration and fusion of low-dose CT images and X-ray thermoacoustic images can combine the advantages of both, thereby achieving higher image quality and providing greater help for medical diagnosis.

[0003] Image fusion generally includes image preprocessing, image registration and image fusion. The key factor affecting the quality of medical image fusion is image registration. Image registration refers to finding one or a series of spatial transformations to make one image consistent with the corresponding points on another image, that is, finding a deformation field to minimize the energy function. In the prior art, common medical image registration methods include traditional iterative methods and deep learning methods. The main idea of the traditional iterative method is: for a certain pair of images, given the energy function and the initial deformation field parameters, different optimization strategies are used to iteratively find better deformation field parameters, so that the minimum value is gradually approached. The main idea of deep learning registration is: first, the neural network layer is used to automatically extract and learn the features of the input image, and then the mapping relationship between the input image pair and the deformation field is obtained, and the deformation field of the unknown image pair is predicted, for example, the Chinese invention patent with document number CN 112102385 B discloses a multi-modal liver magnetic resonance image registration system based on deep learning.

[0004] However, it is difficult to ensure that the deformation field is differentiable homeomorphism by using deep learning to predict the deformation field, that is, it is difficult to ensure that the obtained deformation field has a high degree of smoothness; at the same time, in some fuzzy areas, the neural network model is difficult to predict an accurate deformation field, so the existing deep learning registration method has the disadvantages of poor robustness and low precision. Although the traditional iterative registration method has high precision, it is time-consuming and cannot cope with large-scale image registration. In addition, since most of the current registration research mainly involves non-rigid registration, it is limited by the requirement that the deformation field must remain smooth. Whether it is deep learning or traditional iterative method, most of the current researches take the entire image as input, which leads to too many parameters to be optimized and too large a calculation amount. If it is three-dimensional image registration, the calculation amount will be even larger. SUMMARY

[0005] In view of the problems in the prior art, the present application provides a method for registration and fusion of low-dose CT images and X-ray thermoacoustic images, which aims to solve the problems of low registration accuracy, large calculation amount and poor robustness of the existing image registration and fusion method during registration and fusion. In order to achieve the above-mentioned purpose, the present application provides the following technical solutions.

[0006] A method for registration and fusion of low-dose CT images and X-ray thermoacoustic images, comprising the following steps:

[0007] Step S1, acquiring low-dose CT image and X-ray thermoacoustic image data, and pre-processing the low-dose CT image and X-ray thermoacoustic image data respectively;

[0008] Step S2, constructing a registration model comprising selective large-step patch sampling and initial momentum field patch prediction network, constructing a motion image in the image to be registered based on the registration model, and deforming the motion image;

[0009] Step S3, inputting the deformed motion image and the fixed image into an image fusion network to obtain a fused image.

[0010] Further, the step of acquiring low-dose CT image and X-ray thermoacoustic image data and pre-processing the low-dose CT image and X-ray thermoacoustic image data respectively comprises:

[0011] The pre-processing of the X-ray thermoacoustic image comprises histogram equalization and image normalization processing;

[0012] The pre-processing of the low-dose CT image comprises histogram equalization, contrast stretching and image normalization processing.

[0013] Further, the constructing comprises a registration model of the selective large-step patch sampling and the initial momentum field patch prediction network, and the step of constructing the motion image in the image to be registered and deforming the motion image based on the registration model comprises:

[0014] In step S21, the preprocessed low-dose CT image and the X-ray thermoacoustic image are synchronously sampled multiple times by using a selective large-step patch sampling method to obtain a plurality of pairs of image patches.

[0015] In step S22, the patch data is input into the initial momentum field patch prediction network to obtain a predicted initial momentum field patch, and the initial momentum field patch is integrated to obtain a complete initial momentum field.

[0016] In step S23, the complete initial momentum field is converted into a momentum to obtain a complete smooth deformation field.

[0017] In step S24, the motion image is deformed by using the complete smooth deformation field.

[0018] Further, the step of synchronously sampling multiple times the preprocessed low-dose CT image and the X-ray thermoacoustic image by using the selective large-step patch sampling method to obtain a plurality of pairs of image patches comprises:

[0019] In step S211, the preprocessed low-dose CT image and the X-ray thermoacoustic image are respectively numbered, and the patches at the same position have the same number.

[0020] In step S212, all patch numbers are sequentially synchronously traversed, and the accepted patch numbers are recorded.

[0021] In step S213, based on the accepted patch numbers, the three-dimensional space coordinates of the patches are calculated and divided to obtain a plurality of pairs of image patches in the form of three-dimensional arrays.

[0022] Further, the step of momentum conversion comprises:

[0023] In step S231, the geodesic equation is integrated and smoothed based on the complete initial momentum field to obtain a complete smooth velocity field of the motion image to the fixed image; wherein the geodesic equation is:

[0024] m t +Dmv-Dvm=0 (1)

[0025] Wherein, m represents a momentum field; m t represents the derivative of the momentum field with respect to time t, 0≤t≤1; v represents a smooth velocity field; D is a Jacobian operator.

[0026] Step S232, based on the complete smooth velocity field, the flow equation is integrated to obtain the complete smooth deformation field; wherein, the flow equation is:

[0027]

[0028] Wherein, Φ -1 Indicates the inverse of the smooth deformation field Φ; Indicates Φ -1 The derivative of time t.

[0029] Further, the training of the initial momentum field patch prediction network is based on supervised learning; the true value of the supervised learning is the initial momentum field patch true value obtained by using the improved large deformation differential homeomorphism mapping method to obtain the complete initial momentum field true value and sampling.

[0030] Further, the step of the improved large deformation differential homeomorphism mapping method to obtain the complete initial momentum field true value comprises:

[0031] Given a pair of input images, including a moving image and a fixed image, define the initial value of the complete initial momentum field;

[0032] The complete smooth deformation field is calculated by using momentum conversion, the moving image is deformed by using the complete smooth deformation field, and the similarity of the deformed moving image and the fixed image is calculated;

[0033] Determine whether the similarity reaches the highest, if so, output the complete initial momentum field true value;If not, update the initial momentum after calculating the gradient, continue to optimize until the similarity is the highest, and the complete initial momentum field true value is obtained.

[0034] Further, the image fusion network is a convolutional neural network.

[0035] The beneficial effects of the present application are:

[0036] The application discloses a method for registration and fusion of low-dose CT images and X-ray thermoacoustic images, and belongs to the technical field of medical image processing. BRIEF DESCRIPTION OF DRAWINGS

[0037] Figure 1 is a flowchart of the method for registration and fusion of low-dose CT images and X-ray thermoacoustic images according to the application;

[0038] Figure 2 is a framework diagram of the method for registration and fusion of low-dose CT images and X-ray thermoacoustic images according to the application;

[0039] Figure 3 is a schematic diagram of the initial momentum field patch prediction network according to the application;

[0040] Figure 4 is a schematic diagram of the image fusion convolutional neural network according to the application;

[0041] Figure 5 is a flowchart of the improved large deformation differential homeomorphism metric mapping method used in the application, and the improved large deformation differential homeomorphism metric mapping method refers to the LDDMM-shooting method. DETAILED DESCRIPTION

[0042] The application will be further described in detail below in combination with the drawings and specific embodiments, but the application is not limited to the following embodiments.

[0043] Embodiment one:

[0044] See the attached Figures 1-5 A method for registration and fusion of low-dose CT images and X-ray thermoacoustic images comprises the following steps:

[0045] Step S1, low-dose CT image and X-ray thermoacoustic image data are collected, and low-dose CT image and X-ray thermoacoustic image data are respectively preprocessed;

[0046] Step S2, a registration model including a selective large-step patch sampling and an initial momentum field patch prediction network is constructed, a motion image in the image to be registered is constructed based on the registration model, and the motion image is deformed;

[0047] Step S3, the deformed motion image and the fixed image are input into an image fusion network to obtain a fused image.

[0048] Specifically, image fusion is to synthesize two or more images into a new image by using a specific algorithm. The fusion result can utilize the correlation in time and space and the complementarity of information of two or more images, and make the fused image more comprehensive and clear in describing the scene, thereby being more beneficial to the recognition of human eyes and the automatic detection of machines. Medical image registration fusion has high theoretical and application value. In view of the problems of low contrast and low resolution in the depth direction of low-dose CT images, the X-ray thermoacoustic imaging technology has the advantages of high contrast and high ultrasonic resolution, and the fusion of X-ray thermoacoustic images and low-dose CT images is expected to improve the imaging quality of low-dose CT and provide greater help for medical diagnosis. The method flow of the low-dose CT image and X-ray thermoacoustic image registration fusion of the present application is shown in Figure 1 The specific steps include:

[0049] Step S1, low-dose CT image and X-ray thermoacoustic image data are collected, and the low-dose CT image and X-ray thermoacoustic image data are preprocessed respectively. First, low-dose CT image and X-ray thermoacoustic image data are collected respectively. From the imaging mode, the image registration fusion of the present application belongs to multi-modal registration fusion. Then, the low-dose CT image and X-ray thermoacoustic image data are preprocessed respectively. The main purpose of image preprocessing is to eliminate irrelevant information in the image, restore useful real information, enhance the detectability of relevant information and simplify the data as much as possible to avoid the influence of irrelevant factors on the subsequent model. For example, according to the characteristics of X-ray thermoacoustic image and low-dose CT image, the X-ray thermoacoustic image can be preprocessed by histogram equalization, image normalization, etc., and the low-dose CT image can be preprocessed by histogram equalization, contrast stretching, image normalization, etc.

[0050] Step S2, a registration model including selective large step patch sampling and initial momentum field patch prediction network is constructed, a motion image in the image to be registered is constructed based on the registration model, and the motion image is deformed. Image registration is a prerequisite and key of image fusion, and the accuracy of image registration directly determines the quality of the fusion result. The registration model of the application includes two parts of selective large step patch sampling and initial momentum field patch prediction network. First, a patch-level selective large step sampling method is used to sample two images synchronously multiple times to obtain multiple pairs of smaller size image patches as input data for the subsequent initial momentum field patch prediction network. Compared with the registration method of taking the whole image as input, the application selectively extracts the patches of the image and takes the patches as input data, avoiding the problem of large calculation caused by the large size and dimension of the parameters caused by taking the whole image as input data. Then the patch data is input into the initial momentum field patch prediction network, which will predict a number of small size initial momentum field patches, which are then integrated into a complete initial momentum field. Based on the complete initial momentum field, a smooth deformation field is obtained, and the motion image is deformed using the deformation field to realize the registration of the X-ray thermoacoustic image and the low dose CT image. In the registration problem, one image is regarded as a fixed image and the other image is regarded as a motion image, and the purpose of registration is to align the motion image with the fixed image. When constructing the motion image and the fixed image to be registered, the X-ray thermoacoustic image can be taken as the motion image and the low dose CT image as the fixed image, or the low dose CT image can be taken as the motion image and the X-ray thermoacoustic image as the fixed image.

[0051] Step S3, input the deformed motion image and the fixed image into the image fusion network to obtain the fused image. The deformed motion image and the fixed image are taken as the input of the image fusion network, and the output of the fusion network is the fused image. The image fusion network can adopt a convolutional neural network. The application makes up for the lack of registration research of low dose CT images and X-ray thermoacoustic images, and solves the problem of large calculation and slow speed of existing multi-modal image registration.

[0052] Embodiment two:

[0053] See the attached Figures 1-5On the basis of the embodiment one, the step of collecting low-dose CT image and X-ray thermoacoustic image data and respectively pre-processing the low-dose CT image and the X-ray thermoacoustic image data comprises: the pre-processing of the X-ray thermoacoustic image comprises histogram equalization and image normalization processing; the pre-processing of the low-dose CT image comprises histogram equalization, contrast stretching and image normalization processing. Specifically, the low-dose CT image has the problems of poor contrast and low resolution in the depth direction, so the pre-processing of the low-dose CT image comprises histogram equalization, contrast stretching and image normalization processing, and the pre-processing of the X-ray thermoacoustic image comprises histogram equalization and image normalization processing. Histogram equalization is a method of adjusting the contrast of an image by using an image histogram, by which the brightness can be better distributed on the histogram, and the local contrast is enhanced without affecting the overall contrast. And image normalization is a process of performing a series of standard processing transformations on an image to transform it into a fixed standard form, and after the same parameter image normalization processing, a standard image of the same form can be obtained, which is more conducive to subsequent registration fusion.

[0054] The registration model comprises a selective large-step patch sampling and an initial momentum field patch prediction network, and the step of constructing a motion image in the image to be registered based on the registration model and deforming the motion image comprises:

[0055] In step S21, the pre-processed low-dose CT image and X-ray thermoacoustic image are synchronously sampled multiple times by using a selective large-step patch sampling method to obtain multiple pairs of image patches.

[0056] In step S22, the patch data is input into the initial momentum field patch prediction network to obtain predicted initial momentum field patches, and the initial momentum field patches are integrated to obtain a complete initial momentum field.

[0057] In step S23, the complete initial momentum field is subjected to momentum conversion to obtain a complete smooth deformation field.

[0058] In step S24, the motion image is deformed by using the complete smooth deformation field.

[0059] Specifically, the registration model is used to realize the registration of the low-dose CT image and the X-ray thermoacoustic image, and the specific steps comprise:

[0060] Step S21 involves using a selective large-step patch sampling method to simultaneously sample the preprocessed low-dose CT image and X-ray thermoacoustic image multiple times, resulting in multiple pairs of image patches. Specifically, the preprocessed low-dose CT image and X-ray thermoacoustic image are sampled through a sliding window, synchronously traversed sequentially. When the sampling condition is met, the image patch at the current position is accepted and recorded, and then the traversal continues until all images have been traversed, ultimately resulting in multiple pairs of image patches.

[0061] Step S22: Input the patch data into the initial momentum field patch prediction network to obtain the predicted initial momentum field patch, and integrate the initial momentum field patches to obtain the complete initial momentum field. The initial momentum field patch prediction network is used to predict the corresponding initial momentum field patch for the input image patch data. The process includes: First, defining the complete initial momentum field m0, and setting the element values ​​of the complete initial momentum field m0, i.e., the gray values, to zero, and making its size the same as the moving image and the stationary image; then, inputting the image patch obtained by the selective large stride patch sampling method into the initial momentum field patch prediction network for prediction, the initial momentum field patch prediction network will predict several small-sized initial momentum field patches. Finally, calculate the three-dimensional spatial coordinates of the initial momentum field patch based on the image size and the predicted initial momentum field patch, and set the patch gray value at the corresponding spatial position in m0 to be equal to the gray value of the predicted initial momentum field patch, completing the integration operation of the predicted initial momentum field patch to obtain the complete initial momentum field.

[0062] The initial momentum field patch prediction network can be a pre-trained prediction network or a retrained one. For example, the initial momentum field patch prediction network is as follows: Figure 3 As shown, the initial momentum field patch prediction network consists of an encoder and a decoder. The encoder comprises two parallel sub-encoders with dual channels, learning features independently from moving and stationary image patches, respectively. Each encoder contains three 3×3×3 convolutional layers and two PReLU layers. Compared to ReLU, PReLU avoids zero gradients on negative inputs, effectively improving network performance. It also includes another 2×2×2 convolutional layer and PReLU with a stride of 2. The 2-stride convolutional layer essentially performs pooling, significantly reducing the size of the output patch. The decoder consists of three parallel sub-decoders, each with a structure reversed from the encoder, doubling the number of features. Each sub-decoder receives input from two encoder branches; the three sub-decoders are the x-decoder, y-decoder, and z-decoder. Through transposed convolution operations, the x-decoder, y-decoder, and z-decoder predict the initial momentum field patches in the x, y, and z directions, respectively. The initial momentum field patches in the three directions are vector synthesized to obtain the final predicted value of the initial momentum field patch.

[0063] Step S23, momentum conversion is performed on the complete initial momentum field to obtain a complete smooth deformation field. After obtaining the complete initial momentum field, momentum conversion needs to be performed on the complete initial momentum field to obtain a complete smooth deformation field, which is used for deformation processing of the motion image.

[0064] Step S24, the motion image is deformed using the complete smooth deformation field, and the deformed motion image is the image registered with the fixed image. The present application avoids the problem of large calculation amount caused by the large size and dimension of the parameters due to the input of the entire image as input data by selectively extracting patches of the image and taking the patches as input data. Meanwhile, the complete initial momentum field has the advantage of discreteness, which makes it easy to predict the initial momentum field patch using a neural network.

[0065] Embodiment three:

[0066] See the attached Figures 1-5 On the basis of embodiment two, the step of synchronously sampling the preprocessed low-dose CT image and the X-ray thermoacoustic image multiple times to obtain multiple pairs of image patches includes:

[0067] Step S211, the preprocessed low-dose CT image and the X-ray thermoacoustic image are respectively numbered, and the patches at the same position have the same number;

[0068] Step S212, all patch numbers are sequentially synchronously traversed, and the accepted patch numbers are recorded;

[0069] Step S213, based on the accepted patch numbers, the three-dimensional space coordinates of the patches are calculated and divided to obtain a plurality of pairs of image patches in the form of three-dimensional arrays.

[0070] Specifically, the step of synchronously sampling the preprocessed low-dose CT image and the X-ray thermoacoustic image multiple times to obtain multiple pairs of image patches includes:

[0071] Step S211, the preprocessed low-dose CT image and the X-ray thermoacoustic image are respectively numbered, and the patches at the same position have the same number; for the motion image, the numbers of all patches are calculated according to the size of the image, the sliding window or the three-dimensional space coordinates of the patches, and the same operation is performed on the fixed image. After the numbering operation is completed, the patch numbers at the same positions of the motion image and the fixed image should be the same. The formula for calculating the patch number using the image size and the three-dimensional space coordinates is as follows:

[0072] number=pos_x*width*height+pos_y*height+pos_z

[0073] wherein, number is patch number, pos_x, pos_y and pos_z are three-dimensional space coordinates of the sliding window or patch respectively, width is the width of the image, and height is the height of the image.

[0074] In step S212, all patch numbers are sequentially synchronized and the accepted patch numbers are recorded. Starting from the first number, the patch numbers of the moving image and the fixed image are synchronized, if the patches of the moving image and the fixed image corresponding to a certain number are located in the background area with uniform gray value in the image, then the two patches are meaningless, so the two patches are rejected to be input into the initial momentum field patch prediction network; otherwise, the two patches are accepted to be input into the initial momentum field patch prediction network, and then the following traversal is continued, and the rejected patch numbers and the accepted patch numbers are recorded.

[0075] In step S213, based on the accepted patch numbers, the three-dimensional space coordinates of the patches are calculated and divided to obtain a plurality of image patches in the form of three-dimensional arrays. For the moving image, according to the image size and the accepted patch numbers in the previous step, the three-dimensional space coordinates of the patches can be calculated, the three-dimensional space coordinates are taken as array indexes, the three-dimensional array of the moving image is divided to obtain a plurality of patches in the form of three-dimensional arrays, and the same operation is performed on the fixed image. The formula for calculating the three-dimensional space coordinates of the patches using the image size and the patch number is as follows:

[0076] pos_x = number / (width*height)

[0077] pos_yz = number%(width*height)

[0078] pos_y = pos_yz / height

[0079] pos_z = pos_yz%height

[0080] wherein, number is patch number, pos_x, pos_y, pos_z are three-dimensional space coordinates of the sliding window or patch respectively, pos_yz is a temporary variable, width is the width of the image, and height is the height of the image.

[0081] The selective large stride sampling method significantly improves the registration speed through two ways of selective sliding window and large stride. Firstly, based on the spatial displacement invariance principle of pooling / depooling, selective sampling can be realized by using sliding window to obtain patches. When the patch is located in the above-mentioned fuzzy area, i.e. the background area, the patch is rejected to be input into the initial momentum field patch prediction network, and the prediction work is directly skipped to the subsequent sampling, so as to reduce the number of patches as much as possible. Secondly, according to the existing improved large deformation differential homeomorphism metric mapping theory, i.e. LDDMM-shooting theory, the initial momentum field is discrete, and the continuity constraint is not considered when predicting the initial momentum field patch, both inside and outside the image. Therefore, a larger moving stride can be set for the sliding window, and the large stride operation can minimize the number of patches input into the initial momentum field patch prediction network.

[0082] The step of momentum conversion comprises:

[0083] In step S231, the geodesic equation is integrated and smoothed based on the complete initial momentum field to obtain a complete smoothed velocity field of the moving image to the fixed image; wherein the geodesic equation is:

[0084] m t +Dmv-Dvm=0 (1)

[0085] Wherein, m represents the momentum field; m t represents the derivative of the momentum field with respect to time t, 0≤t≤1; v represents the smoothed velocity field; D is the Jacobian operator; for example, the geodesic equation of momentum and velocity can be integrated first, and then combined with the smoothing operator K to convert the complete initial momentum field into the complete smoothed velocity field.

[0086] In step S232, the flow equation is integrated based on the complete smoothed velocity field to obtain a complete smoothed deformation field; wherein the flow equation is:

[0087]

[0088] Wherein, Φ -1 represents the inverse of the smoothed deformation field Φ; represents the derivative of Φ -1 with respect to time t.

[0089] Embodiment four:

[0090] See the attached Figures 1-5 On the basis of embodiment three, the training of the initial momentum field patch prediction network is based on supervised learning; the true value of the supervised learning is the initial momentum field patch true value obtained by sampling the complete initial momentum field true value obtained by using the improved large deformation differential homeomorphism metric mapping method.

[0091] Specifically, the initial momentum patch prediction network training will affect the final prediction results, the better the training, the more accurate the prediction results. The initial momentum patch prediction network based on supervised learning training needs an initial momentum patch true value as a reference. This value can be obtained in advance, for example, the improved large deformation diffeomorphic metric mapping can be used to obtain the complete initial momentum field, and then the initial momentum patch true value can be obtained by sampling. The initial momentum patch true value can be used to guide the training of the initial momentum patch prediction network.

[0092] The improved large deformation diffeomorphic metric mapping method for obtaining a complete initial momentum field true value comprises the following steps:

[0093] Given a pair of input images, including a moving image and a fixed image, define the initial value of the complete initial momentum field; calculate the complete smooth deformation field using momentum conversion, deform the moving image using the complete smooth deformation field, and calculate the similarity of the deformed moving image and fixed image; determine whether the similarity reaches the highest, if so, output the complete initial momentum field true value; if not, update the initial momentum after calculating the gradient, continue to optimize until the similarity is the highest, and obtain the complete initial momentum field true value.

[0094] Specifically, the large deformation diffeomorphic metric mapping method is a traditional iterative fluid image registration method, and its main idea is to regard the image to be registered as a fluid model that can move, and each position of the fluid has a velocity with different direction and size. The large deformation diffeomorphic metric mapping method deforms the image iteratively by constraining the velocity or displacement of the voxel through a regularization term, and stops iteration when the similarity of the moving image and the reference image is optimal. In the parameter optimization process of the large deformation diffeomorphic metric mapping method, the deformation field or the velocity field is required to be continuous, which will bring many limitations. Therefore, there are variants and improvements of the large deformation diffeomorphic metric mapping method: the LDDMM-shooting method of discrete momentum parameters. The parameters of the LDDMM-shooting energy function are discrete complete initial momentum fields, and the continuity constraint of the image edge does not need to be considered. Its formula is as follows:

[0095] E(m0) = <m0, K m0> + Sim(M, T, Φ), s.t.

[0096] m t +ad *v m=0

[0097] m(0)=m o

[0098]

[0099] Φ -1(0) = id

[0100] m = Lv

[0101] where Sim() denotes a similarity measure, m0is the complete initial momentum field, K is a smoothing operator, M and T denote the moving and fixed images, respectively; Φ denotes a smooth deformation field, ad * is an operator, ad *v m = Dmv - Dvm, m is a momentum field, v is a smooth velocity field, D is a Jacobian operator, m t denotes the derivative of the momentum field with respect to time t, 0≤t≤1; m(0) denotes the momentum field at time 0; Φ -1 denotes the inverse of the smooth deformation field Φ; denotes the derivative of Φ -1 with respect to time t; Φ -1 (0) denotes the inverse of the smooth deformation field at time 0, id denotes the identity matrix, L is the inverse of the smoothing operator K, m is the dual of v, m and v are connected by a definite self-adjoint differential smoothing operator K, i.e., v = Kmand m = Lv. Given the complete initial momentum field m0, the smooth deformation field Φ(x, t) at any time can be calculated. The parameters of the LDDMM-shooting energy function are the discrete complete initial momentum field, and the continuity constraint of the image edge does not need to be considered.

[0102] As shown in Figure 5 , the steps of the LDDMM-shooting method for obtaining the real value of the complete initial momentum field include: given a pair of input images, including a moving image and a fixed image, defining an initial value of the complete initial momentum field; calculating a complete smooth deformation field by using momentum conversion, deforming the moving image by using the complete smooth deformation field, and calculating the similarity of the deformed moving image and the fixed image; judging whether the similarity reaches the highest, if yes, outputting the real value of the complete initial momentum field; if not, updating the initial momentum after calculating the gradient, continuing to optimize until the similarity reaches the highest, and obtaining the real value of the complete initial momentum field.

[0103] The image fusion network is a convolutional neural network. As shown in Figure 4 , the registered images can be fused by the convolutional neural network. The specific process is as follows: first, stack the single-channel medical images into three-channel images as network input; then use convolutional layers, ReLU activation functions and batch normalization to perform feature extraction operations on the images, to obtain multi-channel image features; then use the element maximum fusion rule to perform information fusion processing on the image features; finally, generate a fused image through image reconstruction operation. In order to reduce the loss of image information, the feature extraction performed by the neural network can use a non-downsampling process.

[0104] The application can utilize the no-reference image quality evaluation index, specifically including average gradient, edge intensity, entropy and standard deviation. Among them, the average gradient and the edge intensity reflect the clarity of the image, the higher the value is, the better; the entropy reflects the information quantity of the image, the higher the value is, the better; and the standard deviation reflects the dispersion degree of the image pixel value and the mean value, the higher the value is, the better.

[0105] The application reduces the calculation amount by selecting a large-step patch sampling and taking the patch as the input data; and greatly improves the speed and accuracy of the registration fusion by using the traditional iterative method to guide the initial momentum prediction network training. The application can effectively solve the problems of low registration accuracy, large calculation amount and poor robustness of the existing image registration fusion method during registration fusion; and effectively improves the imaging quality by registration fusion of the low-dose CT image and the X-ray thermoacoustic image.

[0106] The above description is only the preferred embodiment of the application, and does not limit the patent scope of the application, and any equivalent structure or equivalent process transformation, or direct or indirect application in other related technical fields by using the content of the specification and drawings of the application, are also included in the patent protection scope of the application.

Claims

1. A method of low dose CT image and X-ray thermoacoustic image registration fusion, characterized in that, The method comprises the following steps: Step S1, collecting low-dose CT image and X-ray thermoacoustic image data, and respectively preprocessing the low-dose CT image and the X-ray thermoacoustic image data; Step S2, constructing a registration model comprising a selective large-step patch sampling method and an initial momentum field patch prediction network, constructing a motion image in the image to be registered based on the registration model, and deforming the motion image; comprising: Step S21, using the selective large-step patch sampling method to synchronously sample the preprocessed low-dose CT image and X-ray thermoacoustic image multiple times to obtain multiple pairs of image patches; Step S22, inputting the patch data into the initial momentum field patch prediction network to obtain the predicted initial momentum field patch, and integrating the initial momentum field patch to obtain the complete initial momentum field; Step S23, performing momentum conversion on the complete initial momentum field to obtain a complete smooth deformation field; Step S24, deforming the motion image using the complete smooth deformation field; Step S3, inputting the deformed motion image and the fixed image into an image fusion network to obtain a fused image.

2. The method of fusing a low-dose CT image with an X-ray thermoacoustic image registration according to claim 1, characterized in that, The step of collecting low-dose CT image and X-ray thermoacoustic image data, and respectively preprocessing the low-dose CT image and the X-ray thermoacoustic image data comprises: The preprocessing of the X-ray thermoacoustic image comprises histogram equalization and image normalization processing; The preprocessing of the low-dose CT image comprises histogram equalization, contrast stretching and image normalization processing.

3. The method of fusing a low dose CT image with an X-ray thermoacoustic image registration according to claim 1, wherein, The step of using the selective large-step patch sampling method to synchronously sample the preprocessed low-dose CT image and X-ray thermoacoustic image multiple times to obtain multiple pairs of image patches comprises: Step S211, numbering the preprocessed low-dose CT image and X-ray thermoacoustic image respectively, and the patch numbers at the same position are the same; Step S212, synchronously traversing all patch numbers in turn and recording the accepted patch numbers; Step S213, based on the accepted patch numbers, calculating the three-dimensional space coordinates of the patches and dividing them to obtain a plurality of pairs of image patches in the form of three-dimensional arrays.

4. The method of fusing a low dose CT image with an X-ray thermoacoustic image registration according to claim 1, wherein, The step of momentum conversion comprises: Step S231, based on the complete initial momentum field, integrating and smoothing the geodesic equation to obtain a complete smooth velocity field from the motion image to the fixed image; wherein the geodesic equation is: wherein denotes the momentum field; denotes the derivative of the momentum field with respect to time , 0≤t≤1; denotes the smoothed velocity field; is the Jacobian operator; Step S232, based on the complete smooth velocity field, integrating the flow equation to obtain a complete smooth deformation field; wherein the flow equation is: wherein, denotes the inverse of the smooth deformation field denotes the derivative with respect to time .​ 5. The method of registering and fusing a low-dose CT image with an X-ray thermoacoustic image according to claim 1 or 4, characterized in that, The training of the initial momentum field patch prediction network is based on supervised learning; the true value of the supervised learning is the initial momentum field patch true value obtained by using the improved large deformation differential homeomorphism metric mapping method to obtain the complete initial momentum field true value and sampling.

6. The method of fusing a low-dose CT image with an X-ray thermoacoustic image registration according to claim 5, characterized in that, The step of using the improved large deformation differential homeomorphism metric mapping method to obtain the complete initial momentum field true value comprises: Given a pair of input images, including a motion image and a fixed image, define the initial value of the complete initial momentum field; Using the momentum conversion to calculate the complete smooth deformation field, deforming the motion image using the complete smooth deformation field, and calculating the similarity of the deformed motion image and the fixed image; If the similarity reaches the highest, the complete initial momentum field real value is output; if not, the initial momentum is updated after the gradient is calculated, and the optimization is continued until the similarity reaches the highest, and the complete initial momentum field real value is obtained.

7. The method of fusing a low dose CT image with an X-ray thermoacoustic image registration according to claim 1, wherein, The image fusion network is a convolutional neural network.

Images