A parallel magnetic resonance imaging fast reconstruction method based on transform learning and structured low-rank model
By combining transform learning and structured low-rank models, introducing JTL and ADMM and utilizing GPU acceleration, the problems of large computational cost and insufficient quality of parallel magnetic resonance imaging reconstruction models are solved, and fast, high-quality image reconstruction is achieved.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- KUNMING UNIV OF SCI & TECH
- Filing Date
- 2022-11-24
- Publication Date
- 2026-04-10
AI Technical Summary
Existing parallel magnetic resonance imaging reconstruction models have a huge computational load, and the reconstruction speed and quality need to be improved. In particular, the SAKE model relies on low-rank matrix recovery, which leads to excessive computation and blurry artifacts in the images.
By combining transform learning and structured low-rank models, we introduce Joint Sparse Transform Learning (JTL) and Progressive Optimized Gradient Method (ADMM), and utilize GPU acceleration to propose the JTLSAKE method, which iteratively optimizes the image reconstruction process.
It significantly reduces the reconstruction time of magnetic resonance imaging while maintaining good reconstruction quality, image clarity and detail preservation, and the reconstruction speed is 85 times faster than JTL-PLORAKS.
Smart Images
Figure CN115877296B_ABST
Abstract
Description
TECHNICAL FIELD
[0001] The application relates to a parallel magnetic resonance imaging fast reconstruction method based on transform learning and a structured low-rank model, and belongs to the technical field of magnetic resonance imaging. BACKGROUND
[0002] Magnetic resonance (MR) imaging is an image examination technology developed in the 1980s. Because it does not damage the human body and has good soft tissue resolution, it has become one of the important imaging methods for clinical diagnosis. However, the MR imaging scanning time is relatively long, which greatly affects the throughput and imaging quality of the magnetic resonance imaging instrument. Therefore, how to shorten the MR imaging scanning time has become a research hotspot, and the compressed sensing (CS) theory and the parallel imaging (PI) technology provide researchers with some new ideas.
[0003] The compressed sensing theory breaks through the limitation of the Shannon-Nyquist sampling theorem. The theory points out that as long as a signal has a sparse representation in a certain transform, the original signal can be recovered from a small amount of undersampled data through an optimization method. Since MR images have sparse representations in wavelet transforms, finite difference transforms and the like, the application of the compressed sensing theory can greatly improve the MR imaging scanning speed. In order to eliminate the artifacts caused by undersampling, researchers introduce the L1 regularization term of the total variation (TV) and the wavelet transform coefficient to improve the reconstruction quality of the image. However, the fixed sparse transform cannot adapt to the signal. Therefore, some adaptive sparse transform methods have been proposed, such as the methods based on dictionary learning (DL) and transform learning (TL).
[0004] The parallel imaging technology acquires data simultaneously through multiple receiving coils with different sensitivities. The most common reconstruction model of parallel imaging is sensitivity encoding (SENSE), but it requires accurate coil sensitivity information, which limits the application of the model. Subsequently, Lustig et al. proposed the iterative self-consistent parallel imaging reconstruction (SPIRiT) model in 2010 and introduced the L1 norm regularization of the wavelet coefficient, which achieved good reconstruction effect, but still needs to scan additional auto-calibration signals (ACS) for data calibration, which increases the scanning time of MR imaging.
[0005] Shin et al. proposed a structured low-rank parallel imaging reconstruction framework SAKE (Simultaneous Auto-calibrating and K-space Estimation), which can effectively reconstruct magnetic resonance images when ACS is limited, and is more flexible than SENSE and SPIRiT models. Subsequently, Haldar et al. proposed a PLORAKS parallel imaging reconstruction model (low-rank modeling of local k-space neighborhoods with parallel imaging), which can use calibration and non-calibration sampling modes when the undersampling rate is high. However, the above two models both rely on low-rank matrix recovery, resulting in huge computational load. Duan et al. combined the PLORAKS model with TL and proposed TL-PLORAKS and JTL-PLORAKS (joint denoising version), which can obtain better reconstruction quality and edge details, but are too complex and have huge computational load, which is not conducive to practical application. Researchers have proposed some methods implemented using a graphics processing unit (GPU), which can significantly shorten the reconstruction time without reducing the image reconstruction quality.
[0006] In summary, introducing a regularization term in the reconstruction model can effectively improve the reconstruction quality, and using a GPU for parallel optimization can greatly speed up the reconstruction while ensuring the image reconstruction quality. However, the SAKE model relies on low-rank matrix recovery, has huge computational load, and the reconstructed image has obvious blur artifacts, and the reconstruction speed and quality still need to be improved. SUMMARY
[0007] The purpose of the present application is to overcome the shortcomings of the prior art and provide a parallel magnetic resonance imaging fast reconstruction method based on transform learning and structured low-rank model, which can further improve the reconstruction speed of magnetic resonance images.
[0008] The present application is realized by the following technical solutions: a parallel magnetic resonance imaging fast reconstruction method based on transform learning and structured low-rank model. It includes the following steps:
[0009] 1. A parallel magnetic resonance imaging fast reconstruction method based on transform learning and structured low-rank model, comprising the following steps:
[0010] S0: initialization, f 0 = D H y, Z 0 = 0, t 0 =0, k=0;
[0011] Among them, the superscript " 0 " represents the initial value before iteration, This represents the column-vectorized multi-coil k-space data to be reconstructed, where l = 1...L represents the coil index variable, L represents the number of receiving coils used, and the superscript " H "" represents the conjugate transpose operation of a vector or matrix, and the spatial data of the l-th coil k of f is represented as For example: f1 and f L Let f represent the first and Lth coil k-space data of the multi-coil k-space data to be reconstructed, respectively, in column vectorized form; 0 This represents the initial value of f. This represents the multi-coil k-space undersampling operator. I represents the k-space undersampling matrix of a single coil. L Let L×L be the identity matrix. Indicates the Kronecker product. This represents undersampled multi-coil k-space data, N = N h ×N v N represents the number of pixels in the single-coil image to be reconstructed. h N v These represent the number of pixels in the horizontal and vertical directions of the magnetic resonance image, respectively, and M represents the actual number of sampling points for a single coil k-space data. F represents the coil-by-coil Fourier transform operator. h F v N h N v Point Fourier transform matrix, F -1 Indicates the inverse Fourier transform, -1 " indicates the inverse operation of a matrix; Extracting a linear operator for the j-th matrix indicates extracting from each coil image. Image blocks are recombined into The matrix; P1(·) denotes the linear operator extracted from the first matrix, This indicates that the linear operator is extracted from the N2-th matrix. This means combining the extracted image patches into a matrix X, where X... 0 N represents the initial value of X, and N2 represents the number of individual coil images. The number of image patches; express Adaptive sparse transformation matrix of image patch express Point discrete cosine transform matrix, W 0 This represents the initial value of W; Indicates the relationship with WP1(F) -1 f) is the corresponding auxiliary variable. Indicates and Corresponding auxiliary variables, Indicates to WP j (F -1 f) corresponds to the auxiliary variable, j = 1......N2 represents the index variable, B k This represents the value of B obtained in the k-th iteration. B represents the k-th iteration. k The j-th block; Indicates the relationship with WP1(F) -1 f) is the corresponding auxiliary variable. Indicates and Corresponding auxiliary variables, Indicates to WP j (F -1 f) is the corresponding auxiliary variable. This represents the result obtained in the k-th iteration. Indicates the k-th iteration The j-th block, B 0 and They represent B and The initial value; This represents the Lagrange multiplier corresponding to the auxiliary variable B. Indicate u B The j-th block, Indicate u B The first block, Indicate u B The N2nd block, Indicate u B The initial value; Indicates with R C (f) corresponds to the auxiliary variable, R C (f)=[R1(f)......R N1 (f)], operator R C (·): This indicates extracting one from the multi-coil k-space. The matrix, operator R i (·): This indicates extraction from space for each coil k. The blocks are converted into column vectors according to the coil column vectors. The i-th operator, i = 1...N1, where N1 represents the space containing a single coil k. The number of blocks; Indicates with R C (f) corresponds to the auxiliary variable Z 0 and denote the initial values of Z and denote the Lagrange multiplier corresponding to the auxiliary variable Z, denote the initial values of u Z denote the joint hard thresholding function, z denotes the input vector or matrix, θ denotes the threshold, and α > 0, μ2 > 0 are parameters; denote the time acceleration factor, t 0 denote the initial value of t; k denotes the iteration number;
[0012] S1: perform singular value decomposition on , i.e. Γ is a singular value matrix, and Ψ and Φ denote left and right singular matrices, respectively, so that the adaptive sparse transformation matrix W k+1 of the image block in the k+1th iteration can be obtained, and the calculation formula is as follows:
[0013] W k+1 = ΦΨ H
[0014] where X k denotes the matrix composed of the extracted image blocks in the kth iteration, B k is the auxiliary variable in the kth iteration, is the Lagrange multiplier corresponding to B k in the kth iteration;
[0015] S2: calculate the auxiliary variable B in the k+1th iteration, and the calculation formula is as follows:
[0016]
[0017] where f k is the multi-coil k-space data reconstructed in the kth iteration, and are the values of and in the kth iteration.
[0018] S3: update the time acceleration factor t k+1 in the k+1th iteration, and the calculation formula is as follows:
[0019]
[0020] where t k denotes the time acceleration factor in the kth iteration;
[0021] S4: update the auxiliary variable B k+1 in the k+1th iteration, and the calculation formula is as follows:
[0022]
[0023] in, B k The auxiliary variable represents the k-th iteration;
[0024] S5: Yes Perform singular value decomposition, i.e. U and V are the left and right singular matrices, respectively, and Σ = diag(ω), where ω contains... All singular values are obtained by performing low-rank estimation on U, Σ, and V using the low-rank estimation parameter r and then truncating them. as well as Then the auxiliary variable in the (k+1)th iteration for:
[0025]
[0026] in, This represents the auxiliary variable for the k-th iteration.
[0027] S6: Update the low-rank matrix auxiliary variable Z in the (k+1)th iteration. k+1 The calculation formula is as follows:
[0028]
[0029] Among them, Z k The auxiliary variable represents the k-th iteration;
[0030] S7: Calculate the k-space data f of the multi-coil image to be reconstructed in the (k+1)th iteration. k+1 The calculation formula is as follows:
[0031]
[0032] in, To be with R C The adjoint operator corresponding to (·) This forms an LN×LN diagonal matrix. To be with P j The adjoint operator corresponding to (·) Form an LN×LN diagonal matrix, where I is the LN×LN identity matrix, and μ1>0 is a parameter;
[0033] S8: Calculate the matrix X formed by combining the extracted image patches in the (k+1)th iteration. k+1 The calculation formula is as follows:
[0034]
[0035] S9: Calculate B of the k+1th iteration k+1 The corresponding Lagrange multiplier The calculation formula is as follows:
[0036]
[0037] S10: Calculate Z of the k+1th iteration k+1 The corresponding Lagrange multiplier The calculation formula is as follows:
[0038]
[0039] S11: Calculate f k+1 Inverse Fourier transform, and then square root of sum-of-squares (SOS) operation, to obtain a single coil amplitude reconstruction image, the calculation formula is as follows:
[0040]
[0041] Wherein, Indicates the single coil amplitude reconstruction image of the k+1th iteration;
[0042] S12: Calculate x k+1 And the relative error (Relative Error, RE) between x k The calculation formula is as follows:
[0043]
[0044] S13: Determine whether the iteration stopping criterion is met, if the maximum iteration number K (i.e. k=K) is reached or RE <tol is met, go to step S14; otherwise, let k=k+1 and return to S1;
[0045] Wherein, tol is the set tolerance;
[0046] S14: Output the final reconstructed single coil amplitude image x k+1 .
[0047] The beneficial effects of the present application are: in the present application, the joint sparse transform learning (Joint Transform Learning, JTL) is introduced into the SAKE model, ADMM is used for solving, and the optimization gradient method is introduced to speed up the convergence speed, finally GPU is used for acceleration, so that a fast parallel magnetic resonance imaging reconstruction method JTLSAKE is obtained. Theoretical analysis and imaging experiments show that, compared with JTL-PLORAKS, the parallel imaging reconstruction method JTLSAKE proposed in the present application can greatly reduce the reconstruction time of magnetic resonance imaging, and the reconstruction quality is equivalent. BRIEF DESCRIPTION OF DRAWINGS
[0048] Figure 1 Flow chart of the method of the present application;
[0049] Figure 2 Image of the fully sampled dataset1 of the human brain data of a subject using a coil with eight channels;
[0050] Figure 3 Two-dimensional Poisson disc undersampling mask with 4x sampling rate and 24x24 ACS line;
[0051] Figure 4 Figure 5 Images reconstructed from 4x accelerated and 24x24 center-corrected two-dimensional Poisson disc undersampling data using JTLSAKE and JTL-PLORAKS, respectively, under dataset1 are shown.
[0052] Figure 6 Figure 7 Error images corresponding to Figure 4 Figure 5 are shown in DETAILED DESCRIPTION
[0053] The technical solutions of the present application are described in further detail below in conjunction with the drawings.
[0054] Example 1: The present application is a highly efficient reconstruction method based on the SAKE framework.
[0055] Let denote the column vectorized multi-coil k-space data to be reconstructed, l = 1...L denote the coil index variable, L denote the number of receive coils used, “ H ” denote the transpose operation of a matrix, the l-th coil k-space data of f be denoted as For example: f1and f L denote the 1st coil k-space data and the Lth coil k-space data of the column vectorized multi-coil k-space data f to be reconstructed, respectively; f 0 denote the initial value of f, denote the multi-coil k-space undersampling operator, denote the single-coil k-space undersampling matrix, I L be the L x L identity matrix,
[0056] denote the undersampled multi-coil k-space data, N = N h v denote the number of pixels of the single-coil image to be reconstructed, N h v These represent the number of pixels in the horizontal and vertical directions of the magnetic resonance image, respectively, and M represents the actual number of sampling points for a single coil k-space data. Operator R C (·): This indicates extracting one from the multi-coil k-space. The matrix, operator R i (·): This indicates extraction from each coil k space. The blocks are converted into column vectors according to the coil column vectors. The i-th operator, i = 1...N1, where N1 represents the space containing a single coil k. Given the number of blocks and r as the low-rank estimation parameter, the SAKE reconstruction model can be obtained as follows:
[0057]
[0058] Where ||·||2 represents the L2 norm of the vector.
[0059] To improve the reconstruction quality of parallel MRI, this invention combines the JTL regularization term with the SAKE model to propose the JTLSAKE model, thus yielding the following optimization problem:
[0060]
[0061] Where α is the regularization parameter, and N2 represents the number of elements in a single coil image. The number of image patches, express The adaptive sparse transformation matrix of the image patch, P j (·): Extracting a linear operator for the j-th matrix indicates extracting from each coil image. Image blocks are recombined into The matrix, For the Fourier transform operator, F h F v N h N v Point Fourier transform matrix, I L Let L×L be the identity matrix. F represents the Kronecker product. -1 Indicates the inverse Fourier transform, -1 " represents the inverse operation of a matrix, ||·|| 0,2 The joint l0 norm is defined as follows: w represents the input vector or matrix, W H Denotes the conjugate transpose of W, H " represents the conjugate transpose operation of a matrix.
[0062] Introducing auxiliary variables and the corresponding Lagrange multipliers and Using ADMM technology, problem (2) can be transformed into solving subproblems of the following form:
[0063]
[0064]
[0065]
[0066]
[0067]
[0068]
[0069] In formulas (3)-(8), the superscript of the variable " k+1 "and" k "" represents the variables obtained after the (k+1)th and kth iterations, respectively. In formula (3), express The adaptive sparse transformation matrix W of the image patch k+1 This represents the adaptive sparse transformation matrix of the image patch in the (k+1)th iteration, where μ2 > 0 is a parameter. Indicates the relationship with WP1(F) -1 f) is the corresponding auxiliary variable. Indicates and Corresponding auxiliary variables, Indicates to WP j (F -1 f) corresponds to the auxiliary variable, j = 1......N2 represents the index variable, B k This represents the value of B obtained in the k-th iteration. B represents the k-th iteration. k The j-th block, This represents the Lagrange multiplier corresponding to the auxiliary variable B. Indicate u B The j-th block, Indicate u B The first block, Indicate u B The N2nd block, Indicate u B initial value, Indicates the k-th iteration The j-th block, f kk-space data of the multi-coil image to be reconstructed at the kth iteration; in equation (4), Bk+1denotes the multi-coil image at the k+1th iteration, k+1 the jth block of Bk; in equation (5), Zk k+1 denotes the k+1th iteration of the Lagrange multiplier corresponding to R C (f) the corresponding auxiliary variable, denotes the Lagrange multiplier corresponding to Z, denotes the 1st column of the matrix u Z denotes the N1th column of the matrix u Z , i = 1... N1, denotes the kth iteration of the matrix Z k the corresponding dual variable, || · || F denotes the Frobenius norm; in equation (6), f k+1 k-space data of the multi-coil image to be reconstructed at the k+1th iteration; in equation (7), denotes the k+1th iteration of the matrix B the jth block of Bk+1; in equation (8), denotes the k+1th iteration of the matrix Z k+1 the corresponding Lagrange multiplier.
[0070] The above sub-problem is solved next.
[0071] Solving the transformation learning W: since and Problem (3) can be converted into the following form:
[0072]
[0073] where Xk k denotes the matrix composed of the extracted image blocks at the kth iteration, Bk k denotes the auxiliary variable at the kth iteration, denotes the matrix B k the corresponding Lagrange multiplier, || · || F denotes the Frobenius norm.
[0074] Then, problem (9) can be solved using the following method:
[0075]
[0076] W k+1 = ΦΨ H (11)
[0077] where SVD(·) denotes singular value decomposition of a matrix, Γ is a singular value matrix, and Ψ and Φ are left and right singular matrices, respectively.
[0078] Solving sparse coding B j Equation (4) can be simplified as:
[0079]
[0080] Problem (12) can be solved by joint hard thresholding method, i.e.,
[0081]
[0082] where H J (z, θ) is defined as follows:
[0083]
[0084] where z is an input vector or matrix, and θ is a threshold.
[0085] According to and Equation (13) can be further transformed as:
[0086]
[0087] where B k+1 denotes the auxiliary variable of the k+1th iteration.
[0088] Using the OGM fast iterative method to accelerate the sparse coding B, we have:
[0089]
[0090]
[0091]
[0092] where, and denote the auxiliary variables of the k+1th and kth iterations, respectively, t k+1 and t k denote the time acceleration factors of the k+1th and kth iterations, respectively.
[0093] Solving low-rank matrix Z: assuming is the singular value decomposition of , then the solution of problem (5) can be represented as:
[0094]
[0095] where Zk+1 denotes the auxiliary variable of the k+1th iteration, ∑ is the singular value matrix, U and V denote the left and right singular matrices respectively, Σ = diag(ω), ω contains all singular values of Z, and are obtained by using the parameter r to perform low-rank estimation and truncation on U, Σ and V respectively.
[0096] The OGM fast iteration method is used to accelerate the low-rank matrix Z, and then:
[0097]
[0098]
[0099] where, and denote the auxiliary variable of the k+1th iteration, Z k denotes the auxiliary variable of the kth iteration.
[0100] Updating the image f: let the derivative of the objective function of problem (6) be 0, then:
[0101]
[0102] where μ1>0 is a parameter, is the adjoint operator of R C (·), forms a LN×LN identity matrix, and “ * ” denotes the adjoint operator of the matrix; is the adjoint operator of P j (·), forms a LN×LN diagonal matrix, and I is a LN×LN identity matrix. Since W H W = I, then:
[0103]
[0104] Therefore, problem (22) can be further simplified as:
[0105]
[0106] Since all the matrices on the left side of equation (24) are diagonal matrices, their inverses are easy to obtain, and thus the analytical solution of f is:
[0107]
[0108] Substituting f k+1The inverse Fourier transform is performed, and then the SOS (Sum Of Squares) operation is performed to obtain a single amplitude image, and the operation can be expressed as:
[0109]
[0110] Wherein, x k+1 is the finally obtained reconstructed image.
[0111] Solving the Lagrange multiplier u B corresponding to B And The formula (7) can be further expressed as:
[0112]
[0113] Wherein, B k+1 corresponding to the k+1th iteration
[0114] In summary, the application obtains a new parallel magnetic resonance imaging reconstruction method JTLSAKE. Finally, using the relative error (Relative Error, RE) less than the tolerance tol or reaching the maximum iteration number K as the stopping condition of the method, and the definition of RE is as follows:
[0115]
[0116] The specific process is shown in Figure 1 , and the steps are as follows:
[0117] S0: initialization, f 0 =D H y,
[0118]
[0119]
[0120] Z 0 =0, t 0 =0, k=0;
[0121] Wherein, the superscript "initial value before iteration 0 " indicates the initial value of f 0 , f B indicates the initial value of f, and respectively indicate the initial value of and , and indicates the initial value of u denotes the initial value of u Z denotes the initial value of X 0 denotes the initial value of X 0 denotes the initial value of W 0 denotes the initial value of B 0 denotes the initial value of t, k denotes the iteration number
[0122] S1: perform singular value decomposition on , i.e. Γ is a singular value matrix, and Ψ and Φ represent left and right singular matrices respectively, then the adaptive sparse transform matrix W k+1 of the image block in the k+1th iteration can be obtained, as shown in equation (11);
[0123] S2: calculate the auxiliary variable in the k+1th iteration, as shown in equation (16);
[0124] S3: update the time acceleration factor t k+1 in the k+1th iteration, as shown in equation (17);
[0125] S4: update the auxiliary variable B k+1 in the k+1th iteration, as shown in equation (18);
[0126] S5: perform singular value decomposition on , i.e. and use the low-rank estimation parameter r to perform low-rank estimation on U, Σ and V and truncate to obtain and then the auxiliary variable in the k+1th iteration can be obtained, as shown in equation (20);
[0127] S6: update the auxiliary variable Z k+1 of the to-be-reconstructed multi-coil image in the k+1th iteration, as shown in equation (21);
[0128] S7: calculate the k-space data f k+1 of the to-be-reconstructed multi-coil image in the k+1th iteration, as shown in equation (25);
[0129] S8: calculate the matrix X k+1 composed of the extracted image blocks in the k+1th iteration, as shown below:
[0130]
[0131] S9: calculate B k+1 corresponding to the Lagrange multiplier in the k+1th iteration , as shown in equation (27);
[0132] S10: Calculate Z of the k+1th iteration k+1 The corresponding Lagrange multiplier The formula is as shown in equation (8);
[0133] S11: Calculate f k+1 Perform inverse Fourier transform, and then perform square root of sum-of-squares (SOS) operation, to obtain a single-coil amplitude reconstruction image, the formula is as shown in equation (26);
[0134] S12: Calculate x k+1 and the relative error (RE) between x k , the formula is as shown in equation (28);
[0135] S13: Determine whether the iteration stopping criterion is met, if the maximum iteration number K (i.e., k=K) is reached or the RE < tol is met, go to step S14; otherwise, let k=k+1 and return to step S1;
[0136] S14: Output the final reconstructed single-coil amplitude image x k+1 .
[0137] Experimental results:
[0138] In the following experiments, in order to verify the reconstruction performance of the proposed method, the present application compares it with JTL-PLORAKS. All methods are implemented using MATLAB (MathWorks, Natick, MA). All experiments are performed on a desktop computer configured with an Intel(R) Core(TM) i9-12900K@3.20GHz processor, 64GB of memory, and an NVIDIA GeForce RTX 3090 GPU.
[0139] In order to compare the reconstruction performance of each method, the present application uses one human brain image for simulation experiments, named dataset1 (as shown in Figure 2 ). The present application uses a two-dimensional Poisson disc undersampling mask with an acceleration factor of 4 and 24x24 ACS lines to artificially undersample the fully sampled dataset to generate the test dataset, as shown in Figure 3 .
[0140] In the experiment, the signal noise ratio (SNR) is used to evaluate the reconstruction performance of the method. The larger the value of SNR, the higher the quality of image reconstruction. The indicators are calculated in the region of interest. All experiments are obtained by adjusting the parameters of the method to obtain the optimal SNR value. The formula for calculating SNR is as follows:
[0141]
[0142] wherein Var represents the mean square error between the reference image x and the reconstructed image , and MSE is the variance of the reference image x.
[0143] The present application compares the visual presentation of the reconstructed images using each method on the selected data set with an acceleration factor of 4. Figure 4 and Figure 5 The present application shows the reconstructed images using JTLSAKE and JTL-PLORAKS on the data set dataset1. Figure 4 and Figure 5 It can be seen that both JTLSAKE and JTL-PLORAKS have good visual effects and similar reconstruction quality, and the reconstructed images effectively remove artifacts and preserve image details, and have good reconstruction quality.
[0144] In order to more intuitively compare the reconstruction effects of several methods, the present application shows the error images of JTLSAKE and JTL-PLORAKS on the data set dataset1. Figure 6 and Figure 7 It can be seen that the error images of JTL-PLORAKS and JTLSAKE have no significant difference, and the errors of the reconstructed images and the original images are small, indicating that the reconstruction quality of the two methods is similar. Figure 6 and Figure 7 It can be seen that the error images of JTL-PLORAKS and JTLSAKE have no significant difference, and the errors of the reconstructed images and the original images are small, indicating that the reconstruction quality of the two methods is similar.
[0145] In addition, the present application also analyzes the numerical analysis of the SNR of the reconstructed images of the two methods. The SNR of the reconstructed image (such as Figure 4 ) of JTLSAKE is 27.81 dB, and the SNR of the reconstructed image (such as Figure 5 ) of JTL-PLORAKS is 27.74 dB. It can be seen that under the evaluation index, the SNR values of JTLSAKE and JTL-PLORAKS proposed by the present application are not much different, further verifying the results of the visual comparison analysis.
[0146] Finally, the method JTLSAKE and JTL-PLORAKS proposed in this paper have better reconstruction quality, but the reconstruction speed is very slow, so the application calls GPU to speed up the method JTLSAKE to obtain the accelerated version JTLSAKE-GPU of JTLSAKE, which greatly shortens the running time of the method. For dataset1, the reconstruction running time of JTLSAKE-GPU and JTL-PLORAKS is 8.2 seconds and 698.2 seconds respectively, and the reconstruction speed of JTLSAKE-GPU is 85 times that of JTL-PLORAKS, which is more practical.
[0147] The application combines the JTL regularization term with the SAKE model and proposes a new parallel magnetic resonance imaging method, which is called parallel magnetic resonance imaging fast reconstruction method based on transform learning and structured low-rank model (JTLSAKE). After a plurality of experiments, the experimental results show that the method JTLSAKE proposed in the application greatly reduces the reconstruction time of JTL-PLORAKS and has close reconstruction quality.
[0148] The specific embodiments of the application are described in detail above with reference to the drawings, but the application is not limited to the above-mentioned embodiments, and various changes can be made within the knowledge of those skilled in the art without departing from the purpose of the application.
Claims
1. A parallel magnetic resonance imaging fast reconstruction method based on transform learning and structured low-rank model, comprising the following steps: S0: initialization, f 0 = D H y, Z 0 = 0, t 0 = 0, k = 0; Among them, the superscript " 0 " represents the initial value before iteration, This represents the column-vectorized k-space data of the multi-coil receiver to be reconstructed, where l = 1...L represents the coil index variable, L represents the number of receiving coils used, and the superscript " H "" represents the conjugate transpose operation of a vector or matrix, and the spatial data of the l-th coil k of f is represented as f1 and f L Let f represent the first and Lth coil k-space data of the multi-coil k-space data to be reconstructed, respectively, in column vectorized form; 0 This represents the initial value of f. This represents the multi-coil k-space undersampling operator. I represents the k-space undersampling matrix of a single coil. L Let L×L be the identity matrix. Indicates the Kronecker product. This represents undersampled multi-coil k-space data, N = N h ×N v N represents the number of pixels in the single-coil image to be reconstructed. h N v These represent the number of pixels in the horizontal and vertical directions of the magnetic resonance image, respectively, and M represents the actual number of sampling points for a single coil k-space data. F represents the coil-by-coil Fourier transform operator. h F v N h N v Point Fourier transform matrix, F -1 Indicates the inverse Fourier transform, -1 " indicates the inverse operation of a matrix; P j (·): Extracting a linear operator for the j-th matrix indicates extracting from each coil image. Image blocks are recombined into The matrix; P1(·) denotes the linear operator extracted from the first matrix, This indicates that the linear operator is extracted from the N2-th matrix. This means combining the extracted image patches into a matrix X, where X... 0 N represents the initial value of X, and N2 represents the number of individual coil images. The number of image patches; express Adaptive sparse transformation matrix of image patch express Point discrete cosine transform matrix, W 0 This represents the initial value of W; Indicates the relationship with WP1(F) -1 f) is the corresponding auxiliary variable. Indicates and Corresponding auxiliary variables, Indicates to WP j (F -1 f) corresponds to the auxiliary variable, j = 1......N2 represents the index variable, B k This represents the value of B obtained in the k-th iteration. B represents the k-th iteration. j The j-th block; Indicates the relationship with WP1(F) -1 f) is the corresponding auxiliary variable. Indicates and Corresponding auxiliary variables, Indicates to WP j (F -1 f) is the corresponding auxiliary variable. This represents the result obtained in the k-th iteration. Indicates the k-th iteration The j-th block, B 0 and They represent B and The initial value; This represents the Lagrange multiplier corresponding to the auxiliary variable B. Indicate u B The j-th block, Indicate u B The first block, Indicate u B The N2nd block, Indicate u B The initial value; Indicates with R C (f) is the corresponding auxiliary variable. Operator R C (·): This indicates extracting one from the multi-coil k-space. The matrix, operator R i (·): This indicates extraction from space for each coil k. The blocks are converted into column vectors according to the coil column vectors. The i-th operator, i = 1...N1, where N1 represents the space containing a single coil k. The number of blocks; Indicates with R C (f) corresponds to the auxiliary variable Z 0 and Representing Z and The initial value; denotes the Lagrange multiplier corresponding to the auxiliary variable z, denotes u Z the initial value of u; denotes the joint hard thresholding function, z denotes an input vector or matrix, θ denotes a threshold, a > 0, μ2> 0 are parameters; denotes the time acceleration factor, t 0 denotes the initial value of t; k denotes the number of iterations; S1: to perform singular value decomposition, i.e. Γ is a singular value matrix, Ψ and Φ represent left and right singular matrices respectively, and then an adaptive sparse transform matrix W of the image block in the k+1th iteration can be obtained k+1 , and the calculation formula is as follows: W k+1 = ΦΨ H where X k represents the matrix composed of the extracted image blocks at the kth iteration, B k is an auxiliary variable at the kth iteration, is B k at the kth iteration, and λk S2: Calculate the auxiliary variable for the (k+1)th iteration The calculation formula is as follows: where f k is the multi-coil k-space data reconstructed for the kth iteration, and are the values of the kth iteration, respectively. and is the kth iteration. S3: update the time acceleration factor t for the (k+1)th iteration k+1 The calculation formula is as follows: where t k denotes the time acceleration factor for the kth iteration; S4: update the auxiliary variable B for the (k+1)th iteration k+1 The calculation formula is as follows: wherein B k denotes an auxiliary variable of the kth iteration; S5: to perform singular value decomposition, i.e. U and V are left and right singular matrices, respectively, and Σ = diag(ω), ω contains all singular values of A, U, Σ and V are low-rank estimated using a low-rank estimation parameter r and truncated to obtain and the auxiliary variable of the k+1th iteration for: wherein, denotes an auxiliary variable of the kth iteration; S6: update the low-rank matrix auxiliary variable Z of the k+1th iteration k+1 The calculation formula is as follows: wherein Z k denotes an auxiliary variable for the kth iteration; S7: Calculate the k-space data f of the multi-coil image to be reconstructed at the (k+1)th iteration k+1 The calculation formula is as follows: wherein is the adjoint operator corresponding to R C (·) is the adjoint operator corresponding to R forms a diagonal LN x LN matrix, is the adjoint operator corresponding to P j (·) is the adjoint operator corresponding to P forms a diagonal LN x LN matrix, I is the LN x LN identity matrix and μ1 > 0 is a parameter; S8: Compute the matrix X synthesized from the extracted image blocks at the k+1 iteration k+1 The computation formula is as follows: S9: Compute B for the (k+1)th iteration k+1 The corresponding Lagrange multiplier The formula is as follows: S10: Calculate Z for the (k+1)th iteration k+1 the corresponding Lagrange multiplier The calculation formula is as follows: S11: f k+1 The single-coil amplitude reconstruction image is obtained by performing inverse Fourier transform and then square sum of square root (SOS) operation, and the calculation formula is as follows: wherein, represents the single-coil magnitude reconstruction image of the k+1th iteration; S12: Calculate x k+1 and the relative error RE between x k is calculated as follows: RE = ||x k+1 - x k ||2 / ||x k ||2 S13: judging whether an iteration stopping criterion is met, if a maximum iteration number K is reached, i.e. k=K or a condition RE<tol is met, entering step S14; otherwise, setting k=k+1 and returning to S1; wherein tol is a set tolerance; S14: output the final reconstructed single-coil magnitude image x k+1 .