A medical image demons registration method based on adaptive anisotropic regularization
By using the Demons registration method with adaptive anisotropic regularization, combined with wavelet multi-resolution strategy and mutual information-gray-level squared difference fusion similarity measure, the problem of registration accuracy in complex organ motion is solved, high-precision registration at the edge of the target of interest is achieved, and hardware requirements and training time are reduced.
Patent Information
- Application Number
- CN202310725169.2
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2023-06-19
- Publication Date
- 2025-11-28
- Estimated Expiration
- 2043-06-19
AI Technical Summary
Existing medical image registration methods are not accurate enough when dealing with complex organ motion, and deep learning-based methods have problems such as high hardware requirements, long training time and difficulty in data acquisition. Meanwhile, the Demons method based on the optimization framework has poor registration accuracy in complex images.
The Demons registration method with adaptive anisotropic regularization is adopted, which combines wavelet multi-resolution strategy and mutual information-gray-level squared difference fusion similarity measure. The field is updated by smooth deformation through adaptive anisotropic kernel, and the weights and weights wij are optimized for diffusion to achieve diffusion in the local neighborhood of the image. The diffusion is performed by adaptive anisotropic kernel GAA, and the field is updated by smooth deformation through weighted average adaptive anisotropic kernel, and the weights wij are optimized for diffusion to maximize the diffusion in the local neighborhood of the image.
It improves the accuracy and robustness of medical image registration, especially in handling complex deformations at the edges of targets of interest, while reducing hardware requirements and training time.
Smart Images
Figure CN116703988B_ABST
Abstract
Description
TECHNICAL FIELD
[0001] The application belongs to the technical field of image processing, and specifically relates to a medical image Demons registration method based on adaptive anisotropic regularization. BACKGROUND
[0002] Medical image registration is a process of matching corresponding points in space positions of a floating image and a reference image through a spatial transformation, and is the basis of medical image analysis tasks such as medical image three-dimensional reconstruction and medical image fusion, and has important application value in clinical practices such as image-guided surgery, pathological change prediction, tumor growth detection and radiotherapy planning.
[0003] Existing medical image registration methods mainly include an optimization framework-based registration method and a deep learning-based registration method, wherein the deep learning-based registration method can effectively improve the accuracy of medical image non-rigid registration, but still has the following shortcomings: (1) a black box design leads to poor theoretical explanation; (2) a model design is complex, a training process needs high-cost hardware device support and a long training time; and (3) it is inconvenient to deploy in a resource-limited scene, in particular, image data used in medical clinics belongs to important privacy of patients, and data acquisition is difficult and it is difficult and time-consuming to annotate a label suitable for registration due to the need for privacy protection. The optimization framework-based Demons registration method has a good theoretical basis, does not need data and annotation, has low requirements on hardware resources and is convenient to deploy, and can be used as a complementary method of the deep learning-based registration method in a specific application scene, and is widely used in clinical practices.
[0004] Medical image registration method based on optimization framework regards registration process as an optimization problem, and finds the optimal spatial transformation by iteratively optimizing the similarity measure function of reference image and floating image, and scholars have proposed many non-rigid registration models based on optimization framework, for example: FFD model based on B-spline, viscous fluid model and Demons based on optical flow field, wherein, Demons method based on optical flow field is simple to calculate, has good mathematical theoretical basis, and is widely used in medical image non-rigid registration, Demons model regards the registration process as a diffusion process, and utilizes the deformation force between the reference image and the floating image to drive the floating image to deform to the reference image, so as to realize the matching of spatial information. However, Demons method has problems of being unable to estimate complex large-scale deformation and slow convergence speed. In view of the above problems, Vercauteren et al. improved Demons method by keeping topological consistency, and proposed differential homeomorphism Demons method, however, the method uses isotropic Gaussian filter as a regularization term, and the processing mode of each region of the image is kept consistent in the regularization process, so that the registration accuracy is poor in the motion processing of complex medical image tissues and organs, and at the same time, the gray square difference is used as a similarity measure to measure the similarity of images before and after registration, so that it is difficult to accurately measure the similarity for medical images with large gray difference of organs, and registration error is easy to cause. Therefore, it is of great significance to study a medical image registration method which can accurately register the edge region of the target of interest. SUMMARY
[0005] In view of the defects in the prior art, the purpose of the present application is to provide a medical image Demons registration method based on adaptive anisotropic regularization, which can effectively process complex organ motion and obtain accurate registration results.
[0006] In order to achieve the above purpose, the technical scheme adopted by the present application is as follows:
[0007] A medical image Demons registration method based on adaptive anisotropic regularization comprises the following steps:
[0008] 1) input the reference image F and the floating image M to be registered, and initialize the parameters of the adaptive anisotropic model;
[0009] 2) adopt wavelet-based multi-resolution strategy to decompose the reference image F and the floating image M from coarse to fine, and represent them as F i (i=1, 2,..., N) and M i (i=1, 2,..., N), wherein i represents the number of layers of the reference image F and the floating image M;
[0010] 3) initialize the spatial deformation field between the first layer of the reference image F to be registered and the floating image M
[0011] 4) set the layer number i = 1 of the reference image F and the floating image M to be registered, and start from the first layer to register;
[0012] 5) judge whether i≤N is true: if yes, execute step 6); otherwise, jump to step 17);
[0013] 6) initialize the deformation field of the current layer according to the deformation field of the last layer
[0014] 7) set E i as the energy function of the current layer, E i-1 as the energy function of the last layer, judge whether the condition |E i -E i-1 |<0.001 is true: if yes, execute step 8); otherwise, jump to step 16);
[0015] 8) set as the regularization threshold, judge whether the condition is true: if yes, execute step 9); otherwise, jump to step 10);
[0016] 9) dynamically select the maximum regularization kernel radius r max ;
[0017] 10) dynamically select the minimum regularization kernel radius r min ;
[0018] 11) calculate the deformation update field between the reference image F to be registered and the floating image M is expressed as:
[0019]
[0020] In the formula,
[0021] 12) smooth the deformation update field by using an adaptive anisotropic kernel G AA is expressed as:
[0022] 13) update the deformation field:
[0023] 14) smooth the deformation field:
[0024] 15) update the energy function E by using the mutual information-gray square difference fusion similarity measure and the regularization deformation itemi ;
[0025] 16) output the morphing field of the current layer
[0026] 17) output the optimal morphing field is expressed as:
[0027] Further, the step 1) specifically comprises the following steps:
[0028] 1.1) input the reference image F and the floating image M to be registered;
[0029] 1.2) set the regularization threshold
[0030] 1.3) set the maximum regularization kernel radius r max = 8, and the minimum regularization kernel radius r min = 5.
[0031] Further, the step 12) specifically comprises the following steps:
[0032] 12.1) calculate the output of the adaptive anisotropic regularization model, the regularization output image I o is defined as the local linear model of the guide image g i , expressed as:
[0033]
[0034] wherein, and are local linear coefficients, respectively;
[0035] 12.2) the regularization output image I o is regarded as the weighted average of the input image I i , so an explicit adaptive anisotropic kernel G AA is used to describe the output image I o , expressed as:
[0036]
[0037] wherein, the field window is centered at i and has a radius of r;
[0038] 12.3) calculate the adaptive anisotropic kernel G AA , expressed as:
[0039]
[0040] wherein, is the neighborhood window total number of pixels in the image, μ i and respectively, are the guided image g i The mean and variance in the domain window , i is the center of the domain window , k and j represent two different pixel indices, and ε represents a smoothing coefficient.
[0041] Further, the step 12.1) calculates the output of the adaptive anisotropic regularization model by the process of:
[0042] 12.1.1) describes the anisotropic guided regularization model by a local linear equation, so as to find the optimal a i and b i by minimizing the cost function, which is expressed as:
[0043]
[0044] wherein, i is the pixel index in the image, I i and g i respectively represent the input image and the guided image, n represents the total number of pixels in the image, and γ is a regularization constant for preventing a i from being too large, and γ = 0.1. In the structure of the anisotropic guided filter, the variable a i is a scale factor for transferring details, and b i is a bias factor for adjusting the intensity of the region;
[0045] 12.1.2) solves a i and b i by the least square method, which are respectively expressed as:
[0046]
[0047]
[0048] wherein, and and respectively represent the mean of the input image and the guided image, and the diffusion effect of the anisotropic guided filter is related to the covariance between the input and guided regions;
[0049] 12.1.3) normalizes a i by weighted average, so as to obtain strong anisotropy, which is expressed as:
[0050]
[0051]
[0052] wherein, w ijThe weights representing the i-centered j-neighborhood aim to maximize the dilation while preserving the strong edges of the guide image, and are the scale factor and bias factor obtained after the weighted average, respectively;
[0053] 12.1.4) Solving the weights w by the objective function ij , which is expressed as:
[0054]
[0055] where, is the mean value of the guide image g j in the neighborhood window , the obtained weights are:
[0056]
[0057] 12.1.5) Since the weights w ij are not spatially adaptive, the sum of the weights of the central pixel will change, therefore, before obtaining the final weighted parameters and the output image I o , w ij is normalized to 1, and the following are obtained, respectively:
[0058]
[0059]
[0060]
[0061] Further, the step 15) specifically comprises the following steps:
[0062] 15.1) Calculating the mutual information measure MI(F, M) between the reference image F and the floating image M, which is expressed as:
[0063]
[0064] where f and m represent the gray values of the pixels in the reference image and the floating image, respectively, P F (f) and P M (m) are the edge probability distributions of the reference image and the floating image to be registered, respectively, and P FM (f, m) represents the joint probability distribution of the reference image and the floating image to be registered;
[0065] 15.2) Calculating the gray square difference SSD(F, M) between the reference image F and the floating image M, which is expressed as:
[0066]
[0067] In the formula, N represents the total number of pixels of the reference image;
[0068] 15.3) Calculate the mutual information-gray square difference fusion similarity measure MS(F, M) between the reference image F and the floating image M, expressed as:
[0069]
[0070] In the formula, w1 and w2 are the weights of the mutual information measure and the gray square difference measure, w1 is 1, and w2 is 8;
[0071] 15.4) Calculate the regularization of the registration , expressed as:
[0072]
[0073] 15.5) Update the energy function E of the registration i , expressed as:
[0074]
[0075] Compared with the prior art, the present application has the following technical effects:
[0076] The medical image Demons registration method based on adaptive anisotropic regularization provided by the present application firstly adopts an adaptive anisotropic model as a registration regularization term, optimizes the weight by the local neighborhood variance to obtain strong anisotropic regularization, and realizes diffusion maximization by using weighted average, so as to truly simulate the complex motion of organs; secondly, the parameter regularization parameter selection is adaptively performed according to the image content, so as to reduce the influence of parameter selection on the registration accuracy; then, the mutual information-gray square difference fusion similarity measure can effectively solve the problem of gray intensity difference between the images to be registered, and improve the accuracy and robustness of the registration result; in addition, the multi-resolution strategy based on wavelet is introduced to perform coarse-to-fine registration, so as to further improve the execution speed of the registration while ensuring the registration accuracy. It can be seen that the present application has high registration accuracy for medical images with complex deformation at the edge of the target of interest. BRIEF DESCRIPTION OF DRAWINGS
[0077] Figure 1 The flow chart of the method of the present application is shown in the figure;
[0078] Figure 2 The comparison figure of the registration result of the registration method based on the optimization framework is shown in the figure;
[0079] Figure 3 The difference visual result comparison figure of the present application and the registration method based on the optimization framework is shown in the figure;
[0080] Figure 4 Comparison chart of registration results of the present application and registration method based on deep learning;
[0081] Figure 5 Comparison chart of differential visual results of the present application and registration method based on deep learning DETAILED DESCRIPTION
[0082] The specific content of the present application is further explained and described in detail below in combination with embodiments.
[0083] As shown in the figure, a medical image Demons registration method based on adaptive anisotropic regularization, characterized in that, comprising the following steps: Figure 1
[0084] 1) input the reference image F to be registered and the floating image M, and initialize the parameters of the adaptive anisotropic model, specifically including the following steps:
[0085] 1.1) input the reference image F to be registered and the floating image M;
[0086] 1.2) set the regularization threshold
[0087] 1.3) set the maximum regularization kernel radius r max = 8, the minimum regularization kernel radius r min = 5;
[0088] 2) use the wavelet-based multi-resolution strategy to decompose the reference image F and the floating image M from coarse to fine, respectively represented as: F i (i = 1, 2,..., N) and M i (i = 1, 2,..., N), wherein i represents the number of layers of the reference image F and the floating image M;
[0089] 3) initialize the spatial deformation field
[0090] 4) set the number of layers i = 1 of the reference image F and the floating image M to be registered, and start registration from the first layer;
[0091] 5) judge whether i ≤ N is true: if yes, execute step 6); otherwise, jump to step 17);
[0092] 6) initialize the deformation field of the current layer according to the deformation field
[0093] 7) set E i as the energy function of the current layer, Ei-1 is the energy function of the last layer, the condition |E i i-1 is not true, go to step 16); otherwise, go to step 8);
[0094] 8) set is the regularization threshold, the condition is not true, go to step 10); otherwise, go to step 9);
[0095] 9) dynamically select the maximum regularization kernel radius r max ;
[0096] 10) dynamically select the minimum regularization kernel radius r min ;
[0097] 11) calculate the deformation update field between the reference image F and the floating image M to be registered is expressed as:
[0098]
[0099] wherein,
[0100] 12) smooth the deformation update field using the adaptive anisotropic kernel G AA is expressed as: Specifically, the following steps are included:
[0101] 12.1) first calculate the output of the adaptive anisotropic regularization model, and the specific process is as follows:
[0102] 12.1.1) describe the anisotropic guided regularization model using a local linear equation, so as to find the optimal a i and b i by minimizing the cost function, which is expressed as:
[0103]
[0104] wherein, i is the pixel index in the image, I i and g i represent the input image and the guide image respectively, n represents the total number of pixels in the image, γ is a regularization constant for preventing a i from being too large, and γ = 0.1, in the structure of the anisotropic guided filter, the variable a i is a scale factor for transferring details, and b i is a bias factor for adjusting the intensity of the region;
[0105] 12.1.2) solve ai and b i , respectively, are expressed as:
[0106]
[0107]
[0108] wherein, and respectively represent the mean of the input image and the guided image, the diffusion effect of the anisotropic guided filtering is related to the covariance between the input and the guided region, when the structure of the input and the guided region is inconsistent, the solved covariance is very small, so that a i value can be obtained, and strong anisotropic diffusion is further realized;
[0109] 12.1.3) the weighted average specification a i , so that strong anisotropy can be obtained, which is expressed as:
[0110]
[0111]
[0112] wherein, w ij represents the weight of the j-neighborhood centered on i, the goal is to realize the maximization of expansion while maintaining the strong edges of the guided image, and are the scale factor and the bias factor obtained after the weighted average, respectively;
[0113] 12.1.4) solving the weight w ij by the objective function, the objective function is expressed as:
[0114]
[0115] wherein, is the mean of the guided image g j in the neighborhood window , the solved weight is:
[0116]
[0117] 12.1.5) since the weight w ij is not spatially adaptive, the sum of the weights of the center pixels will change, therefore, before obtaining the final weighted parameter and the output image I o , w ij is normalized to 1, respectively, to obtain:
[0118]
[0119]
[0120]
[0121] The regularized output image I o is defined as a guided image g i with a local linear model, denoted as:
[0122]
[0123] where and are local linear coefficients, respectively.
[0124] 12.2) The regularized output image I o is considered as a weighted average of the input image I i , so an explicit adaptive anisotropic kernel G AA is adopted to describe the output image I o , denoted as:
[0125]
[0126] where the field window is centered at i with radius r.
[0127] 12.3) The adaptive anisotropic kernel G AA is calculated, denoted as:
[0128]
[0129] where is the total number of pixels in the neighborhood window , μ i and are the mean and variance of the guided image g i in the field window , i is the center of the field window , k and j represent two different pixel indices, and ε represents a smoothing coefficient.
[0130] 13) Update the deformation field:
[0131] 14) Smooth the deformation field:
[0132] 15) Update the energy function E i by fusing the similarity measure and the regularized deformation term using mutual information-gray square difference, which includes the following steps:
[0133] 15.1) Calculate the mutual information measure MI(F,M) between the reference image F and the floating image M, which is represented as:
[0134]
[0135] where f and m represent the gray scale values of the pixels in the reference image and the floating image respectively, P F (f) and P M (m) are the edge probability distributions of the reference image and the floating image to be registered respectively, P FM (f,m) represents the joint probability distribution of the reference image and the floating image to be registered; the greater the mutual information MI(F,M) between the reference image F and the floating image M, the higher the similarity between the two images; if the mutual information MI(F,M) between the two images is 0, it means that the two images are independent of each other and have no similarity;
[0136] 15.2) Calculate the gray scale square difference SSD(F,M) between the reference image F and the floating image M, which is represented as:
[0137]
[0138] where N represents the total number of pixels in the reference image; the smaller the value of SSD, the better the registration result;
[0139] 15.3) Calculate the mutual information-gray scale square difference fusion similarity measure MS(F,M) between the reference image F and the floating image M, which is represented as:
[0140]
[0141] where w1 and w2 are the weights of the mutual information measure and the gray scale square difference measure respectively, w1 takes the value of 1 and w2 takes the value of 8;
[0142] 15.4) Calculate the regularization of the registration which is represented as:
[0143]
[0144] 15.5) Update the energy function E i of the registration, which is represented as:
[0145]
[0146] 16) Output the deformation field of the current layer
[0147] 17) Output the optimal deformation field which is represented as:
[0148] In order to verify the registration performance of the medical image Demons registration method based on adaptive anisotropic regularization (AADemons) proposed in the embodiment, the registration performance is evaluated from three aspects of comparison analysis with the registration method based on the optimization framework, comparison analysis with the registration method based on deep learning and registration efficiency analysis. The verification environment adopted is that the CPU is IntelCore(R) i9-12900H, the main frequency is 4800MHz, the memory is 16GB, the operating system is Windows 11, and Matlab2022a is used for programming.
[0149] In order to evaluate the performance of the AADemons method, the objective quantitative evaluation and the subjective visual evaluation are used to evaluate the pros and cons of the registration results. The objective quantitative evaluation refers to the calculation of the gray absolute difference SAD, the normalized mutual information NMI, the cross correlation CC, the root mean square error RMSE and the precision improvement rate P e Five ways are used to measure the performance of the registration method, wherein:
[0150] (1) The gray absolute difference SAD is calculated directly by using the gray information between the images to be registered, and the calculation formula is represented as:
[0151]
[0152] The greater the SAD is, the lower the registration accuracy is, and vice versa, the higher the registration accuracy is;
[0153] (2) The normalized mutual information NMI is represented as:
[0154]
[0155] In the formula, H(F) and H(M) respectively represent the Shannon entropy of the reference image F and the floating image M, and H(F, M) represents the joint entropy between the reference image F and the floating image M. The statistical correlation between the two images is calculated by the entropy value, and the greater the NMI value is, the better the registration result is;
[0156] (3) The cross correlation CC is represented as:
[0157]
[0158] In the formula, (i,j) represents the pixel point in the image, m*n represents the number of image pixels, And respectively represent the average value of the pixels in the reference image F and the floating image M. The CC measurement value is between 0 and 1, and the closer the CC value is to 1, the higher the registration accuracy is;
[0159] (4) The root mean square error RMSE is represented as:
[0160]
[0161] The registration effect is measured by the error between the images to be registered. The smaller the RMSE value, the better the registration effect.
[0162] (5) In order to measure the difference between the AADemons method and other methods in terms of registration accuracy, the accuracy improvement rate P e is measured, P e is expressed as:
[0163]
[0164] In the formula, RMSE base represents the root mean square error of the comparison method, RMSE improved represents the root mean square error of the AADemons method.
[0165] The subjective visual evaluation is intuitively evaluated by observing the registered images and the difference results before and after registration. When the registered image is closer to the reference image, that is, the difference result is closer to black, it indicates that the registration accuracy is higher, and vice versa.
[0166] In order to verify the registration performance of AADemons, different modal (CT, MRI) medical images of brain and abdomen were randomly selected in the public data set for a large number of registration experiments. The registration results of five groups of images were randomly selected for display, and were sequentially recorded as Image1 to Image5. The AADemons method of the present embodiment was compared with four methods, namely Active Demons (AD), differential homeomorphism Demons (DD), Active Demons method based on balance coefficient (MIAD), and deformable image registration method (GD). All comparison methods maintain the same parameters, and the maximum number of iterations is set to 200. At the same time, DD, MIAD, GD and AADemons all use a three-layer multi-resolution strategy for registration experiments. The parameters of the MIAD method are selected in the optimal interval, that is, a = 1.5, k = 0.5. The parameters of GD are set to the optimal parameters: r = 5, ε = 0.1, K = 3750.
[0167] Figure 2 The registration results of AADemons and the registration method based on the optimization framework are shown. (a) is the reference image, (b) is the floating image, (c) to (g) are the result images obtained by using AD, DD, MIAD, GD, and AADemons for registration, respectively.
[0168] Figure 3The differential visualization results of AADemons and the registration methods based on optimization framework are shown, wherein: (a) is the differential visualization result before registration, (b)-(f) are respectively the differential visualization results of registration using AD, DD, MIAD, GD and AADemons, from Figure 3 The performance of the AADemons method can be more intuitively judged by observing the accuracy of the registration result through differential visualization display of the reference image and the spatially deformed floating image;
[0169] Comparison Figure 2 and Figure 3 It is found that, compared with the other four methods, the AADemons method proposed in the embodiment can obtain good registration results in the region of medical image tissue organs and other detailed information rich regions, and from Figure 3 It can be seen from the differential visualization display that, compared with the other three methods, the black area of the differential image obtained by the AADemons method is more obvious, and it can be seen that the AADemons method obtains a spatial deformation closer to the motion of complex organs in the registration process, so that the spatial correspondence relationship between the floating image after registration and the reference image is more matched.
[0170] In order to ensure the reliability of the registration result, SAD, NMI, CC, RMSE and P e Five evaluation indexes are used to objectively and quantitatively evaluate the registration results of the above five methods, as shown in Table 1:
[0171] Table 1 Analysis of registration results using five evaluation indexes
[0172]
[0173]
[0174] It can be seen from Table 1 that, compared with the other five methods, the AADemons method has the smallest gray absolute difference SAD and root mean square error RMSE of the registration result, the largest normalized mutual information NMI and cross-correlation CC of the images before and after registration, and the registration accuracy of the AADemons method is improved to a certain extent compared with the differential homeomorphism Demons (DD).
[0175] In summary, SAD, NMI, CC, RMSE and P e The conclusion obtained by the objective and quantitative evaluation of the registration result by the five evaluation indexes is consistent with the observation result of the above subjective visualization, which shows that the AADemons has better registration accuracy compared with the other four methods in terms of registration accuracy.
[0176] In the comparative experiment, ActiveDemons (AD) uses isotropic regularization model for registration, which ignores the spatial anisotropy and deformation information of image gray scale change, resulting in a decrease in registration accuracy in areas rich in detail information such as target edges; the DD method uses isotropic Gaussian to process the regularization process of registration, and uses gray square difference SSD to calculate the similarity of registration, which makes the DD method not suitable for processing complex medical images, thereby reducing the registration accuracy; in the registration process of the MIAD method, due to the non-fixed parameter selection of the balance coefficient k and the homogenization coefficient α, different parameters selected for different medical images have a great influence on the registration accuracy, thereby making the MIAD method not have strong robustness; in the registration process of the GD method, the pseudo-contour generated by SLIC is used to guide the regularization process, and due to the fact that the SLIC pseudo-segmentation method is not suitable for complex medical image processing, the GD method cannot obtain accurate registration results in areas where tissue organs slide or image detail information is rich; the AADemons of the embodiment uses adaptive anisotropic regularization for registration processing, can use local neighborhood variance to optimize the weight to obtain strong anisotropic regularization, and uses weighted average to realize maximum diffusion, which can truly simulate the complex motion of organs, at the same time, the AADemons adaptively selects the regularization parameter according to the image content, effectively reduces the influence of the registration parameter on the registration performance, and uses the mutual information-gray square difference fusion improvement method as the similarity measure of registration, which can effectively solve the problem of gray intensity difference between the reference image and the floating image, and improve the accuracy and robustness of the registration result.
[0177] To further verify the performance of AADemons, it is compared with the registration method based on deep learning VoxelMorph, wherein the training and testing processes of VoxelMorph are performed on Linux configured with two NVIDIA Tesla P100 graphics cards, and the parameter batch size is set to 4, Epoch is set to 5*2000, and the initial value of the training learning rate is set to 10 -4 , which is reduced to 5*10 4 after 6*10 -5 , which is reduced to 2.5*10 4 after 8*10 -5 ;
[0178] The LPBA data set of the brain is selected for registration result comparison. Since the data set is based on three-dimensional images, first, VoxelMorph is used for registration test, and second, the three-dimensional images used in the test process are grouped and sliced, wherein each group of sliced images includes a reference image, a floating image, and a registered image. Then, the AADemons method is used to register the reference image and the floating image in sequence according to the image slice order. After registration, the registration accuracy of the VoxelMorph-registered slice and the AADemons-registered result is calculated with respect to the corresponding reference image. For performance evaluation of the method, subjective visual evaluation and objective quantitative evaluation are used for evaluation. At the same time, a group of images randomly selected in the above experimental process are used for experimental visualization display and registration quantitative index analysis. The visual results of the registration experiment comparison are shown in Figure 4 The difference visual results are shown in Figure 5 Wherein:
[0179] Figure 4 The registration results of AADemons and the deep learning-based registration method are shown, wherein (a) represents a reference image, (b) is a floating image, (c) is a VoxelMorph registration result image, and (d) is an AADemons registration result image.
[0180] Figure 5 The difference visual results of AADemons and the deep learning-based registration method are shown, wherein (a) is the difference visualization of the reference image and the floating image, (b) is the VoxelMorph registration difference visualization result, and (c) is the AADemons registration difference visualization result. Figure 5 The registered floating image and the reference image are difference visualized to show the registration results for more intuitive analysis of the registration results.
[0181] It can be seen from the observation of Figure 4 and Figure 5 that the black area of the AADemons registration result is much more than that of VoxelMorph, especially in the tissue organ sliding motion area, which is more obvious. Therefore, compared with VoxelMorph, the registration accuracy of AADemons is higher.
[0182] In order to further analyze the accuracy of the AADemons registration method, the Figure 4The first to the eighth images were selected for evaluation, and were recorded as data1-data8, respectively. Four indexes, SAD, NMI, CC and RMSE, were used to evaluate the registration performance of VoxelMorph and AADemons. The registration quantitative analysis results are shown in Table 2. As can be seen from Table 2, the SAD and RMSE of the AADemons registration result are lower than those of VoxelMorph, and the NMI and CC are higher than those of VoxelMorph. Therefore, the AADemons registration accuracy is higher, and a better registration result can be obtained. This is consistent with the subjective visualization result obtained by observing the registration result of the eight slices in Fig. 6. Figure 4 and Figure 5
[0183] Table 2 Analysis of registration results using four evaluation indexes
[0184]
[0185] In order to ensure the fairness of comparison, four groups of medical images from different modalities and different tissues and organs were selected for registration in the same experimental environment. In order to eliminate system error, all methods were executed 10 times for each group of images to be registered, and then the efficiency of the registration method was analyzed by calculating the average running time. The selected images were brain MRI (size: 442x429), brain CT (size: 617x482), abdominal CT_1 (size: 844x844) and abdominal CT_2 (size: 544x546). The parameter settings of the six registration methods were consistent with the above accuracy evaluation settings. Table 3 shows the average running time of the five methods for registering the four groups of images to be registered.
[0186] Table 3 Average registration time of six methods on four groups of medical images
[0187]
[0188] It can be seen from Table 3 that the registration time of AADemons is lower than that of the other four methods. Among them, compared with other methods, the registration efficiency of brain MRI image is improved by more than 5%, the registration time of brain CT image is shortened by more than 58%, and the registration efficiency of two abdominal CT images is improved by more than 17% and 20% respectively. This is because AD does not use multi-resolution strategy in the registration process, resulting in registration process in the original size of the image, which increases the registration time. In the registration process of GD, SLIC is used to obtain the pseudo-contour of the reference image to guide the regularization process. Since this process significantly increases the computational load, and the number of pseudo-segmented blocks required is large, the registration time is longer. In the registration process of DD and MIAD, gray mean square error and increased balance coefficient are used for registration driving respectively. Compared with DD and MIAD, AADemons uses adaptive anisotropic regularization model for registration, which can obtain a deformation field more consistent with the organ movement, thus better driving the registration process. At the same time, the mutual information-gray square difference fusion similarity measure is used to measure the similarity of registration, which further strengthens the driving force of the registration process, realizes faster convergence speed, and thus the registration efficiency of AADemons is effectively improved compared with AD, DD, MIAD and GD.
Claims
1. A method of medical image Demons registration based on adaptive anisotropic regularization, characterized in that, The method comprises the following steps: 1) inputting a reference image F to be registered and a floating image M, and initializing parameters of an adaptive anisotropic model; 2) using a wavelet-based multi-resolution strategy to decompose the reference image F and the floating image M from coarse to fine, respectively represented as: F i (i = 1, 2,..., N) and M i (i = 1, 2,..., N), wherein i represents the number of layers of the reference image F and the floating image M; 3) initializing a spatial morphing field between the first layer of reference images F to be registered and the floating image M 4) setting a layer number i=1 of the reference image F and the floating image M to be registered currently, and starting registration from the first layer; 5) judging whether i≤N is true: if yes, executing step 6); otherwise, jumping to step 17); 6) the deformation field of the previous layer initialize the deformation field of the current layer 7) Set E i E is the energy function of the current layer i-1 E is the energy function of the previous layer i E i-1 | < 0.001 is true: if yes, go to step 8); otherwise, go to step 16); 8) Set the regularization threshold, the condition is true: if yes, go to step 9); otherwise, go to step 10); 9) dynamically select the largest regularizing kernel radius r max ; 10) dynamically select the smallest regularized kernel radius r min ; 11) computing a deformation update field between the reference image F to be registered and the floating image M is represented as: In the formulae, 12) Utilizing an adaptive anisotropic kernel G AA Smooth deformation update field is represented as: 13) update the deformation field: 14) Smooth deformation field: 15) updating the energy function Ei by using a mutual information-gray square difference fusion similarity measure and a regularization deformation item; 16) output the morphing field of the current layer 17) output optimal deformation field is represented as:
2. The adaptive anisotropic regularization based medical image Demons registration method according to claim 1, characterized in that, The step 1) specifically comprises the following steps: 1.1) inputting a reference image F to be registered and a floating image M; 1.2) Set regularization threshold 1.3) Set maximum regularization kernel radius r max = 8, minimum regularization kernel radius r min = 5.
3. The adaptive anisotropic regularization based medical image Demons registration method according to claim 1, characterized in that, The step 12) specifically comprises the following steps: 12.1) computing the output of an adaptive anisotropic regularization model, regularizing the output image I o defined as a local linear model of the guiding image g i is represented as: wherein and are local linear coefficients; 12.2) regularizing the output image I o is considered as a weighted average of the input image I i , explicit adaptive anisotropic kernels G AA are used to regularize the output image I o is described by: where the field window centered at i with radius r; 12.3) Compute adaptive anisotropic kernel G AA is denoted as: In the formula, For neighborhood windows The total number of pixels in μ i and The guide image g i In the domain window The mean and variance in the range, where i is the neighborhood window. The center of the image is where k and j represent two different pixel indices, and ε represents the smoothing coefficient.
4. The adaptive anisotropic regularization based medical image Demons registration method according to claim 3, characterized in that, The step 12.1) calculates an output of the adaptive anisotropic regularization model, and the process is as follows: 12.1.1) An anisotropic directionally biased regularization model is adopted by using local linear equations to describe the anisotropy, so as to find the optimal a i and b i are expressed as: where i is the pixel index in the image, I i and g i denote the input image and the guide image, respectively, n represents the total number of pixels in the image, and γ is a regularization constant to prevent a i too large a i is a scale factor for transferring details, and b i is a bias factor for adjusting the region intensity. 12.1.2) Solving for a by least squares i and b i are given by: wherein and respectively denote the mean of the input image and the guiding image, the diffusion effect of the anisotropic guided filtering is related to the covariance between the input and guiding regions. 12.1.3) Specification a by weighted average i Thus, a strong anisotropy is obtained, expressed as: where w ij representing the weight of the j-neighborhood centered at i, aims at maximizing the dilation while preserving the strong edges of the guiding image, and are the scale factor and bias factor obtained after the weighted averaging, respectively; 12.1.4) Solving for the weights w by the objective function ij The objective function is expressed as: wherein To guide the image g j In the neighborhood window The resulting weight is solved as: 12.1.5) due to the weights w ij The sum of the weights of the central pixels changes, so the final weighting parameters and the output image I o are obtained by normalizing the weights w ij to 1, respectively:
5. The adaptive anisotropic regularization based medical image Demons registration method according to claim 1, characterized in that, The step 15) specifically comprises the following steps: 15.1) calculating a mutual information measure MI(F, M) between the reference image F and the floating image M, and the calculation is as follows: where f and m represent the gray scale values of the pixels in the reference and floating images, respectively, O F (f) and P M (n) are the marginal probability distributions of the edges of the reference and floating images, respectively, P FM (f, m) represents the joint probability distribution of the reference and floating images to be registered. 15.2) calculating a gray square difference SSD(F, M) between the reference image F and the floating image M, and the calculation is as follows: In the formula, N represents a total pixel number of the reference image; 15.3) calculating a mutual information-gray square difference fusion similarity measure MS(F, M) between the reference image F and the floating image M, and the calculation is as follows: In the formula, w1 and w2 are weights of the mutual information measure and the gray square difference measure respectively, w1 takes a value of 1, and w2 takes a value of 8; 15.4) Calculation of the regularization of the registration is represented as: 15.5) updating the energy function E of the registration i is expressed as: