Rapid magnetic resonance image reconstruction method based on reweighted shrinkage threshold algorithm

Abstract

The invention discloses a rapid magnetic resonance image reconstruction method based on a reweighted shrinkage threshold algorithm. The method comprises the following steps: step 1, data acquisition and preprocessing; step 2, sparse transformation and coding matrix construction; step 3, constructing an optimization model; step 4, an iterative reconstruction process; and 5, evaluating and analyzing a result. According to the method, the problem that the balance between the convergence speed and the reconstruction precision is affected due to the fact that regularization parameters are difficult to adjust in real time according to the convergence state in iteration of traditional ISTA is solved; the image detail reduction capability is further enhanced by combining the RISTAC algorithm of contourlet transform, effective improvement on traditional wavelet transform is proved, and beneficial reference is provided for rapid MRI imaging research and practical application.

Description

Technical Field

[0001] The invention relates to a fast magnetic resonance image reconstruction method based on a reweighted shrinkage threshold algorithm, and belongs to the field of fast magnetic resonance image reconstruction. Background Art

[0002] The Nyquist-Shannon sampling theorem, a basic guiding principle for traditional image acquisition systems, dictates that the sampling frequency be at least twice the highest frequency component of the signal to ensure that the original signal can be reconstructed without distortion. As a result, traditional imaging methods often require capturing raw image data at excessively high sampling rates, a process that is both time-consuming and requires sophisticated equipment capabilities. At the same time, traditional magnetic resonance imaging (MRI) procedures are time-consuming, susceptible to motion-induced artifacts, and may lead to reduced patient compliance. In addition, in imaging methods that require the use of contrast agents, extended sampling times may reduce the effectiveness of the contrast agents, resulting in poor image quality or unclear diagnostic results.

[0003] To address these diagnostic challenges, the theory of compressed sensing (CS) has been introduced, which allows the reconstruction of images from sparse representations and randomly sampled data. Sparse representation techniques achieve sparse representation of images by transforming the original image into a suitable transform domain, where the image has only a few non-zero coefficients. This can be achieved by minimizing a cost function based on the transform domain. Compressed sampling is obtained by undersampling the original signal, which significantly reduces the required bandwidth and the burden of image acquisition.

[0004] In the field of contemporary fast magnetic resonance imaging, compressed sensing technology has become a mature method. The iterative shrinkage threshold algorithm (ISTA) has become a major agent-based algorithm due to its computational efficiency in reconstructing magnetic resonance images. However, traditional ISTA-type algorithms lack the ability to dynamically adjust the regularization parameters, which affects the convergence speed.

[0005] Therefore, the present application proposes a fast magnetic resonance image reconstruction method based on a reweighted shrinkage threshold algorithm, which corrects the regularization weight by dynamically adjusting the weights of two adjacent iterations, thereby improving the convergence efficiency while retaining the shrinkage characteristics. Summary of the invention

[0006] The purpose of the present invention is to propose a fast magnetic resonance image reconstruction method based on a reweighted shrinkage threshold algorithm to solve the problem that it is difficult for traditional ISTA to adjust the regularization parameters in real time according to the convergence state during iteration, which affects the balance between convergence speed and reconstruction accuracy.

[0007] To achieve the above object, the technical solution adopted by the present invention is: a fast magnetic resonance image reconstruction method based on a reweighted shrinkage threshold algorithm, which comprises the following steps:

[0008] Step 1: Data collection and preprocessing;

[0009] Step 2: Sparse transformation and coding matrix construction;

[0010] Step 3: Build an optimization model;

[0011] Step 4, iterative reconstruction process;

[0012] Step 5: Result evaluation and analysis.

[0013] Furthermore, in step 1, undersampling is performed using a variable density spiral trajectory, and the sampling rates are set to 20%, 30%, 40% and 50%, corresponding to acceleration factors of 5×, 3.3×, 2.5× and 2×, respectively; and the original k-space data is undersampled to obtain undersampled observation data y.

[0014] Furthermore, the step 2 comprises:

[0015] (1) Select sparse transform: Use contour wavelet transform as sparse representation method;

[0016] (2) Apply contour wavelet transform: Perform contour wavelet decomposition on the undersampled observation data y to obtain a sparse coefficient matrix θ, which is as follows:

[0017] θ=Ψ -1 y

[0018] Where Ψ represents the contour wavelet transform matrix, Ψ -1 is its inverse transform, used to convert the sparse coefficients back to the image domain;

[0019] Construct the encoding matrix A, combined with the two-dimensional discrete Fourier transform and undersampling mode, defined as:

[0020] A=F u Ψ

[0021] Among them, F u represents the two-dimensional discrete Fourier transform matrix combined with the undersampling mode, Ψ is the contour wavelet transform matrix; the encoding matrix A is used to map the sparse coefficients to the observed data.

[0022] Furthermore, the step 3 comprises:

[0023] (1) Optimization problem definition: Based on the compressed sensing theory, a convex optimization problem with the L1 norm as the core is constructed to achieve sparse reconstruction. The optimization objective function is defined as follows:

[0024]

[0025] Among them, λ is the regularization parameter, which is used to balance the data consistency term and the sparsity constraint;

[0026] (2) Adding a reweighting mechanism: To further improve the effect of sparse representation, a reweighting mechanism is introduced to dynamically adjust the regularization parameter μ k , balancing iteration speed and reconstruction accuracy.

[0027] Furthermore, the Lagrange multiplier method is used to constrain the optimization objective function and obtain the optimal solution:

[0028]

[0029] Furthermore, the step 4 includes:

[0030] (1) Initialization:

[0031] Set the initial sparse coefficient θ 0 =0;

[0032] Set the initial number of iterations k = 0 and the initial update parameter u 0 =0;

[0033] (2) Iterative update:

[0034] Update the regularization parameter μ k , the formula is as follows:

[0035]

[0036] Update consistency item G k (θ k ,θ k-1 ), the formula is as follows:

[0037]

[0038] Introducing factor μ k Dynamically adjust α and perform soft threshold shrinkage. The formula is as follows:

[0039]

[0040] Combined with the p-threshold function, the shrinkage function formula is as follows:

[0041]

[0042] Update number of iterations:

[0043] k=k+1

[0044] (3) Termination conditions:

[0045] The convergence condition of the iteration is determined by the relative error between the values ​​of two adjacent iterations, and the termination condition of the image iteration is given by the following formula:

[0046]

[0047] The termination threshold is set to ε, and the iteration is stopped if it is met;

[0048] (4) Output the reconstructed image:

[0049] The final sparse coefficient θ is transformed by contour wavelet inverse to obtain the reconstructed MR image u, the formula is as follows:

[0050] u=Ψ-1θ

[0051] Among them, -1 represents the contour wavelet inverse transform matrix, which is used to transform the sparse coefficients θ back to the image domain to obtain the reconstructed MR image u.

[0052] Further, the step 5 comprises:

[0053] (1) Reconstructed image display: Compare the reconstruction effects of different algorithms at different sampling rates, and display the reconstructed image and reconstruction error image;

[0054] (2) Performance index calculation: The reconstruction quality is quantitatively evaluated using three indicators: mean square error, peak signal-to-noise ratio, and structural similarity index.

[0055] The beneficial effect of the present invention is that the RISTA proposed by the present invention significantly improves the performance in fast MRI reconstruction through adaptive regularization weight adjustment. The RISTAC algorithm combined with contourlet transform further enhances the image detail restoration capability, confirms the effective improvement of traditional wavelet transform, and provides a useful reference for fast MRI imaging research and practical application. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] Figure 1 The flowchart of the method for fast magnetic resonance image reconstruction based on the reweighted shrinkage threshold algorithm is shown;

[0057] Figure 2 Reconstructed image of human neck at 30% sampling rate;

[0058] Figure 3 This is the error image of human neck reconstruction at 30% sampling rate;

[0059] Figure 4 To compare the effects of neck reconstruction with different evaluation indicators;

[0060] Figure 5 Reconstructed image of human knee at 30% sampling rate;

[0061] Figure 6 Reconstruction error image of human knee at 30% sampling rate;

[0062] Figure 7To compare the knee reconstruction effects with different evaluation indicators;

[0063] Figure 8 Reconstructed image of human foot at 30% sampling rate;

[0064] Fig. 9 This is the reconstruction error image of the human foot at a sampling rate of 30%;

[0065] Fig.10 To compare the effects of foot reconstruction using different evaluation indicators. DETAILED DESCRIPTION

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

[0067] This embodiment discloses a fast magnetic resonance image reconstruction method based on a reweighted shrinkage threshold algorithm, and describes the specific steps and experimental results of applying the reweighted iterative threshold shrinkage algorithm (RISTA) and its variant algorithm combined with contour wavelet transform (RISTAC) to accelerate the reconstruction of magnetic resonance imaging (MRI) under the compressed sensing framework. The specific implementation steps are as follows:

[0068] Step 1: Data collection and preprocessing

[0069] 1. Data Collection

[0070] Three MRI images depicting different anatomical regions of the human body (neck, knee, and foot) were selected for testing.

[0071] The original MRI was acquired using a 1.5T Philips Achieva scanner with an 8-channel receive coil.

[0072] The acquisition parameters are shown in Table 1.

[0073] Table 1 MRI acquisition parameters of test images

[0074]

[0075] 2. Undersampled k-space data:

[0076] A variable density spiral trajectory was used for undersampling, and the sampling rates were set to 20%, 30%, 40% and 50%, corresponding to acceleration factors of 5×, 3.3×, 2.5× and 2×, respectively.

[0077] The original k-space data is undersampled to obtain undersampled observation data y.

[0078] Step 2: Sparse transformation and coding matrix construction

[0079] 1. Choose sparse transformation

[0080] Contour wavelet transform is used as a sparse representation method to overcome the limitations of traditional wavelet transform in edge and contour detection. Contour wavelet transform can capture details and edges in images more effectively and reduce artifacts and detail loss.

[0081] 2. Apply contour wavelet transform:

[0082] Perform contour wavelet decomposition on the undersampled observation data y to obtain the sparse coefficient matrix θ, which is as follows:

[0083] θ=Ψ -1 y (1)

[0084] Where Ψ represents the contour wavelet transform matrix, Ψ -1 is its inverse transform, which is used to convert the sparse coefficients back to the image domain.

[0085] Construct the encoding matrix A, combined with the two-dimensional discrete Fourier transform and undersampling mode, defined as:

[0086] A=F u Ψ (2)

[0087] Among them, F u represents the two-dimensional discrete Fourier transform matrix combined with the undersampling mode, Ψ is the contour wavelet transform matrix. The encoding matrix A is used to map the sparse coefficients to the observed data.

[0088] Step 3: Build an optimization model

[0089] 1. Optimization problem definition:

[0090] According to the compressed sensing theory, a convex optimization problem with L1 norm as the core is constructed to achieve sparse reconstruction.

[0091] The optimization objective function is defined as follows:

[0092]

[0093] Here, λ is a regularization parameter used to balance the data consistency term and the sparsity constraint.

[0094] 2. Add a re-weighting mechanism

[0095] In order to further improve the effect of sparse representation, a reweighting mechanism is introduced to dynamically adjust the regularization parameter μ k , balancing iteration speed and reconstruction accuracy.

[0096] Step 4: Iterative Reconstruction Process

[0097] 1. Initialization:

[0098] Set the initial sparse coefficient θ 0 =0.

[0099] Set the initial number of iterations k = 0 and the initial update parameter u 0 =0.

[0100] 2. Iterative update:

[0101] Update the regularization parameter μ k , the formula is as follows:

[0102]

[0103] Update consistency item G k (θ k ,θ k-1 ), the formula is as follows:

[0104]

[0105] Introducing factor μ k Dynamically adjust α and perform soft threshold shrinkage. The formula is as follows:

[0106]

[0107] Combined with the p-threshold function, the shrinkage function formula is as follows:

[0108]

[0109] Update number of iterations:

[0110] k=k+1 (8)

[0111] 3. Termination conditions:

[0112] The convergence condition of the iteration is determined by the relative error between the values ​​of two adjacent iterations, and the termination condition of the image iteration is given by the following formula:

[0113]

[0114] The termination threshold is set to ε, and the iteration is stopped if it is met.

[0115] 4. Output the reconstructed image:

[0116] The final sparse coefficient θ is transformed by contour wavelet inverse to obtain the reconstructed MR image u. The formula is as follows:

[0117] u=Ψ-1θ (10)

[0118] Among them, -1 represents the contour wavelet inverse transform matrix, which is used to transform the sparse coefficients θ back to the image domain to obtain the reconstructed MR image u.

[0119] Step 5: Result evaluation and analysis

[0120] 1. Reconstructed image display:

[0121] The reconstruction effects of five different algorithms (ISTA, FISTA, GTIA, RISTA and RISTAC) at different sampling rates (20%, 30%, 40%, 50%) are compared, and the reconstructed images and reconstruction error images are displayed.

[0122] 2. Performance index calculation:

[0123] The reconstruction quality is quantitatively evaluated using three indicators: mean square error (MSE), peak signal-to-noise ratio (PSNR) and structural similarity index (SSIM).

[0124] In summary, the image reconstruction principle of the present invention mainly includes the following contents:

[0125] Content 1: Magnetic Resonance Imaging Model under Compressed Sensing

[0126] In order to apply the principle of compressed sensing to accelerate magnetic resonance imaging, the preliminary steps include using a subset of k-space data obtained by undersampling as a known constraint. Subsequently, using the prior knowledge of the sparsity of magnetic resonance images in a specific transform domain, a convex optimization problem centered on the L1 norm is constructed based on the compressed sensing theory. Under the premise of incomplete measurement, the reconstruction of magnetic resonance (MR) images involves solving the following optimization problem:

[0127]

[0128] Where m∈C N Represents the MR image to be reconstructed, y∈C M is the undersampled MR image data in k-space, is the observation matrix, expressed as:

[0129] F u =kF

[0130] Where k represents the diagonal matrix of the undersampling mask in k-space and represents the two-dimensional discrete Fourier transform.

[0131] Equation F u It can be transformed into a minimum norm optimization problem:

[0132]

[0133] The MR image m can be sparsely represented as:

[0134] m=Ψθ

[0135] Where Ψ=[ψ 1 ,ψ 2 ,K,ψ N ]∈RN×N is the sparse transformation matrix, and θ is the projection coefficient in the sparse transformation domain. Therefore, equation F u can be rewritten as:

[0136]

[0137] Among them, F u Ψ T represents the encoding matrix.

[0138] For an uncertain system, taking into account the noise in the magnetic resonance scanning process, the above equation can be expressed as follows:

[0139]

[0140] MR image reconstruction involves both sparsity and consistency of k-space data, which is solved by the following Lagrangian constraint equation:

[0141]

[0142] Among them, λ is the regularization parameter, which is used to balance the weights of the two items. The above equation is constrained by the Lagrange multiplier method to obtain the optimal solution:

[0143]

[0144] Content 2: Reweighted Iterative Shrinkage Threshold Algorithm (RISTA)

[0145] The main concept of RISTA is to track the iteration process and automatically determine the iteration stage. By analyzing the results of the first two iterations, RISTA automatically adjusts the iteration weight λ to balance the iteration speed and accuracy. The specific process is as follows:

[0146] 1. Perform soft threshold shrinkage

[0147]

[0148] 2. Introducing factor μ k Dynamically adjusting α, the above equation can be rewritten as:

[0149]

[0150] 3. Regularization parameter μ k The update function is as follows:

[0151]

[0152] During the iteration process, θ k-1 and θ k-2 represents the reconstruction results of the first two iterations. As the iterations proceed, μ kGradually tends to 1. In order to optimize the trade-off between convergence speed and reconstruction accuracy, the regularization parameter λ is dynamically adjusted. In the early stages of the iteration, more emphasis is placed on the fidelity term to promote rapid convergence; while in the later iterations, the regularization term dominates to ensure the accuracy of reconstruction. The parameter ξ is used to control the rate of change of the regularization parameter, thereby fine-tuning the iteration process.

[0153] 4. Combined with the p-threshold function, the shrinkage function becomes:

[0154]

[0155] where p is a parameter used to adjust the speed at which different elements in the matrix θ shrink.

[0156] 5. Iteration termination conditions:

[0157] The convergence condition of the iteration is determined by the relative error between the values ​​of two adjacent iterations, and the termination condition of the image iteration is given by the following formula:

[0158]

[0159] The termination threshold is set to ε, and the iteration is stopped if it is met.

[0160] 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.

[0161] The parts not involved in the present invention are the same as the prior art or can be implemented by using the prior art.

Claims

1. A fast magnetic resonance image reconstruction method based on a reweighted shrinkage threshold algorithm, characterized in that: The steps include: Step 1: Data collection and preprocessing; Step 2: Sparse transformation and coding matrix construction; Step 3: Build an optimization model; Step 4, iterative reconstruction process; Step 5: Result evaluation and analysis.

2. A fast magnetic resonance image reconstruction method based on a reweighted shrinkage threshold algorithm according to claim 1, characterized in that: In step 1, undersampling is performed using a variable density spiral trajectory, and the sampling rates are set to 20%, 30%, 40% and 50%, corresponding to acceleration factors of 5×, 3.3×, 2.5× and 2×, respectively; and the original k-space data is undersampled to obtain undersampled observation data y.

3. The method for rapid magnetic resonance image reconstruction based on a reweighted shrinkage threshold algorithm according to claim 1, characterized in that: The step 2 comprises: (1) Select sparse transform: Use contour wavelet transform as sparse representation method; (2) Apply contour wavelet transform: Perform contour wavelet decomposition on the undersampled observation data y to obtain a sparse coefficient matrix θ, which is as follows: θ=Ψ -1 y Where Ψ represents the contour wavelet transform matrix, Ψ -1 is its inverse transform, used to convert the sparse coefficients back to the image domain; Construct the encoding matrix A, combined with the two-dimensional discrete Fourier transform and undersampling mode, defined as: A=F u P Among them, F u represents the two-dimensional discrete Fourier transform matrix combined with the undersampling mode, Ψ is the contour wavelet transform matrix; the encoding matrix A is used to map the sparse coefficients to the observed data.

4. The method for rapid magnetic resonance image reconstruction based on a reweighted shrinkage threshold algorithm according to claim 1, characterized in that: The step 3 comprises: (1) Optimization problem definition: Based on the compressed sensing theory, a convex optimization problem with the L1 norm as the core is constructed to achieve sparse reconstruction. The optimization objective function is defined as follows: Among them, λ is the regularization parameter, which is used to balance the data consistency term and the sparsity constraint; (2) Adding a reweighting mechanism: To further improve the effect of sparse representation, a reweighting mechanism is introduced to dynamically adjust the regularization parameter μ k , balancing iteration speed and reconstruction accuracy.

5. The method for rapid magnetic resonance image reconstruction based on a reweighted shrinkage threshold algorithm according to claim 4, characterized in that: The Lagrange multiplier method is used to constrain the optimization objective function and obtain the optimal solution:

6. The method for rapid magnetic resonance image reconstruction based on a reweighted shrinkage threshold algorithm according to claim 5, characterized in that: The step 4 comprises: (1) Initialization: Set the initial sparse coefficient θ 0 =0; Set the initial number of iterations k = 0 and the initial update parameter u 0 =0; (2) Iterative update: Update the regularization parameter μ k , the formula is as follows: Update consistency item G k (θ k ,θ k-1 ), the formula is as follows: Introducing factor μ k Dynamically adjust α and perform soft threshold shrinkage. The formula is as follows: Combined with the p-threshold function, the shrinkage function formula is as follows: Update number of iterations: k=k+1 (3) Termination conditions: The convergence condition of the iteration is determined by the relative error between the values ​​of two adjacent iterations, and the termination condition of the image iteration is given by the following formula: The termination threshold is set to ε, and the iteration is stopped if it is met; (4) Output the reconstructed image: The final sparse coefficient θ is transformed by contour wavelet inverse to obtain the reconstructed MR image u, the formula is as follows: u=Ψ-1θ Among them, -1 represents the contour wavelet inverse transform matrix, which is used to transform the sparse coefficients θ back to the image domain to obtain the reconstructed MR image u.

7. The method for rapid magnetic resonance image reconstruction based on a reweighted shrinkage threshold algorithm according to claim 1, characterized in that: The step 5 comprises: (1) Reconstructed image display: Compare the reconstruction effects of different algorithms at different sampling rates, and display the reconstructed image and reconstruction error image; (2) Performance index calculation: The reconstruction quality is quantitatively evaluated using three indicators: mean square error, peak signal-to-noise ratio, and structural similarity index.

Images