A Non-Rigid Image Registration Method Based on Sparse Prior

Abstract

The present invention discloses a non-rigid image registration method based on sparse prior, comprising the following steps: constructing a bounded deformation function of the displacement field between a floating image and a fixed image, and using the deformation function of the displacement field as a regularization term of a variational model; establishing a sparse image registration model by using the data term of the L0 norm according to the sparsity characteristics of the grayscale difference between the registered image and the fixed image; and further adopting the L p The sparse approximation of the pseudo-norm is approximately L p A sparse image registration model is proposed; a primal-dual iterative weighted algorithm is used to solve the registration model and obtain a registration result. The invention is more robust to noise, thereby obtaining a more effective registration result.

Description

Technical Field

[0001] The present invention relates to the technical field of image processing and computer vision, and in particular to a non-rigid registration method for medical images. Background Art

[0002] Image registration involves spatially aligning images acquired from different sensors, at different times, or under different conditions for subsequent analysis and processing. It has been widely used in many fields such as medical image analysis and computer vision. The essence of image registration is to find a spatial geometric transformation between two or more images. Since the imaging objects of medical images are mostly human soft tissues or soft organs, they undergo non-rigid deformation when affected by breathing, external pressure, or lesions such as tumors. Therefore, medical image registration usually adopts non-rigid registration methods [1]. From a mathematical point of view, medical images can be regarded as a collection of measurement values of a physical quantity (such as ultrasound echo signals or X-ray radiation absorption) in a specific two-dimensional or three-dimensional space. These measurement values are usually called image gray values or intensity values. These specific measurement positions are called pixels and voxels in two-dimensional and three-dimensional situations, respectively. Generally, this measurement space is denoted as d = 2, 3 represents the dimension of the space. The image I can be viewed as a function from the image domain Ω to the range R.

[0003] Image registration generally involves two or more images. This paper focuses on the non-rigid registration between two medical images. For two medical images to be registered, one is generally called a fixed image F and the other is called a floating image M, denoted as F and M:Ω→R respectively. They can be regarded as two scalar functions defined on the image domain Ω. The goal of registration is to find the displacement field u:Ω→R between the corresponding pixel coordinates of the two images to be registered. d ,u=(u1,…,u d ),u1,…,u d is the component of the displacement field u[1]. Adding this displacement field to the floating image M will yield the registered image M. * =M(x+u(x)), so that the registered image M * The pixel points in the fixed image F correspond to the same point of the object, that is,

[0004] Here, the accuracy of the correspondence between M* and F is often expressed by a similarity metric calculated based on the pixel values (voxel values) of the image. After the similarity metric is given, the non-rigid registration problem can be regarded as an optimization problem. Many existing non-rigid image registration methods transform the non-rigid registration problem into a corresponding variational model under the framework of functional analysis or directly obtain the corresponding partial differential equation system from the physical mechanics model. This type of method usually regards the displacement field u between the images to be registered as a function in a specific function space. In the framework of functional analysis, the non-rigid registration problem can generally be written as [2]:

[0005] min u E(u)=λ·Data(M,F,u)+Regu(u),

[0006] That is, in the form of data term + regularization term. Here, the data term is a similarity measure, and λ is a weight coefficient that balances the data term and the regularization term.

[0007] For data items, the more applicable and commonly used similarity metrics [3] include normalized mutual information (NMI), the L1 norm of image difference (SAD), the L2 norm of image difference (SSD), and the local correlation coefficient (LCC). However, these similarity metrics are designed based on the statistics of pixel values, are sensitive to noise in the image, and cannot characterize the sparsity of the image difference after registration. The L0 norm is a classic metric function used to characterize data sparsity and is widely used in computer vision and machine learning [4].

[0008] Regarding the regularization term, existing variational models can be roughly divided into two categories. The first category is to find smooth displacement fields [2,3], such as elastic models, fluid models, optical flow models, diffusion models, and differential homeomorphism models. The second category is to find non-smooth or even discontinuous displacement fields, such as registration methods based on total variation [5]. Bounded variation functions (BV functions) have been widely used in many fields of image processing and analysis and have been proven to be a very effective type of function for image modeling [6]. However, the components of the displacement field are treated as independent functions and are processed and regularized separately. This is not reasonable because the components of the displacement field are actually coupled. In our previous work, medical images were regarded as elastic-plastic materials [7], allowing the displacement field to be locally anisotropic and discontinuous, and bounded deformation functions (BD functions) were used to characterize the displacement field that may have discontinuous points. We used the SSD data term for the same-modality medical image registration problem and achieved good results in the non-rigid registration task of medical images.

[0009] The cited references are as follows:

[0010] [1] J. Modersitzki. Image registration for radiology: A review. International Journal of Computer Assisted Radiology and Surgery[J], 4(2): 111 - 121, 2009.

[0011] [2] Sotiras A, Davatzikos C, Paragios N. Deformable medical image registration: A survey. [J]. IEEE Transactions on Medical Imaging, 32(7): 1153 - 1190, 2013.

[0012] [3] Modersitzki J. FAIR: Flexible Algorithms for Image Registration[M]. Vol.

[0013] 6. Philadelphia: SIAM, 2009.

[0014] [4] Krig S. Computer Vision Metrics: Textbook Edition[M]. Springer Cham, 2019.

[0015] [5] Vishnevskiy V, Gass T, Székely G, et al. Total Variation Regularization of Displacements in Parametric Image Registration[J]. IEEE Transactions on Medical Imaging, 36(2): 385 - 395, 2017.

[0016] [6] Pragliola M. et al. On and Beyond Total Variation Regularization in Imaging: The Role of Space Variance[J]. SIAM Review 65(3): 601 - 685, 2023.

[0017] [7]Nie,Ziwei and Yang,Xiaoping.Deformable Image Registration Using Functions of Bounded Deformation[J].IEEE Transactions on Medical Imaging,38(6)2019:1488-1500. Summary of the Invention

[0018] Objective of the Invention: Aiming at the defects and deficiencies of the prior art, the present invention proposes a non-rigid image registration method based on sparse prior, which is robust to image noise and allows the displacement field between images to be discontinuous locally, thus achieving better registration results.

[0019] Technical Solution: Aiming at the sparsity of the gray difference between the registered image and the fixed image and the robustness to noise, the present invention uses the L0 norm of the difference between the registered images as the data term in the non-rigid registration problem, and at the same time regularizes the displacement field into a bounded deformation function, and gives the corresponding iterative weighted primal-dual algorithm to solve the established sparse non-rigid registration model. Specifically, a non-rigid image registration method based on sparse prior includes the following steps:

[0020] (l) Construct a bounded deformation function of the displacement field u between the images to be registered;

[0021] (2) According to the sparsity of the gray difference between the registered image and the fixed image, use the L0 norm as the data term:

[0022] (3) Take the bounded deformation function of the displacement field u constructed in step (1) as the regularization term of the variational model, and combine the L0 sparse data term defined in step (2) to establish a sparse image registration model;

[0023] (4) Sparsely approximate the sparse image registration model established in step (3) using the L p (0 < p < 1) quasi-norm to obtain the L p sparse image registration model;

[0024] (5) Linearize the data term in the L p sparse image registration model in step (4), and use the primal-dual iterative weighted algorithm to solve it to obtain the registration result.

[0025] Beneficial effects: The present invention utilizes bounded deformation regularization terms to characterize non-smooth or even discontinuous displacement fields, and utilizes the advantage of L0 norm's robustness to noise to establish data items of the registration model. The established model has lower requirements for the regularity of the displacement field, and the sparse prior of the data items is more in line with the actual situation of the difference between medical images after registration, that is, the grayscale difference between the registered image and the fixed image has a lot of sparsity, thereby being able to obtain better registration effects while satisfying the bounded deformation.

[0026] Contents of attached figure

[0027] Figure 1 This is an overall flow chart of the image registration method of the present invention;

[0028] Figure 2 A schematic diagram comparing the registration results of a pair of noisy artificial images between the present invention and the prior art;

[0029] Figure 3 Schematic diagram of the displacement field of the present invention and the prior art in registering a pair of artificial images containing noise;

[0030] Figure 4 Schematic diagram comparing the prostate MR image registration results of the present invention and the prior art.

[0031] Figure 5 Schematic diagram of the displacement field of prostate MR image registration between the present invention and the prior art. DETAILED DESCRIPTION

[0032] In order to make the objectives, technical solutions and advantages of the present invention more clear, the embodiments of the present invention will be further described in detail below with reference to the accompanying drawings.

[0033] like Figure 1 As shown, the invention provides a non-rigid image registration method for processing discontinuous deformation, comprising the following steps:

[0034] (1) Consider the displacement field u between the floating image M and the fixed image F as a bounded deformation function, and define the deformation TD(u) of the displacement field u:

[0035]

[0036] in express is a first-order differentiable symmetric tensor with compact support, express The divergence of Ω represents the d-dimensional image area, d = 2 or 3, respectively, represents the 2-dimensional or 3-dimensional image area, Sym 2 (R d ) represents the second-order d-dimensional symmetric tensor space, that is, all d×d symmetric matrices.

[0037] Correspondingly, the bounded deformation function space BD(Ω) is defined as:

[0038] {u ∈ L 1 (Ω) | TD(u) < ∞}

[0039] where L 1 (Ω) represents the integrable space.

[0040] (2) According to the gray - level distribution characteristics between the images to be registered, that is, the sparsity of the gray - level difference between the registered image M * and the fixed image F, the L0 data term is adopted. Specifically, it means that the floating image M in the images to be registered is subtracted from the fixed image F in the images to be registered after the action of the displacement field u, and the L0 norm is used as the data term of the variational model. The specific form is:

[0041] Data(M, F, u) = ‖M(x + u(x)) - F(x)‖0.

[0042] where x is the pixel point in the image region.

[0043] (3) Take TD(u) defined in step (1) as the regularization term of the variational model, and combine with the L0 sparse data term defined in step (2) to establish a sparse image registration model:

[0044] u * = argmin u∈BD(Ω) {λ · Data(M, F, u)+TD(u)},

[0045] The feasible solution space BD(Ω) in the above registration model is the bounded deformation function space.

[0046] (4) Due to the introduction of the L0 norm, the sparse image registration model is an NP - hard optimization problem. Therefore, the data term in the sparse image registration model in step (3) is approximated by the L p (0 < p < 1) quasi - norm, and the L p sparse image registration model is obtained as follows:

[0047]

[0048] (5) M(x + u(x)) in the L p sparse image registration model has high non - convexity, making the optimization problem difficult to solve. Considering the case where the displacement field is small enough and the floating image M is smooth, through the first - order Taylor expansion, the L p sparse image registration model is linearized as:

[0049]

[0050] (6) After linearization, L p The sparse image registration model is discretized and an iterative weighted algorithm is used to obtain u k The update rule is:

[0051]

[0052] in

[0053] (7) Using Fenchel duality, according to where X and Y are dual spaces, ‖·‖1 and ‖·‖ ∞ represents the 1-norm and ∞-norm in the tensor sense, L p The sparse image registration model can be transformed into the following saddle point problem:

[0054]

[0055] in

[0056] (8) The saddle point problem in step (4) is solved using the adaptive primal-dual algorithm to obtain the registration result.

[0057] In the present invention, λ=15, p=0.5 is used in the numerical experiments. The following primal-dual iterative weighted algorithm is given:

[0058] Initialization: Input the images to be registered M and F, let u 0,0 =0,q 0,0 =0, U * =u 0,0 ;

[0059] Select: step size parameter τ, σ>0, θ∈[0,1];

[0060] Iterate double loop: perform outer loop and inner loop calculations as follows until convergence or the maximum number of iterations is reached, and return the displacement field U obtained by alignment * and the registered image M * ;

[0061] (8.1) k-fold outer loop: Use the update rule of the iterative weighted algorithm to obtain the weights Updates;

[0062] (8.2) n-fold inner loop: Update the dual variable q using the update rule of the primal-dual algorithm k,n+1 , the original variable u k,n+1 ,

[0063] (8.3) Let n = n + 1 and enter the next round of iteration until the condition for stopping the inner loop is met;

[0064] (8.4) Interpolate the image, that is, M k+1 =M k (x+u k,n+1 );

[0065] (8.5) Let k = k + 1 and enter the next round of iteration until the condition for stopping the outer loop is met.

[0066] The specific steps are shown in the following table:

[0067]

[0068]

[0069] In the above algorithm, E represents the energy functional, m inner Indicates the maximum number of iterations of the inner loop, for fixed k, n k Indicates the number of iterations that meet the inner stopping condition, m outer In the numerical experiments of all the examples in this invention, τ=σ=0.1,θ=1,m inner =40,m outer =400.

[0070] The above describes the detailed implementation steps of the non-rigid image registration method proposed in the present invention. In order to more clearly understand the advantages of the present invention, the following two examples are used to compare the registration results of the present invention with those of other methods.

[0071] Example 1:

[0072] This experiment is to register a pair of artificial images with noise, such as Figure 2 As shown, (a) is a fixed image, (b) is a floating image, (c) is the L2TV model (SSD data item, TV regularization item) registration result, (d) is the L2BD model (SSD data item, BD regularization item) registration result, (e) is the L2Diff model (SSD data item, differential homeomorphism regularization item) registration result, (f) is the LpBD model registration result of the sparse prior registration proposed in this invention. Figure 3 As shown, (a) is the displacement field of the L2TV model registration, (b) is the displacement field of the L2BD model registration, (c) is the displacement field of the L2Diff model registration, and (d) is the displacement field of the LpBD model of the present invention.

[0073] Furthermore, in order to compare different non-rigid registration models, we use some commonly used evaluation indicators to evaluate the registration results of each model. Here we use the classic mean structure similarity index measure (mSSIM), mean normalized correlation coefficient (NCC) and mean square error (MSE) to evaluate the similarity between the image registration result and the fixed image. Among them, the values of the two evaluation indicators mSSIM and mNCC are both in [0,1]. The larger the value, the higher the similarity, that is, the better the registration result. The value range of MSE is [0,1]. The smaller the value, the smaller the registration error, that is, the better the registration result.

[0074] Table 1 compares the registration results obtained using the L2TV model, L2BD model, L2Diff model, and the LpBD model of the present invention for Example 1, using three quantitative metrics: mean structural similarity (mSSIM), mean normalized correlation coefficient (mNCC), and mean square error (MSE). The optimal value is highlighted in bold. This demonstrates the effectiveness of the LpBD model proposed in this paper.

[0075] Table 1 Quantitative indicators of manual image registration results

[0076] Not registered L2TV L2BD L2Diff LpD mSSIM 0.158437 0.77603 0.84661 0.884362 0.989993 mNCC 0.725041 0.994241 0.989857 0.995458 0.998635 MSE 0.083413 0.001839 0.00331 0.001433 0.000429

[0077] Example 2:

[0078] This example applies the sparse prior-based non-rigid medical image registration LpBD model to the registration of prostate MR images. Figure 4 As shown, (a) is a fixed image, (b) is a floating image, (c) is the L2TV model registration result, (d) is the L2BD model registration result, (e) is the L2Diff registration result, and (f) is the LpBD model registration result proposed in this invention. Figure 5 As shown, (a) is the displacement field of the L2TV model registration, (b) is the displacement field of the L2BD model registration, (c) is the displacement field of the L2Diff model registration, and (d) is the displacement field of the LpBD model proposed in this invention.

[0079] Table 2 compares the registration results for Example 2 using three quantitative metrics: structural similarity (mSSIM), mean normalized correlation coefficient (mNCC), and mean square error (MSE). The L2TV model, L2BD model, L2Diff model, and the proposed LpBD model are shown in bold. The optimal values are shown in bold. Table 2 demonstrates the effectiveness of the proposed LpBD model.

[0080] Table 2 Quantitative indicators of prostate MR image registration results

[0081] Not registered L2TV L2BD L2Diff LpD SSIM 0.771722 0.830713 0.833644 0.774724 0.949622 NCC 0.898121 0.958826 0.959946 0.901464 0.991738 MSE 0.010484 0.003709 0.003604 0.010111 0.00073

[0082] In summary, the LpBD model based on sparse prior proposed in the present invention can obtain more effective registration results.

Claims

1. A non-rigid image registration method based on sparse prior, characterized in that: The following steps are involved: (1) Constructing a bounded deformation function of the displacement field u between the images to be registered; (2) According to the characteristics of the consistent grayscale distribution of the image to be registered and the sparsity of the grayscale difference between the registered image and the fixed image, the L0 norm is used as the data item of the variational model; (3) Using the bounded deformation function of the displacement field u constructed in step (1) as the regularization term of the variational model, combined with the data term of the L0 norm defined in step (2), a sparse image registration model is established; (4) According to the sparse approximation of the L p quasi-norm, where 0 < p < 1, approximate the sparse image registration model established in step (3) as the L p sparse image registration model; (5) For L in step (4) p The sparse image registration model is locally linearized, and the primal-dual iterative weighted algorithm is used to solve L in step (4). p Sparse image registration model to obtain the registration result; The bounded deformation function form of the displacement field u in step (1) is as follows: in express is a first-order differentiable symmetric tensor with compact support, express The divergence of Ω represents the d-dimensional image area, d = 2 or 3, respectively, represents the 2-dimensional or 3-dimensional image area, Sym 2 (R d ) represents the second-order d-dimensional symmetric tensor space, that is, all d×d symmetric matrices.

2. The non-rigid image registration method based on sparse prior according to claim 1, characterized in that: The specific form of the variational model data item in step (2) is: Data(M,F,u)=‖M(x+u(x))-F(x)‖0, Where M is the floating image in the image to be registered, F is the fixed image in the image to be registered, and x is the pixel point of the image area.

3. The non-rigid image registration method based on sparse prior according to claim 2, characterized in that: The sparse image registration model in step (3) is: he * =argmin u∈BD(Ω) {E(u)=λ·Data(M,F,u)+TD(u)}, The feasible solution space BD(Ω) in the above registration model is {u∈L 1 (Ω)|TD(u)<∞} is a bounded deformation function space, L 1 (Ω) represents the space of integrable functions.

4. The non-rigid image registration method based on sparse prior according to claim 3, characterized in that: The step (4) is based on the L0 norm using L p Quasi-norm approximation, the sparse image registration model is approximated as L p The sparse image registration model is as follows: u * =argmin u∈BD(Ω) {E(u)=λ∫ Ω |M(x+u(x))-F(x)| p dx+TD(u)}。 5. The non-rigid image registration method based on sparse prior according to claim 4, characterized in that: The step (5) comprises: (5.1) for L p The sparse image registration model is linearized as: (5.2) For the linearized L p The sparse image registration model is discretized and an iterative weighted algorithm is performed to obtain u k The update rule is: in (5.3) According to X and Y are dual spaces, ‖·‖1 and ‖·‖ ∞ represents the 1-norm and ∞-norm in the tensor sense, L p The sparse image registration model is transformed into the following saddle point problem: in 6. The non-rigid image registration method based on sparse prior according to claim 5, characterized in that: The primal-dual iterative weighted algorithm comprises the following steps: Initialization: Input the images to be registered M and F, let u 0,0 =0,q 0,0 =0, U * =u 0,0 ; Select: step size parameter τ, σ>0, θ∈[0,1]; Iterate double loop: perform outer loop and inner loop calculations as follows until convergence or the maximum number of iterations is reached, and return the displacement field U obtained by alignment * and the registered image M * ; (7.1) k-fold outer loop: Use the update rule of the iterative weighted algorithm to obtain the weights Updates; (7.2) n-fold inner loop: Update the dual variable q using the update rule of the primal-dual algorithm k,n+1 , the original variable u k,n+1 , (7.3) Let n = n + 1 and enter the next round of iteration until the condition for stopping the inner loop is met; (7.4) Interpolate the image, that is, M k+1 =M k (x+u k,n+1 ); (7.5) Let k = k + 1 and enter the next round of iteration until the condition for stopping the outer loop is met.

Images