A CS-MRI Image Reconstruction Method Based on Fractional Norm Low-Rank Constraint

Abstract

The present invention discloses a CS-MRI image reconstruction method based on fractional norm low-rank constraint. It belongs to the technical field of digital image processing. It is an image reconstruction method that uses the correlation between similar image patches to perform low-rank constraint on the matrix composed of them. First, similar image patches are found according to the Euclidean distance and a similar image patch matrix is constructed. Then, a fractional norm is used to constrain the low-rank property of the similar image patch matrix. Finally, the alternating direction multiplier method is used to optimize and solve the established image reconstruction model. The present invention applies a low-rank constraint of fractional norm to the similar image patch matrix, effectively extracts the structural similarity of the image, and at the same time accurately solves each sub-problem of the image reconstruction model, so that the obtained image effectively suppresses the artifacts caused by undersampling, restores a large number of detailed textures, and is closer to the real image in terms of visual effect. Therefore, it can be used for the reconstruction of medical images.

Description

Technical Field

[0001] The present invention belongs to the technical field of digital image processing, and particularly relates to performing fractional norm low-rank constraint on a similar image block matrix to realize CS-MRI image reconstruction for high-quality medical image restoration. Background Art

[0002] Magnetic resonance imaging (MRI) has the advantages of no ionizing radiation and high resolution, and is widely used in clinical diagnosis and medical research. However, its long scanning time is likely to cause motion artifacts in the image or discomfort to the patient. The compressive sensing (CS) theory proposed in recent years points out that when a signal is sparse, the signal can be accurately recovered from a small amount of random linear measurement data. Therefore, magnetic resonance imaging based on the compressive sensing theory allows for significantly shortening the scanning time by randomly undersampling K-space data. On this basis, establishing a sparse representation and reconstruction model adapted to the image and designing a corresponding optimization algorithm are the keys to obtaining high-quality reconstructed images.

[0003] Since there are a large number of similar structures in the image, compared with sparse representation of a single image block, constraining the low-rank property of a similar image block matrix as an object can make more full use of the characteristics of the image itself to improve the reconstruction quality of the image. However, existing CS-MRI reconstruction methods based on low-rank constraint of a similar image block matrix generally use unbounded functions such as the nuclear norm or the logarithm of the determinant of the matrix to approximately replace the rank of the matrix. Therefore, the low-rank constraint on the similar image block matrix is restricted to a certain extent, and at the same time, it is difficult to avoid local optimal solutions in the optimization and solution process of the low-rank matrix, which seriously affects the final image reconstruction performance.

[0004] In fact, the fractional norm of the singular value of a matrix has continuous and bounded properties, and its approximation degree to the rank of the matrix is better than that of existing low-rank regular functions. Therefore, the fractional norm can be used to effectively constrain the low-rank property of the similar image block matrix, and then establish an image reconstruction model and perform accurate solution to enhance the restoration accuracy of the similar image block and make the obtained reconstructed image closer to the real image. Summary of the Invention

[0005] The purpose of the present invention is to make full use of the correlation between similar image blocks in the image and propose a CS-MRI image reconstruction method based on fractional norm low-rank constraint. This method uses a fractional norm to effectively constrain the low-rank property of the similar image block matrix, and establishes a corresponding image reconstruction model. At the same time, in the process of solving the reconstruction model using the alternating direction method of multipliers, by accurately solving the sub-problems of each variable, the obtained reconstructed image can effectively suppress the artifact phenomenon and restore a large amount of detailed information. Specifically, it includes the following steps:

[0006] (1) Input the K-space data y of an MRI and the sampling template, and perform pre-reconstruction using traditional reconstruction methods to obtain the initial reconstructed image x (0) ;

[0007] (2) For each -sized target image patch, respectively delimit a fixed-sized rectangular region centered on it in the initial reconstructed image x (0) , and find m image patches (including the target image patch itself) with the smallest Euclidean distance from it within this region to form a similar image patch matrix, represents the i-th similar image patch matrix, where x is the image to be reconstructed, and R i,j is the extraction matrix of the j-th similar image patch in X i , and is the extraction operator of the entire similar image patch matrix;

[0008] (3) Use the L a,ε fractional norm to constrain its low-rank property, represents the L i fractional norm of X a,ε , and its definition formula is:

[0009]

[0010] where, σ k is the k-th singular value of X i , a > 0 is the exponential factor and its value is not greater than 1, and ε > 0 is a small constant;

[0011] (4) On the basis of using the L a,ε fractional norm to constrain the low-rank property of the similar image patch matrix, establish a reconstruction model for the image and the similar image patch matrix:

[0012]

[0013] where, represents the square of the two-norm of the vector, and F u = UF is the undersampled Fourier encoding matrix, U is the undersampled matrix, F is the orthogonal Fourier transform matrix, and λ is the regularization parameter of the K-space data fidelity term;

[0014] (5) Introduce the Lagrange multiplier B i and the penalty parameter β, and solve each optimization variable in the reconstruction model by the alternating direction method of multipliers:

[0015] (5a) Fix x, B i and β, and solve the sub-problem for the similar image patch matrix X i :

[0016]

[0017] Among them, represents the square of the Frobenius norm of the matrix;

[0018] (5b) Fix X i , B i and β, and solve the sub-problem of the image x:

[0019]

[0020] (6) After obtaining the optimal solutions of the variables in the reconstruction model, the update formulas for the Lagrange multiplier B i and the penalty parameter β are:

[0021]

[0022] β = γβ

[0023] where γ > 1 is the growth factor. Repeat steps (5) - (6) until the reconstructed image meets the conditions or the number of iterations reaches the preset upper limit.

[0024] The innovation of the present invention is to take the matrix composed of similar image patches found in the image as a whole processing object, and realize the low-rank property constraint of the similar image patch matrix through the L a,ε fractional norm; use the alternating direction multiplier method to solve the corresponding image reconstruction model, and calculate the optimal solution for the sub-problem of the similar image patch matrix through an iterative method under this optimization framework, so that the low-rank property and restoration accuracy of the obtained similar image patch matrix are significantly improved. At the same time, the sub-problem of the image is solved by fast approximation to obtain a high-quality reconstruction result.

[0025] The beneficial effects of the present invention: Taking the similar image patch matrix as the processing object, making full use of the correlation between similar image patches; using the L a,ε fractional norm to constrain the low-rank property of the similar image patch matrix, enhancing the strength of the low-rank constraint compared with the existing low-rank regular functions; under the optimization framework of the alternating direction multiplier method, accurately solving the sub-problems of each variable in the image reconstruction model. Therefore, the reconstructed image restores a large amount of details and effectively suppresses artifacts, and has a good visual effect.

[0026] The present invention is mainly verified by simulation experiments, and all steps and conclusions are verified correctly on MATLAB8.0. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] Figure 1 is the flowchart of the working process of the present invention;

[0028] Figure 2is the water balloon MRI image used in the simulation of the present invention;

[0029] Figure 3 are the reconstruction results of the water balloon MRI image with a sampling rate of 10% using different methods;

[0030] Figure 4 is the error corresponding to the reconstruction results of the water balloon MRI image with a sampling rate of 10% using different methods. Detailed implementation manners

[0031] Referring to Figure 1 , the present invention is a CS-MRI image reconstruction method based on fractional norm low-rank constraint, and the specific steps are as follows:

[0032] Step 1, perform pre-reconstruction on the image and construct a similar image block matrix for the target image block.

[0033] (1a) Input a piece of K-space data y of MRI and a sampling template, and perform pre-reconstruction using the traditional image reconstruction method with sparse constraint in the wavelet domain to obtain an initial reconstructed image x (0) ;

[0034] (1b) Extract target image blocks of size (0) at equal distances along the vertical and horizontal directions in the initial reconstructed image x ;

[0035] (1c) For each target image block, respectively delimit a fixed-size rectangular area centered on it in the initial reconstructed image x (0) , and find m image blocks (including the target image block itself) with the smallest Euclidean distance from it in this area to form a similar image block matrix, denotes the i-th similar image block matrix, where x is the image to be reconstructed, R i,j is the extraction matrix of the j-th similar image block in X i , and is the extraction operator of the entire similar image block matrix.

[0036] Step 2, according to the correlation between similar image blocks, adopt the low-rank characteristic of L a,ε fractional norm constraint for the constructed similar image block matrix, so as to establish an image reconstruction model.

[0037] (2a) Define the L a,ε fractional norm of the similar image block matrix:

[0038]

[0039] where σ k is for X iThe k-th singular value, where a > 0 is an exponential factor with a value not greater than 1, and ε > 0 is a small constant;

[0040] (2b) Based on the low-rank property of the fractional norm-constrained similar image patch matrix, a reconstruction model for the image and the similar image patch matrix is established: a,ε The fractional norm constraints the low-rank property of the similar image patch matrix, and a reconstruction model regarding the image and the similar image patch matrix is established:

[0041]

[0042] Wherein, denotes the square of the two-norm of the vector, and F u = UF is the undersampled Fourier encoding matrix, U is the undersampling matrix, F is the orthogonal Fourier transform matrix, and λ is the regularization parameter of the K-space data fidelity term.

[0043] Step 3, Solve the established reconstruction model by the alternating direction method of multipliers.

[0044] (3a) Since the established reconstruction model is a multivariate minimization problem with equality constraints, the corresponding augmented Lagrangian function has the following expression:

[0045]

[0046] Wherein, B i denotes the Lagrange multiplier, β denotes the penalty parameter, Re(·) denotes the real part function of a complex number, tr(·) denotes the trace function of a matrix, is the conjugate transpose matrix of the Lagrange multiplier B i , denotes the square of the Frobenius norm of the matrix. Further, according to and the alternating direction method of multipliers, each optimization variable in the reconstruction model is solved;

[0047] (3b) Fix x, B i and β, and solve the sub-problem regarding the similar image patch matrix X i :

[0048]

[0049] This sub-problem can be decomposed into a scalar optimization problem regarding the singular values of X i and is solved according to the following steps:

[0050] (3b1) Perform singular value decomposition on the matrix to obtain:

[0051]

[0052] Wherein, G i is the left singular vector matrix, and Vi is the right singular vector matrix, δ1, δ2, …, δ m is the matrix 's singular values, diag(δ1, δ2, …, δ m ) represents the diagonal matrix with δ1, δ2, …, δ m as diagonal elements;

[0053] (3b2) decomposes this sub-problem into scalar optimization problems regarding each singular value of X i , where the scalar optimization problem regarding the k-th singular value σ i of X k is:

[0054]

[0055] (3b3) Solves the optimal solution of the scalar optimization problem regarding σ k according to the fixed-point iteration method. The specific steps are to set the initial value of σ k to and calculate the iteration formula:

[0056]

[0057] v = v + 1

[0058] where v represents the number of iterations, represents the calculated value of σ k in the v-th iteration. Repeat calculating this iteration formula until v reaches the preset upper limit;

[0059] (3b4) After obtaining the calculated value of the final iteration, the discriminant condition for determining whether is the optimal solution is:

[0060]

[0061] When satisfies the discriminant condition, let When does not satisfy the discriminant condition, let σ k = 0;

[0062] (3b5) After solving all the singular values of X i , the expression of the optimal solution X i of this sub-problem is:

[0063] X i = G i diag(σ1, σ2, …, σ m )V i H

[0064] Among them, diag(σ1, σ2, …, σ m ) represents a diagonal matrix with σ1, σ2, …, σ m as diagonal elements.

[0065] (3c) Fix X i , B i and β, and solve the sub-problem regarding the image x:

[0066]

[0067] This sub-problem is a large-scale least squares problem, and its approximate solution can be quickly calculated according to the following steps:

[0068] (3c1) Introduce an intermediate image variable z and a small constant τ > 0. The calculation formula for the intermediate image variable z is:

[0069]

[0070] Among them, Ι is the identity matrix, is the conjugate transpose matrix of the extraction matrix R i,j , x last represents the reconstructed image in the previous alternating direction multiplier method iteration, (X i +B i / β) j represents the j-th column of the matrix X i +B i / β;

[0071] (3c2) After obtaining the intermediate image variable z, the expression for the approximate solution x of this sub-problem is:

[0072] x = F H (λU H U + βτΙ) -1 (λU H y + βτFz)

[0073] Among them, U H is the conjugate transpose matrix of the undersampling matrix U, and F H is the conjugate transpose matrix of the orthogonal Fourier transform matrix F.

[0074] Step 4, after obtaining the optimal solutions of the variables in the reconstruction model, the update formulas for the Lagrange multiplier B i and the penalty parameter β are:

[0075]

[0076] β = γβ

[0077] where γ > 1 is the growth factor, and steps (3) to (4) are repeated until the reconstructed image meets the condition or the number of iterations reaches the preset upper limit.

[0078] The effects of the present invention can be further illustrated by the following simulation experiments:

[0079] I. Experimental conditions and content

[0080] Experimental conditions: A two-dimensional random sampling template with a sampling rate of 10% is used in the experiment; the experimental image adopts a real water balloon MRI image as Figure 2 shown; the peak signal-to-noise ratio (PSNR) and the structural similarity (SSIM) are used as objective evaluation indicators for the reconstruction results, and their definitions are as follows:

[0081]

[0082]

[0083] where x is the original image, is the reconstructed image, MAX is the maximum pixel value in the original image, N is the number of pixels in the image, μ x and are the means of x and respectively, σ x and are the standard deviations of x and respectively, is the covariance between x and , and C1 and C2 are constants to ensure the numerical stability of SSIM. The higher the values of PSNR and SSIM, the better the quality of the reconstructed image.

[0084] Experimental content: Under the above experimental conditions, the representative FDLCP method, PANO method, NLR method in the field of CS-MRI image reconstruction and the method of the present invention are selected for comparison.

[0085] Experiment 1: The K-space data of the MRI image shown in Figure 2 is reconstructed by the method of the present invention, the FDLCP method, the PANO method and the NLR method respectively. Among them, the FDLCP method realizes the sparse representation of each target image block by learning an orthogonal dictionary from various classified target image blocks and using the l1 norm as the sparse constraint term, and its reconstruction result is Figure 3 (a), and the reconstruction error is Figure 4 (a); the PANO method performs a two-dimensional wavelet transform on similar image blocks respectively, and then performs a wavelet transform on the wavelet coefficients at the same sub-band positions of the similar image blocks to obtain three-dimensional wavelet coefficients and constrain their sparsity, and its reconstruction result is Figure 3 (b), and the reconstruction error is Figure 4(b); The NLR method forms a similar image patch matrix with the similar image patches found in the image, and uses the logarithm of the determinant of the matrix to constrain its low-rank property. Its reconstruction result is Figure 3 (c), and the reconstruction error is Figure 4 (c); In the experiment, the method of the present invention sets the image patch size The number of image patches m in the similar image patch matrix is 32, L a,ε The exponential factor a = 0.5 and the small constant ε = 10 in the fractional norm -3 , the regularization parameter λ of the K-space data fidelity term is 10 6 , the small constant τ = 0.01, the growth factor γ = 1.2 in the calculation of the approximate solution of x, and the maximum number of iterations is 100 times. The reconstruction result of the method of the present invention is Figure 3 (d), and the reconstruction error is Figure 4 (d).

[0086] Figure 3 The small rectangular area is the selected image content to be enlarged, and the large rectangular area is the enlarged image content. From Figure 3 The reconstruction results, it can be seen that the FDLCP method introduces a large amount of foggy artifacts in the background and smooth areas of the image; the PANO method also results in obvious foggy artifacts and the restored texture details are relatively blurred; the NLR method makes the whole smooth area turn white and causes over-smoothing in the texture area; the method of the present invention effectively suppresses the foggy artifacts, while maintaining the original gray level in the smooth area and obtaining high-contrast image details. From Figure 4 The reconstruction error graph, it can be seen that the reconstruction error of the method of the present invention is significantly smaller than that of other methods, so the corresponding reconstruction result is closest to the original image.

[0087] Table 1 PSNR indexes of different reconstruction methods

[0088] Image FDLCP method PANO method NLR method Method of the present invention Water polo diagram 26.38 27.92 29.64 31.65

[0089] Table 1 gives Figure 3 the PSNR index situation of the reconstruction results of each method in, it can be seen that the PSNR value of the method of the present invention is greatly improved compared with other methods, indicating that the reconstruction quality obtained by the method of the present invention is the best, and this result is consistent with the reconstruction effect diagram.

[0090] Table 2 SSIM indexes of different reconstruction methods

[0091] Image FDLCP method PANO method NLR method Method of the present invention Water polo diagram 0.6402 0.7434 0.7362 0.7695

[0092] Table 2 gives Figure 3 the SSIM index situation of the reconstruction results of each method in, it can be seen that the SSIM value of the method of the present invention is the highest, indicating that the reconstruction result of the method of the present invention has a high degree of approximation to the original image, and this result is consistent with the reconstruction effect diagram.

[0093] The above experiments show that the reconstructed images obtained by the method of the present invention visually remove a large number of artifacts and restore rich detail information, and the corresponding objective evaluation indexes are better than those of the comparative method. Therefore, it can be seen that the present invention is effective for medical image reconstruction.

Claims

1. A CS-MRI image reconstruction method based on fractional norm low-rank constraint, comprising the following steps: (1) Input the K-space data y of an MRI and the sampling template, and perform pre-reconstruction using a traditional reconstruction method to obtain the initial reconstructed image x (0) ; (2) For each target image block of a certain size, a rectangular area of a fixed size centered on it is delimited in the initial reconstructed image x (0) respectively, and m image blocks with the smallest Euclidean distance from it are searched within this area, and a similar image block matrix including the target image block itself is formed. denotes the i-th similar image block matrix, where x is the image to be reconstructed, and R i,j is the extraction matrix of the j-th similar image block in X i , and is the extraction operator of the entire similar image block matrix. (3) Use the L a,ε fractional norm to constrain its low-rank property, representing X i 's L a,ε fractional norm, and its definition formula is: where, σ k is the k-th singular value of X i , a > 0 is an exponential factor and its value is not greater than 1, and ε > 0 is a small constant; (4) On the basis of utilizing the low-rank property of the matrix of similar image patches with fractional norm constraints, a reconstruction model for the image and the matrix of similar image patches is established: a,ε ​ Among them, represents the square of the two-norm of the vector, F u = UF is the undersampled Fourier encoding matrix, U is the undersampling matrix, F is the orthogonal Fourier transform matrix, and λ is the regularization parameter of the K-space data fidelity term; (5) Introduce the Lagrange multiplier B i and the penalty parameter β, and solve each optimization variable in the reconstruction model by the alternating direction method of multipliers: (5a) Fix x and B i and β, and solve the sub-problem regarding the similar image patch matrix X i : Among them, represents the square of the Frobenius norm of the matrix; (5b) Fix X i , B i and β, solve the sub-problem for the image x: After obtaining the optimal solutions of the variables in the reconstructed model, regarding the Lagrange multiplier B i and the penalty parameter The update formula for β is: β = γβ where γ > 1 is the growth factor, and steps (5) to (6) are repeated until the reconstructed image meets the conditions or the number of iterations reaches the preset upper limit.

2. A CS-MRI image reconstruction method based on fractional norm low-rank constraint according to claim 1, characterized in that Solving the sub-problem of the similar image patch matrix X in step (5a) i can be carried out according to the following steps: (5a1)Perform singular value decomposition on the matrix to obtain: Among them, G i is the left singular vector matrix, V i is the right singular vector matrix, δ1, δ2, …, δ m are the singular values of the matrix , and diag(δ1, δ2, …, δ m ) represents the diagonal matrix with δ1, δ2, …, δ m as diagonal elements; (5a2) Iteratively solve the similar image patch matrix X i for all singular values of X i For the k-th singular value σ k of X, let its initial value be and calculate the iterative formula: v = v + 1 where v represents the number of iterations, represents the calculated value of σ in the v-th iteration, k and this iterative formula is repeatedly calculated until v reaches a preset upper limit; (5a3) After obtaining the calculated value of the final iteration then judge The discriminant condition for whether it is the optimal solution is: When meets the discrimination condition, let When does not meet the discrimination condition, let σ k = 0; (5a4) After obtaining all the singular values of X i , the expression of the optimal solution X i of this sub-problem is: Among them, diag(σ1, σ2, …, σ m ) represents a diagonal matrix with σ1, σ2, …, σ m as diagonal elements.

3. A CS-MRI image reconstruction method based on fractional norm low-rank constraint according to claim 1, characterized in that, In step (5b), to solve the sub-problem regarding the image x, the following steps can be taken: (5b1) This sub-problem is a large-scale least squares problem. To quickly calculate its approximate solution, an intermediate image variable z and a small constant τ > 0 are introduced. The calculation formula for the intermediate image variable z is: where \(I\) is the identity matrix, is the extraction matrix \(R\) i,j of the conjugate transpose matrix, \(x\) last represents the reconstructed image in the previous alternating direction method of multipliers iteration, \((X\) i +B i / \(\beta)\) j represents the \(j\)-th column of the matrix \(X\) i +B i / \(\beta\); (5b2) After obtaining the intermediate image variable z, the expression for the approximate solution x of this sub-problem is: x = F H (λU H U + βτΙ) -1 (λU H y + βτFz) Among them, U H is the conjugate transpose matrix of the undersampling matrix U, and F H is the conjugate transpose matrix of the orthogonal Fourier transform matrix F.

Images