Method for sparse reconstruction of magnetic resonance images based on multi-dimensional proximal space total variation

Abstract

The application discloses a kind of magnetic resonance image sparse reconstruction methods based on multidimensional adjacent space total variation, it includes the following steps: step 1, data acquisition and pretreatment;Step 2, image reconstruction and gradient calculation: using data to carry out inverse Fourier transform, obtain reconstruction matrix M recovery , then calculate the gradient of reconstruction image;Step 3, the establishment of MPSTV regularization constraint: construct regularization term ‖m‖ MPSTV Based on MPSTV;Step 4, optimization and weight adjustment of constraint condition: in combination with the characteristics of data and gradient information in reconstruction process, by dynamically adjusting the weight λ1 And λ2 Of regularization term;Step 5, the solution of objective function and image reconstruction: construct optimization objective function, and solve the reconstruction of magnetic resonance image m by minimizing the function.This application not only solves the problem of step effect caused by traditional total variation method, but also improves the quality and precision of reconstruction image, provides a new solution for the fast reconstruction of magnetic resonance image.

Description

Technical Field

[0001] This invention belongs to the field of image reconstruction technology, and relates to fast magnetic resonance imaging reconstruction technology, specifically a sparse reconstruction method for magnetic resonance images based on multidimensional neighborhood space total variation. Background Technology

[0002] Compressed sensing (CS) is a signal processing technique designed to recover the original signal from a small amount of measurement data. It leverages the sparsity and redundancy of signals to explore their underlying structure and content globally, sampling signals at a rate far lower than traditional sampling theories (such as the Nyquist-Shannon sampling theorem) while maintaining low information loss, thus achieving more efficient signal acquisition and reconstruction. Therefore, CS is widely used in fast magnetic resonance imaging (MRI) reconstruction. It not only accelerates the acquisition process of traditional and parallel MRI data but can also be combined with deep learning algorithms as a prior model for deep learning-based MRI reconstruction. Current reconstruction algorithms often use Total Variation (TV) as a regularization term to smooth the reconstructed image. However, traditional TV methods primarily focus on the variational information of the horizontal and vertical sub-bands of the image, which may lead to a staircase effect in the reconstructed image; the TV regularization weights are often manually adjusted empirically, lacking objectivity and limiting the reconstruction results. Therefore, how to effectively mitigate or eliminate the step effect while maintaining the advantages of TV, and how to achieve automatic and forward-looking selection of regularization weights, has become a research hotspot in the field of fast 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 by proposing a sparse reconstruction method for magnetic resonance images based on Multidimensional Proximity Space Total Variation (MPSTV). This method is based on 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 technology fusion, while optimizing the image detail and signal-to-noise ratio.

[0004] To achieve the above objectives, the technical solution adopted by this invention is as follows: a sparse reconstruction method for magnetic resonance images based on multidimensional neighborhood space total variation, comprising the following steps:

[0005] Step 1, Data Acquisition and Preprocessing: 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 using the filtered K-space data to obtain the reconstruction matrix M. recovery Next, the gradient of the reconstructed image is calculated, and the feature information of the image is extracted using singular value decomposition.

[0007] Step 3: Establishing MPSTV regularization constraints: Constructing the MPSTV-based regularization term ‖m‖ MPSTV This term comprehensively considers variational information from different directions:

[0008]

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

[0010] Step 4: Optimization of constraints and adjustment of weights: Combining the characteristics of K-space data and gradient information in the reconstruction process, the weights of the regularization term are dynamically adjusted to ensure that edge details are preserved while suppressing image noise.

[0011] Step 5, Solving the objective function and image reconstruction: Constructing an optimization objective function and solving for the reconstruction of the magnetic resonance image m by minimizing this function:

[0012]

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

[0014] Furthermore, in step 1, firstly, K-space data f is acquired using a magnetic resonance imaging device. This data contains the frequency domain information of the image. Then, a Gaussian filter G(x,y) is applied to perform low-pass filtering on the K-space data to separate low-frequency and high-frequency components, and the filtered K-space data K is calculated. filtered (x, y) is used for subsequent image reconstruction:

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

[0016] Where K(x,y) is the original K-space data, It is a Gaussian filter.

[0017] Furthermore, in step 2, the reconstruction matrix M recovery Represented as:

[0018]

[0019] Furthermore, the threshold component T of the image gradient is calculated. low and T high :

[0020]

[0021] Where max(|M recovery |) represents the maximum absolute value of a pixel in the reconstructed image;

[0022] Next, the image reconstruction matrix M mentioned above is used. 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 with diagonal elements σ i For the singular values ​​of the image, extract these singular values ​​to form a vector σ = {σ1, σ2, σ3, ..., σ...} n The Sobel operator, combined with singular value weighting, is used to compute the gradients of the image in the x and y directions:

[0025]

[0026] Where, σ x and σ y These are singular values ​​that are related to the x and y directions.

[0027] Furthermore, using gradient components, the gradient magnitude G of the image is calculated:

[0028]

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

[0030] Furthermore, step 4 involves setting a low-frequency T. low and high frequency T high To achieve this:

[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 for edge regions, used to preserve image details. smooth λ is the regularization weight for smooth regions, used to suppress image noise. transition These are the regularization weights for the transition region, used to ensure a smooth image transition; T low and T high The low-frequency and high-frequency threshold components are obtained by Gaussian filtering in K-space.

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

[0034] The beneficial effects of this invention are as follows: By increasing MPSTV constraints and optimizing the weights of MPSTV regularization terms, this invention significantly improves the reconstruction accuracy of magnetic resonance images, reduces image information loss and staircase effect, effectively reduces the staircase effect in traditional TV constraints, further enhances image reconstruction performance, retains more global and local information of the reconstructed image, and significantly improves the accuracy and economic benefits of hospital diagnosis. Furthermore, this invention is suitable for rapid reconstruction needs in practical applications, significantly improves computational efficiency, reduces patient waiting time after MRI examinations, and promotes related scientific research, driving the advancement of medical imaging technology and having significant implications for the development of the medical industry. Attached Figure Description

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

[0036] Figure 2 The image shown is a test image from an example embodiment.

[0037] Figure 3 The image shown is a k-space undersampling trajectory diagram when the imaging acceleration factor AF = 5 in the embodiment. Detailed Implementation

[0038] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments.

[0039] A sparse reconstruction method for magnetic resonance images based on multidimensional neighborhood space total variation is proposed. This method achieves higher reconstruction accuracy by adding constraints to the magnetic resonance image reconstruction process, and reduces the staircase effect by adding variational information from two oblique sub-bands within the same layer and variational information from the vertical projection positions of adjacent layers, thereby improving the reconstruction performance of magnetic resonance imaging. Experiments show that when the three variational pairs are isotropic internally and anisotropic between each other, the staircase effect can be effectively reduced, effectively achieving the expected goal. The specific details of this method are as follows:

[0040] 1. Theory and Methodology

[0041] (1) Sparse reconstruction model of magnetic resonance images

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

[0043]

[0044] This formula integrates three objectives in the optimization process of reconstructing the nuclear magnetic resonance image m, namely, through the MPSTV regularization term ‖m‖ MPSTV To preserve the edge details of the image, a regularization term l1, ||Φm||1, is used to promote the sparsity of the solution in the transform domain, and then an error term l2 is used. Ensure that the reconstructed image has the minimum difference from the observed data f.

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

[0046] The reconstructed magnetic resonance image can be obtained by solving equation (1).

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

[0048] Traditional total variation constraints, which only include variational information in the x and y directions, often lead to staircase effects. Based on this and by exploring the characteristics of magnetic resonance images, a novel total variation norm constraint is designed to reduce staircase effects, thereby further improving the performance of multi-constraint model magnetic resonance reconstruction algorithms. The specific implementation is as follows:

[0049] Due to the global correlation of reconstructed magnetic resonance images, from the perspective of multidimensional proximity correlation space, the value of each voxel in the magnetic resonance image is similar to the values ​​in its neighborhood. This implies that the magnetic resonance image exhibits a locally smooth structural pattern relationship across both temporal and spatial scales. Furthermore, most differences between adjacent frames in the temporal domain are zero; these can be considered statistical differences in the magnetic resonance image, suggesting considerable sparsity in the gradient domain. Based on these characteristics of reconstructed magnetic resonance images, we have invented a magnetic resonance image reconstruction model that integrates multidimensional conjugate spatial factors and temporal factors, based on the similarity of neighboring magnetic resonance images—a sparse reconstruction method for magnetic resonance images 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: the weight coefficient ω and the gradient coefficient. The weighting coefficients in the formula are mainly used to control the contribution of gradients in each direction to the total variation. Different weights reflect the importance assessment 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, k in different spatial dimensions.

[0054] This formula comprehensively considers the gradient changes of 3D data in different directions and calculates the total variational value through weighted summation. The formula primarily reduces the staircase effect by increasing the variational information of the two oblique sub-bands within the same layer and the variational information of the vertical projection positions of adjacent layers. i, j, and k represent rows, columns, and different layers within the same layer, respectively. ω1 is the weighting parameter used to calculate the variational information in the vertical and horizontal directions within the same layer; ω2 is the weighting parameter used to calculate the variational information in the two oblique directions within the same layer; and ω3 is the weighting parameter used to calculate the variational information in the vertical projection position directions of different layers. Furthermore, our preliminary experiments show that when the three variational pairs are isotropic within each pair and anisotropic between each pair, the staircase effect can be effectively reduced, improving reconstruction performance.

[0055] (3) MPSTV Regularization Term Weight Optimization Method

[0056] In MRI, K-space represents the frequency domain information of an image. To better balance edge detail preservation and noise suppression during image reconstruction, the applicant employed a method that uses a Gaussian filter to determine the frequency threshold component T.low and T high It combines image gradient and Singular Value Decomposition (SVD) with dynamically adjusted regularization weights λ. ε The method is used to optimize the MPSTV regularization term weights. For the MPSTV regularization term λ1‖m‖ in formula (1) MPSTV , 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, separating the low-frequency and high-frequency components. The formula is as follows:

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

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

[0060]

[0061] Here, σ determines the bandwidth of the filter. Low-frequency components represent the global information of the image, while high-frequency components represent the details and edge information. The total energy of the K-space data after passing through the Gaussian filter is then calculated, and a threshold is set based on a certain proportion of this energy, thereby determining the threshold components T for low and high frequencies. low and T high The energy E of the filtered K-space data filtered Threshold energy E of low and high frequencies low and E high Represented by the following formulas respectively:

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

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

[0064] Here, α and β are proportional parameters used to control the low-frequency and high-frequency thresholds, respectively (0 < α < β < 1). The corresponding image reconstruction matrix M is then 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] Where max(|M recovery |) represents the maximum absolute value of the pixel values ​​in the reconstructed image. Next, using the image reconstruction matrix M described above... 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 with diagonal elements σ i For the singular values ​​of the image, extract these singular values ​​to form a vector σ = {σ1, σ2, σ3, ..., σ...} n The Sobel operator, combined with singular value weighting, is used to compute the gradients of the image in the x and y directions:

[0069]

[0070] Where, σ x and σ y These are singular values ​​related to the x and y directions. After calculation using the above formula, each component of the image gradient incorporates the influence of singular values, allowing for a more comprehensive consideration of the image's feature information in different directions.

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

[0072]

[0073] The gradient magnitude G incorporates singular value information and reflects the intensity changes at each pixel in the image, improving the accuracy of edge features in different directions. Subsequently, based on the calculated gradient magnitude G and low-frequency / high-frequency threshold T... low and T high and singular value σ i For regularization weight λ ε Adjustments have been made, and the specific formula is as follows:

[0074]

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

[0076] In this way, by combining the image gradient calculation method with singular values, the allocation of regularization weights can be effectively optimized, thereby better balancing the needs of edge detail preservation and noise suppression during image reconstruction.

[0077] 2. Experimental Verification

[0078] To test the performance of the proposed algorithm, the applicant selected 176 consecutive MRI images of the head. Images from the 60th and 140th slices 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, downsampled to 256×256. Test images are shown below. Figure 2 As shown (from left to right, these are the MR images of layer 60 and layer 140, respectively).

[0079] To better compare the reconstruction performance of the algorithms, the applicant selected a 5x imaging acceleration factor, i.e., AF = 5. The k-space undersampled trajectory when AF = 5 is as follows: Figure 3 As shown.

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

[0081] Table 1. Performance comparison of MR image data from the 60th and 140th layers of the head at 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 λ. ε Compare them.

[0084] Table 2. Comparison of regularization weight performance during the reconstruction of MR images at layers 60 and 140 of the head with AF=5.

[0085]

[0086] The foregoing has shown and described 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 scope of protection of the present invention in any way, and all technical solutions obtained by equivalent substitution or other means fall within the scope of protection of the present invention. Parts not covered in this invention are the same as or can be implemented using existing technology.

Claims

1. A sparse reconstruction method for magnetic resonance images based on multidimensional neighborhood space total variation, characterized in that, Includes the following steps: Step 1, Data Acquisition and Preprocessing: 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 using the filtered K-space data to obtain the reconstruction matrix M. recovery Next, the gradient of the reconstructed image is calculated, and the feature information of the image is extracted using singular value decomposition. Step 3: Establishing MPSTV regularization constraints: Constructing the MPSTV-based regularization term ‖m‖ MPSTV This term comprehensively considers variational information from different directions: in, Represents the gradient of the image in the plane. This represents the gradient in the diagonal direction. ω1 represents the gradient between vertical layers; i, j, k represent the row, column and different layers of the same layer, 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 diagonal directions of the same layer; and ω is the weight parameter used to calculate the variational information in the vertical projection position direction of different layers. Step 4: Optimization of constraints and adjustment of weights: Combining the characteristics of K-space data and gradient information in the reconstruction process, the weights of the regularization term are dynamically adjusted to ensure that edge details are preserved while suppressing image noise. Step 5, Solving the objective function and image reconstruction: Constructing an optimization objective function and solving for the reconstruction of the magnetic resonance image m by minimizing this function: Where Φ is the discrete wavelet transform matrix; ||m|| MPSTV For MPSTV regularization, ||Φm||1 is for l1 regularization. Let λ1 be the weight of the l2 error term, λ2 be the weight of the MPSTV regularization term, and λ3 be the weight of the l1 regularization term. The model's error term is defined as follows: μ 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. The MPSTV regularization term λ1‖m‖ in the formula is also included. MPSTV Further, ∑λ ε ||m ε || MPSTV ; In step 2, the reconstruction matrix M recovery Represented as: The threshold component T of the image gradient is obtained through calculation. low and T high : Where max(|M recovery |) is the maximum absolute value of the pixel value in the reconstructed image; α and β are the proportional parameters used to control the low-frequency and high-frequency thresholds, respectively; Next, the image reconstruction matrix M mentioned above is used. 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 with diagonal elements σ i For the singular values ​​of the image, extract these singular values ​​to form a vector σ = {σ1, σ2, σ3, ..., σ...} n The Sobel operator, combined with singular value weighting, is used to compute the gradients of the image in the x and y directions: Where, σ x and σ y These are singular values ​​that relate to the x and y directions; Calculate the gradient magnitude G of the image using gradient components: Among them, G x and G y These represent the gradients in the x and y directions, respectively. Step 4 involves setting a low frequency T. low and high frequency T high To achieve this: Among them, G ε λ is the gradient magnitude of the ε-th pixel. ε λ is the regularization weight of the ε-th pixel; edge λ is the regularization weight for edge regions, used to preserve image details. smooth λ is the regularization weight for smooth regions, used to suppress image noise. transition These are the regularization weights for the transition region, used to ensure a smooth image transition; T low and T high The low-frequency and high-frequency threshold components are obtained by Gaussian filtering in K-space.

2. The sparse reconstruction method for magnetic resonance images based on multidimensional neighborhood space total variation as described in claim 1, characterized in that, In step 1, firstly, K-space data f is acquired using a magnetic resonance imaging (MRI) device. This data contains the frequency domain information of the image. Then, a Gaussian filter G(x,y) is applied to the K-space data for low-pass filtering to separate low-frequency and high-frequency components, and the filtered K-space data K is calculated. filtered (x, y) is used for subsequent image reconstruction: K filtered (x,y)=K(x,y)·G(x,y) Where K(x,y) is the original K-space data, It is a Gaussian filter.

3. The sparse reconstruction method for magnetic resonance images based on multidimensional neighborhood space total variation as described in claim 1, characterized in that, In step 5, the optimization process for reconstructing the MRI image m integrates three objectives, namely, through the MPSTV regularization term ‖m‖ MPSTV To preserve the edge details of the image, a regularization term l1, ||Φm||1, is used to promote the sparsity of the solution in the transform domain, and then an error term l2 is used. Ensure that the reconstructed image has the minimum difference from the observed data f.

Images