Non-rigid image registration method based on sparse prior

Abstract

The invention discloses a non-rigid image registration method based on sparse prior, and the method comprises the following steps: constructing a bounded deformation function of a displacement field between a floating image and a fixed image, and taking the deformation function of the displacement field as a regular term of a variational model; establishing a sparse image registration model by adopting a data item of an L0 norm according to the sparsity characteristic of the gray difference between the registered image and the fixed image; further carrying out sparse approximation by adopting an Lp quasi-norm to obtain an Lp sparse image registration model through approximation; and solving the registration model by using a primal dual iterative weighting algorithm to obtain a registration result. According to the method, the robustness to noise is higher, so that a more effective registration result is obtained.

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 regarded as a function from the image domain Ω to the range R.

[0003] Image registration generally involves two or more images. The present invention 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, which are denoted as F, 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 gives the registered image M * =M(x+u(x)), so that the registered image M * The pixels in the image and the pixels 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 measure calculated based on the pixel values ​​(voxel values) of the image. After the similarity measure 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, which is used to balance 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), image difference L1 norm (SAD), image difference L2 norm (SSD) and 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] For the regularization term, the existing variational models can be roughly divided into two categories. The first category is to find a smooth displacement field [2,3], such as elastic models, fluid models, optical flow models, diffusion models and differential homeomorphism models. The second category is to find a non-smooth or even discontinuous displacement field, such as the registration method 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 class of functions 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. SSD data terms were used for the same-modality medical image registration problem, and good results were achieved 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 have discontinuities locally, thus achieving better registration results.

[0019] Technical Solution: Considering the sparsity of the gray - level difference between the registered image and the fixed image and its 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. At the same time, the displacement field is regularized into a function of bounded deformation, and a corresponding iterative weighted primal - dual algorithm is given 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] (1) Construct a function of bounded deformation for the displacement field u between the images to be registered;

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

[0022] (3) Take the function of bounded deformation of the displacement field u constructed 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;

[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 solve it using the primal - dual iterative weighted algorithm to obtain the registration result.

[0025] Beneficial effect: 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 on 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 in 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 premise of bounded deformation.

[0026] Figure content

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

[0028] Figure 2 A schematic diagram showing the comparison of the registration results of a pair of artificial images containing noise between the present invention and the prior art;

[0029] Figure 3 A 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 It is a schematic diagram comparing the prostate MR image registration results of the present invention and the prior art.

[0031] Figure 5 The figure is a schematic diagram of the displacement field of the 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, a non-rigid image registration method for processing discontinuous deformations of the invention comprises the following steps:

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

[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] Accordingly, 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-scale distribution characteristics between the images to be registered, that is, the sparsity of the gray-scale 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 image to be registered is added with the action of the displacement field u and then subtracted from the fixed image F in the image to be registered, 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) Taking TD(u) defined in step (1) as the regularization term of the variational model and combining the L0 sparse data term defined in step (2), a sparse image registration model is established:

[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 an adaptive primal-dual algorithm to obtain the registration result.

[0057] In the present invention, λ=15, p=0.5 are used in 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] Iteration 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 the present 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 registration results of the present invention and other methods are compared and explained through two examples below.

[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 the present invention. Figure 3 As shown, (a) is the displacement field of L2TV model registration, (b) is the displacement field of L2BD model registration, (c) is the displacement field of L2Diff model registration, and (d) is the displacement field of 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 results 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 by using the L2TV model, L2BD model, L2Diff model and the LpBD model of the present invention for Example 1 using three quantitative indicators: mean structural similarity mSSIM, mean normalized correlation coefficient mNCC and mean square error MSE. The optimal value is bolded. It can be seen from this that the LpBD model proposed in the present invention is effective.

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

[0076] Not registered L2TV L2BD L2Diff LqCy 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 by the present invention. Figure 5 As shown, (a) is the displacement field of L2TV model registration, (b) is the displacement field of L2BD model registration, (c) is the displacement field of L2Diff model registration, and (d) is the displacement field of LpBD model proposed in the present invention.

[0079] Table 2 compares the registration results of Example 2 under the L2TV model, L2BD model, L2Diff model and the LpBD model proposed by the present invention using three quantitative indicators: structural similarity mSSIM, mean normalized correlation coefficient mNCC and mean square error MSE. The optimal value is bolded. Table 2 illustrates the effectiveness of the LpBD model proposed by the present invention.

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

[0081] Not registered L2TV L2BD L2Diff LqCy 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 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), to establish a sparse image registration model; (4) According to the sparse approximation of the L p quasinorm, 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, and obtain the registration result.

2. The non-rigid image registration method based on sparse prior according to claim 1, characterized in that: 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.

3. The non-rigid image registration method based on sparse prior according to claim 2, 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.

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

5. The non-rigid image registration method based on sparse prior according to claim 4, characterized in that: The step (4) is based on the L0 norm using L p Quasi-norm approximation, approximating the sparse image registration model 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)}。 6. The non-rigid image registration method based on sparse prior according to claim 5, 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 7. The non-rigid image registration method based on sparse prior according to claim 6, 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]; Iteration 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 the inner loop to stop 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