Magnetic resonance image sparse reconstruction method based on multi-dimensional adjacent space total variation

Abstract

The invention discloses a magnetic resonance image sparse reconstruction method based on multi-dimensional adjacent space total variation. The method comprises the following steps: step 1, data acquisition and preprocessing; step 2, performing image reconstruction and gradient calculation: performing inverse Fourier transform by using the data to obtain a reconstruction matrix Mrecovery, and then calculating the gradient of a reconstructed image; step 3, establishment of MPSTV regularization constraint: constructing a regularization item | | m | | MPSTV based on the MPSTV; step 4, optimization of constraint conditions and weight adjustment: combining data characteristics and gradient information in a reconstruction process, and dynamically adjusting weights lambda 1 and lambda 2 of a regularization item; and step 5, objective function solving and image reconstruction: constructing an optimized objective function, and solving reconstruction of the magnetic resonance image m by minimizing the function. According to the method, the step effect problem caused by a traditional total variation method is solved, meanwhile, the quality and precision of the reconstructed image are improved, and a new solution is provided for rapid reconstruction of the magnetic resonance image.

Description

Technical Field

[0001] The invention belongs to the technical field of image reconstruction, and relates to a fast magnetic resonance imaging reconstruction technology, in particular to a magnetic resonance image sparse reconstruction method based on multi-dimensional neighboring space total variation. Background Art

[0002] Compressed sensing (CS) is a signal processing technique that aims to recover the original signal from a small amount of measurement data. It uses the sparsity and redundancy of the signal to explore the basic structure and content of the signal from a global perspective, and samples the signal at a rate far lower than the traditional sampling theory (such as the Nyquist–Shannon sampling theorem) while ensuring a low information loss rate, thereby performing more efficient signal acquisition and reconstruction. Therefore, CS is also widely used in the field of fast magnetic resonance imaging reconstruction. It can not only accelerate the acquisition process of traditional and parallel magnetic resonance data, but also be combined with deep learning algorithms and applied as a priori models in the field of deep learning magnetic resonance imaging reconstruction. Current reconstruction algorithms often use total variation (TV) as a regularization term to smooth the reconstructed image. However, the traditional total variation method mainly focuses on the horizontal and vertical sub-band variation information of the image, which may lead to a step effect in the reconstructed image; the total variation regularization weights are often manually adjusted based on experience, which lacks objectivity and restricts the reconstruction results. Therefore, how to effectively reduce or eliminate the step effect while maintaining the advantages of TV and realize automatic and prospective selection of regularization weights has become a research hotspot in the field of rapid magnetic resonance imaging reconstruction in recent years. Summary of the invention

[0003] The purpose of this patent is to address the problems existing in the background technology and propose a sparse reconstruction method for magnetic resonance images based on multidimensional proximity space total variation (MPSTV). This method is based on the compressed sensing sparse reconstruction theory, constructs a multidimensional proximity space total variation strategy, and dynamically adjusts the regularization weights. It aims to improve the reconstruction speed and quality of magnetic resonance images through technical fusion, while optimizing the image detail performance and signal-to-noise ratio.

[0004] To achieve the above object, the technical solution adopted by the present invention is: a sparse reconstruction method of magnetic resonance images based on multi-dimensional neighboring space total variation, which comprises the following steps:

[0005] Step 1: Data acquisition and preprocessing: The filtered K-space data K is obtained through acquisition and processing. filtered (x,y);

[0006] Step 2: Image reconstruction and gradient calculation: Perform inverse Fourier transform on the filtered K-space data to obtain the reconstruction matrix M recovery ,Then, the gradient of the reconstructed image is calculated, and the feature information of the image is extracted using singular value decomposition;

[0007] Step 3. Establishment of MPSTV regularization constraints: Construct the regularization term ‖m‖ based on MPSTV MPSTV , which comprehensively considers the variational information in different directions:

[0008]

[0009] in, represents the gradient of the image in the plane, represents the gradient in the diagonal direction, represents the gradient between vertical layers;

[0010] Step 4: Optimization of constraint conditions and weight adjustment: Based on the characteristics of K-space data and the gradient information during reconstruction, the weights λ1 and λ2 of the regularization term are dynamically adjusted to ensure that image noise is suppressed while retaining edge details.

[0011] Step 5: Solving the objective function and image reconstruction: Construct an optimization objective function and solve the reconstruction of the magnetic resonance image m by minimizing the function:

[0012]

[0013] where Φ is the discrete wavelet transform matrix; specifically, ‖m‖ MPSTV is the MPSTV regularization term, ‖Φm‖1 is the l1 regularization term, is the l2 error term, λ1 is the weight of the MPSTV regularization term, λ2 is the weight of the l1 regularization term, is the error term of the reconstruction model, μ is the error term weight, R is the K-space phase encoding matrix, F is the two-dimensional Fourier transform matrix, and f is the undersampled K-space data.

[0014] Furthermore, in step 1, first, the K-space data f is obtained by a magnetic resonance imaging device, and the data contains the frequency domain information of the image; then, a Gaussian filter G(x, y) is applied to low-pass filter the K-space data to separate the low-frequency and high-frequency components, and the filtered K-space data K is calculated. filtered (x,y) for subsequent image reconstruction:

[0015] K filtered (x,y)=K(x,y)·G(x,y)

[0016] Among them, K(x,y) is the original K space data, is a Gaussian filter.

[0017] Furthermore, in step 2, the matrix M is reconstructed recovery It is expressed as:

[0018]

[0019] Furthermore, the threshold component T of the image gradient is obtained by calculation low and T high :

[0020]

[0021] Among them, max(|M recovery |) is the maximum absolute value of the pixel value in the reconstructed image;

[0022] Secondly, the image reconstruction matrix M recovery , perform SVD:

[0023] M recovery =USV T

[0024] Where U is the left singular vector matrix, V is the right singular vector matrix, and S is a diagonal matrix whose diagonal elements σ i are the singular values ​​of the image, and these singular values ​​are extracted to form a vector σ={σ1,σ2,σ3,……,σ n}, and use the Sobel operator combined with singular value weighting to calculate the gradient of the image in the x and y directions:

[0025]

[0026] Among them, σ x and σ y are the singular values ​​associated with the x and y directions.

[0027] Furthermore, the gradient component is used to calculate the gradient magnitude G of the image:

[0028]

[0029] Among them, G x and G y Represent the gradients in the x and y directions respectively.

[0030] Furthermore, the step 4 is performed by setting the low frequency T low and high frequency T high Threshold to achieve:

[0031]

[0032] Among them, Gε is the gradient magnitude of the εth pixel, λ ε is the regularization weight of the ε-th pixel; λ edge is the regularization weight of the edge region, used to preserve image details, λ smooth is the regularization weight of the smooth region, used to suppress image noise, λ transition is the regularization weight of the transition region, which is used to ensure smooth transition of the image; T low and T high It is the low-frequency and high-frequency threshold components calculated by Gaussian filtering in K space.

[0033] Furthermore, in step 5, three objectives are integrated in the optimization process of reconstructing the MRI image m, namely, the MPSTV regularization term ‖m‖ MPSTV To maintain the edge details of the image, the l1 regularization term ‖Φm‖1 is used to promote the sparsity of the solution in the transform domain, and then the l2 error term Ensure that the reconstructed image has minimal difference from the observed data f.

[0034] The beneficial effects of the present invention are as follows: by adding MPSTV constraints and optimizing the weight of the MPSTV regularization term, the present invention significantly improves the reconstruction accuracy of magnetic resonance images, reduces image information loss and step effects, effectively reduces the step effects in traditional TV constraints, further improves image reconstruction performance, retains more global and local information of reconstructed images, and significantly improves hospital diagnosis accuracy and economic benefits. In addition, the present invention is suitable for the needs of rapid reconstruction in practical applications, significantly improves computational efficiency, and reduces the waiting time of patients after MRI examinations; it also promotes the development of related scientific research, promotes the advancement of medical imaging technology, and is of great significance to the development of the medical industry. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 This is a structural diagram of the present invention reconstructed based on the MPSTV magnetic resonance model;

[0036] Figure 2 This is a test image in the embodiment;

[0037] Figure 3 1 is a k-space undersampling trajectory diagram when the imaging acceleration factor AF=5 in the embodiment. DETAILED DESCRIPTION

[0038] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0039] A sparse reconstruction method for magnetic resonance images based on total variation of multi-dimensional neighboring space is proposed. It obtains higher reconstruction accuracy by adding constraints in magnetic resonance image reconstruction, and reduces the staircase effect by increasing the variational information of two oblique subbands in the same layer and the variational information of the vertical projection position of adjacent layers to improve the reconstruction performance of magnetic resonance imaging. Experiments show that when the three variational pairs are isotropic internally and anisotropic between the three variational pairs, the staircase effect can be effectively reduced and the expected goal can be effectively achieved. The specific contents of this method are as follows:

[0040] 1. Theory and Methods

[0041] (1) Magnetic resonance image sparse reconstruction model

[0042] The basic mathematical model of sparse reconstruction of magnetic resonance images based on compressed sensing is as follows:

[0043]

[0044] This formula combines three objectives in the optimization process of reconstructing the MRI image m, namely, through the MPSTV regularization term ‖m‖ MPSTV To maintain the edge details of the image, the l1 regularization term ‖Φm‖1 is used to promote the sparsity of the solution in the transform domain, and then the l2 error term Ensure that the reconstructed image has minimal difference from the observed data f.

[0045] This formula has two constraints, where λ1 and λ2 are the weights of the MPSTV regularization term and the l1 regularization term, respectively. Φ represents the discrete wavelet transform matrix, which can be used to perform l1 regularization on m in the discrete wavelet transform domain. represents the error term of the reconstruction model, μ represents the weight of the error term, which can control the influence of the error term in the overall objective function, R represents the K-space phase encoding matrix, F represents the two-dimensional Fourier transform matrix, and f represents the undersampled K-space data.

[0046] By solving equation (1), the reconstructed magnetic resonance image can be obtained.

[0047] (2) Multidimensional Neighborhood Space Total Variation Optimization Method

[0048] In the traditional total variation constraint, since only the variation information in the x and y directions is included, the staircase effect often occurs. Based on this and exploring the characteristics of magnetic resonance images, a new total variation norm constraint is designed to reduce the staircase effect, thereby further improving the performance of the magnetic resonance reconstruction algorithm of the multi-constraint model. The specific implementation method is as follows:

[0049] Since there is global correlation in MRI reconstruction images, from the perspective of multidimensional neighboring correlation space, the value of each voxel in the MRI is similar to the value in its neighborhood, which means that there is a local smooth structural pattern relationship in the MRI in terms of time and space. In addition, most of the differences between adjacent frames in the time domain are 0, which can be considered as the difference information statistics in the MRI, which means that there is considerable sparsity between the gradient domains of the MRI. Based on the above characteristics of MRI reconstruction images, we invented a MRI reconstruction model based on the similarity of neighboring MRI images, which integrates multidimensional conjugate spatial factors and time factors, that is, a sparse reconstruction method for MRI based on multidimensional proximity space total variation (MPSTV).

[0050] ‖m‖ MPSTV The core idea can be expressed as follows:

[0051]

[0052]

[0053] This formula has two main parts, namely the weight coefficient ω and the gradient coefficient The weight coefficient in the formula is mainly used to control the contribution of gradients in various directions to the total variation. Different weights reflect the importance of gradients in different directions when calculating the total variation. The gradient norm is mainly used to represent the rate of change of data points i, j, and k in different spatial dimensions.

[0054] This formula comprehensively considers the gradient changes of three-dimensional data in different directions, and calculates the total variation value through weighted summation. The above formula mainly reduces the staircase effect by adding the variational information of the two oblique subbands in the same layer and the variational information of the vertical projection position of the adjacent layer. i, j, k represent the rows of the same layer, the columns of the same layer and the different layers respectively. ω1 is the weight parameter used to calculate the variational information in the vertical and horizontal directions of the same layer. ω2 is the weight parameter used to calculate the variational information in the two oblique directions of the same layer. ω3 is the weight parameter used to calculate the variational information in the vertical projection position direction of different layers. In addition, we found through preliminary experiments that when the three variational pairs are isotropic internally and anisotropic between the three variational pairs, the occurrence of staircase effect can be effectively reduced and the reconstruction performance can be improved.

[0055] (3) MPSTV regularization term weight optimization method

[0056] In MRI, K space represents the frequency domain information of the image. In order to better balance edge detail preservation and noise suppression during image reconstruction, the applicant adopts a method of determining the frequency threshold component T using a Gaussian filter.low and T high The regularization weight λ is adjusted dynamically based on image gradient and singular value decomposition (SVD). ε The method is used to optimize the MPSTV regularization term weight. For the MPSTV regularization term λ1‖m‖ in formula (1) MPSTV , which can be further expressed as ε ‖m ε ‖ MPSTV .

[0057] First, for the K-space data K(x,y), a Gaussian filter G(x,y) is applied to perform low-pass filtering to separate low-frequency and high-frequency components. The formula is as follows:

[0058] K filtered (x,y)=K(x,y)·G(x,y) (6)

[0059] Among them, G(x,y) is a standard Gaussian filter, defined as:

[0060]

[0061] Among them, σ determines the bandwidth of the filter. The low-frequency component represents the global information of the image, while the high-frequency component represents the details and edge information of the image. Then the total energy of the K-space data after passing through the Gaussian filter is calculated, and the threshold is set at a certain ratio of the energy to determine the threshold components T of the low-frequency and high-frequency. low and T high The energy E of the filtered K-space data filtered , low-frequency and high-frequency threshold energy E low and E high They are expressed in the following formulas:

[0062] E filtered =Σ x,y |K filtered (x,y)| 2 (8)

[0063] E low =α·E filtered ,E high =β·E filtered (9)

[0064] Among them, α and β are proportional parameters used to control the low-frequency and high-frequency thresholds respectively (0<α<β<1). Then, the corresponding image reconstruction matrix M is obtained by calculating the inverse Fourier transform of the filtered K-space data. recovery , and calculate the threshold component T of the image gradient low and Thigh :

[0065]

[0066] Among them, max(|M recovery |) is the maximum absolute value of the pixel value in the reconstructed image. recovery , perform SVD:

[0067] M recovery =USV T (12)

[0068] Where U is the left singular vector matrix, V is the right singular vector matrix, and S is a diagonal matrix whose diagonal elements σ i are the singular values ​​of the image, and these singular values ​​are extracted to form a vector σ={σ1,σ2,σ3,……,σ n}, and use the Sobel operator combined with singular value weighting to calculate the gradient of the image in the x and y directions:

[0069]

[0070] Among them, σ x and σ y It is the singular value related to the x and y directions. After calculation according to the above formula, each component of the image gradient combines the influence of the singular value, which can more completely consider the characteristic information of the image in different directions.

[0071] Using the above gradient components, calculate the gradient magnitude G of the image:

[0072]

[0073] Among them, the gradient amplitude G combines the singular value information and reflects the intensity change at each pixel of the image, which improves the accuracy of edge features in different directions of the image. Then, according to the calculated gradient amplitude G, low-frequency and high-frequency thresholds T low and T high And the singular values ​​σ i For the regularization weight λ ε Make adjustments, the specific formula is as follows:

[0074]

[0075] Among them, G ε is the gradient magnitude of the εth pixel, λ ε is the regularization weight of the ε-th pixel. edge is the regularization weight of the edge region (i.e., the region with a large gradient magnitude), which is used to preserve image details. smoothis the regularization weight of the smooth area (i.e. the area with small gradient amplitude), which is used to suppress image noise. transition is the regularization weight in the transition region, which is used to ensure smooth transition of the image. low and T high It is the low-frequency and high-frequency threshold components calculated by Gaussian filtering in K space.

[0076] In this way, the image gradient calculation method combined with singular values ​​can effectively optimize the distribution of regularization weights, thereby better balancing the needs of edge detail preservation and noise suppression during image reconstruction.

[0077] 2. Experimental demonstration

[0078] In order to test the performance of the proposed algorithm, the applicant selected 176 consecutive MRI images of the head. The 60th and 140th MR images were randomly selected for detailed analysis. The test dataset was acquired on a 3.0T Siemens Trio Tim scanner at Jiangsu Provincial People's Hospital. The original image resolution was 1376×1376, and the resolution was downsampled to 256×256. The test images are as follows: Figure 2 As shown (from left to right are the 60th and 140th layer MR images respectively).

[0079] In order to better compare the reconstruction performance of the algorithm, the applicant selected a 5-fold imaging acceleration factor, i.e., AF = 5. The k-space undersampling trajectory when AF = 5 is shown in Figure 3 shown.

[0080] In this embodiment, the applicant selected three evaluation parameters to compare the image reconstruction performance of different algorithms, including mean square error (MSE), peak signal-to-noise ratio (PSNR), and normalized mutual information (NMI).

[0081] Table 1. Comparison of performance of 60th and 140th MR image data of the head when AF=5

[0082]

[0083] In this embodiment, the applicant selected an empirical regularization weight λ1=0.01, which is different from the automatically determined regularization weight λ ε Make a comparison.

[0084] Table 2. Comparison of the performance of regularized weights during reconstruction of 60th and 140th head MR images when AF = 5

[0085]

[0086] The above shows and describes the basic principles, main features and advantages of the present invention. Those skilled in the art should understand that the above embodiments do not limit the protection scope of the present invention in any form, and all technical solutions obtained by equivalent replacement and the like fall within the protection scope of the present invention. The parts not involved in the present invention are the same as the prior art or can be implemented by the prior art.

Claims

1. A method for sparse reconstruction of magnetic resonance images based on multidimensional neighboring space total variation, characterized in that: The steps include: Step 1: Data acquisition and preprocessing: The filtered K-space data K is obtained through acquisition and processing. filtered (x,y); Step 2: Image reconstruction and gradient calculation: Perform inverse Fourier transform on the filtered K-space data to obtain the reconstruction matrix M recovery ,Then, the gradient of the reconstructed image is calculated, and the feature information of the image is extracted using singular value decomposition; Step 3. Establishment of MPSTV regularization constraints: Construct the regularization term ‖m‖ based on MPSTV MPSTV , which comprehensively considers the variational information in different directions: in, represents the gradient of the image in the plane, represents the gradient in the diagonal direction, represents the gradient between vertical layers; Step 4: Optimization of constraint conditions and weight adjustment: Based on the characteristics of K-space data and the gradient information during reconstruction, the weights λ1 and λ2 of the regularization term are dynamically adjusted to ensure that image noise is suppressed while retaining edge details. Step 5: Solving the objective function and image reconstruction: Construct an optimization objective function and solve the reconstruction of the magnetic resonance image m by minimizing the function: Where Φ is the discrete wavelet transform matrix; specifically, ‖m‖ MPSTV is the MPSTV regularization term, ‖Φm‖1 is the l1 regularization term, is the l2 error term, λ1 is the weight of the MPSTV regularization term, λ2 is the weight of the l1 regularization term, is the error term of the reconstruction model, μ is the error term weight, R is the K-space phase encoding matrix, F is the two-dimensional Fourier transform matrix, and f is the undersampled K-space data.

2. The method for sparse reconstruction of magnetic resonance images based on multi-dimensional neighboring space total variation according to claim 1, characterized in that: In step 1, first, the K-space data f is obtained by a magnetic resonance imaging device, and the data contains the frequency domain information of the image; then, a Gaussian filter G(x, y) is applied to low-pass filter the K-space data to separate the low-frequency and high-frequency components, and the filtered K-space data K is calculated. filtered (x,y) for subsequent image reconstruction: K filtered (x,y)=K(x,y)·G(x,y) Among them, K(x,y) is the original K space data, is a Gaussian filter.

3. The method for sparse reconstruction of magnetic resonance images based on multi-dimensional neighboring space total variation according to claim 1, characterized in that: In step 2, the matrix M is reconstructed recovery It is expressed as:

4. The method for sparse reconstruction of magnetic resonance images based on multi-dimensional neighboring space total variation according to claim 3, characterized in that: The threshold component T of the image gradient is obtained by calculation low and T high : Among them, max(|M recovery |) is the maximum absolute value of the pixel value in the reconstructed image; Secondly, the image reconstruction matrix M recovery , perform SVD: M recovery =USV T Where U is the left singular vector matrix, V is the right singular vector matrix, and S is a diagonal matrix whose diagonal elements σ i are the singular values ​​of the image, and these singular values ​​are extracted to form a vector σ={σ1,σ2,σ3,……,σ n }, and use the Sobel operator combined with singular value weighting to calculate the gradient of the image in the x and y directions: Among them, σ x and σ y are the singular values ​​associated with the x and y directions.

5. The method for sparse reconstruction of magnetic resonance images based on multi-dimensional neighboring space total variation according to claim 4, characterized in that: Using the gradient component, calculate the gradient magnitude G of the image: Among them, G x and G y Represent the gradients in the x and y directions respectively.

6. The method for sparse reconstruction of magnetic resonance images based on multi-dimensional neighboring space total variation according to claim 5, characterized in that: Step 4 is to set the low frequency T low and high frequency T high Threshold to achieve: Among them, G ε is the gradient magnitude of the εth pixel, λ ε is the regularization weight of the ε-th pixel; λ edge is the regularization weight of the edge region, used to preserve image details, λ smooth is the regularization weight of the smooth region, used to suppress image noise, λ transition is the regularization weight of the transition region, which is used to ensure smooth transition of the image; T low and T high It is the low-frequency and high-frequency threshold components calculated by Gaussian filtering in K space.

7. The method for sparse reconstruction of magnetic resonance images based on multi-dimensional neighboring space total variation according to claim 1, characterized in that: In step 5, three objectives are integrated in the optimization process of reconstructing the MRI image m, namely, through the MPSTV regularization term ‖m‖ MPSTV To maintain the edge details of the image, the l1 regularization term ‖Φm‖1 is used to promote the sparsity of the solution in the transform domain, and then the l2 error term Ensure that the reconstructed image has minimal difference from the observed data f.

Images