A multi-slice parallel magnetic resonance imaging fast reconstruction method based on transform learning

Abstract

The present application relates to a kind of based on the fast reconstruction method of multi-slice parallel magnetic resonance imaging of transform learning, belong to magnetic resonance imaging technical field.The present application is based on multi-slice low rank tensor completion model (Multi-slice Hankel Tensor Completion, MS-HTC), combined with joint sparse transform learning (Joint Transform Learning, JTL) regular term, propose a kind of multi-slice parallel magnetic resonance imaging fast reconstruction model, named FMS-JTLHTC.The proposed method uses Alternating Direction Method of Multipliers (Alternating Direction Method of Multipliers, ADMM) to solve, then introduce Fast Iterative Shrinkage-Thresholding Algorithm (Fast Iterative Shrinkage-Thresholding Algorithm, FISTA) to accelerate convergence, and using Graphics Processing Unit (Graphics Processing Unit, GPU) to speed up algorithm.Experiment compares MS-HTC algorithm, experimental results show that the image reconstructed by FMS-JTLHTC algorithm has better reconstruction quality, and reconstruction speed is improved by 15 times.

Description

Technical Field

[0001] This invention relates to a fast reconstruction method for multi-slice parallel magnetic resonance imaging based on transform learning, belonging to the field of magnetic resonance imaging technology. Background Technology

[0002] Magnetic Resonance Imaging (MRI) is one of the modern medical imaging diagnostic tools. It does not emit ionizing radiation and can provide high-resolution images of soft tissues such as the brain, heart, liver, and blood vessels. In 2D multi-slice MRI, when the gaps between slices are very small (usually zero or a few millimeters), adjacent 2D slices are typically highly correlated. Utilizing the data redundancy between slices can effectively improve the reconstruction quality. However, since MRI data is acquired line by line along the phase-encoding direction in K-space, obtaining full-sample data takes a long time. Therefore, reducing MRI imaging time is of great practical significance. Parallel imaging techniques and the advent of compressed sensing (CS) have significantly shortened MR imaging scan time.

[0003] To reconstruct high-quality MR images by utilizing the redundancy between slices, researchers have proposed algorithms for joint multi-slice image reconstruction. These algorithms, in a non-uniformly undersampled slice sequence containing both highly undersampled slices (H-slices) and slightly undersampled slices (L-slices), estimate the missing K-space data of the H-slice from adjacent L-slices using interpolation. Then, based on the CS method, the interpolated H-slices and L-slices are reconstructed separately to obtain a multi-slice image. These algorithms are commonly referred to as interpolated compressed sensing (iCS) algorithms.

[0004] While algorithms like iCS offer excellent reconstruction performance and can achieve lower sampling rates, they do not utilize the multi-coil channel information of the image. Therefore, Deka and Datta proposed a CS pMRI reconstruction algorithm based on wavelet domain forest sparsity and Joint Total Variation (JTV), which utilizes both multi-slice and multi-coil channel information and performs parallel computation on multi-core CPUs and GP-GPUs. Liu et al. proposed the Multi-slice Hankel Tensor Completion (MS-HTC) model, which converts K-space data into block-wise Hankel tensors to utilize the similar coil sensitivity and similar image information between adjacent slices in multi-slice multi-coil MR data, reconstructing multi-slice MR images in a calibration-free manner.

[0005] While high-quality multi-slice MR images can be reconstructed by simultaneously utilizing information from data slices and between slices, significant aliasing artifacts exist in the reconstructed multi-slice images under high undersampling rates. Furthermore, the tensor decomposition and the conversion between K-space data and block Hankel tensors in the MS-HTC algorithm are computationally complex, resulting in long image reconstruction times and noticeable blur artifacts in the reconstructed images. Therefore, the reconstruction speed and quality of multi-slice MR images still need improvement. Summary of the Invention

[0006] The purpose of this invention is to overcome the shortcomings of existing multi-slice reconstruction techniques and provide a fast reconstruction method for multi-slice parallel magnetic resonance imaging based on transform learning, which can further improve the reconstruction quality and reconstruction speed of multi-slice magnetic resonance images.

[0007] This invention is achieved through the following technical solution: a fast reconstruction method for multi-slice parallel magnetic resonance imaging based on transform learning. It includes the following steps:

[0008] 1) Determine a multi-slice parallel magnetic resonance imaging reconstruction model based on transform learning;

[0009] 2) The alternating direction multiplier method is used to solve the reconstruction model. Specifically, singular value decomposition is used to solve the transformation learning matrix, joint hard thresholding is used to solve the sparse coding matrix, sequential truncated high-order singular value decomposition is used to solve the approximate low-rank tensor, fast iterative thresholding is used to accelerate the solved multi-slice parallel magnetic resonance K-space data, and GPU is used to accelerate the entire reconstruction model solution process.

[0010] 3) The square root of the sum of squares of the solved multi-slice parallel magnetic resonance K-space data is used to obtain the multi-slice parallel magnetic resonance image.

[0011] The specific steps are as follows:

[0012] S0: Initialization, f 0 =D H y, t 0 =0, k=0,

[0013] Among them, the superscript " 0 " represents the initial value before iteration, This represents the column-vectorized multi-coil K-space data to be reconstructed, i = 1...N C This represents the coil index variable, s = 1...N S N represents the slice index variable. CN represents the number of receiving coils used. S Indicates the number of slices. Superscript " H "" represents the conjugate transpose operation of a vector or matrix, and the spatial data of the i-th coil of the s-th slice of f is represented as: f 11 and Let f represent the first coil of the first slice and the Nth coil of the first slice of the column-vectorized multi-slice multi-coil K-space data f to be reconstructed. C K-space data of each coil, and Represent the Nth column of the multi-slice, multi-coil K-space data f to be reconstructed, which is vectorized into columns. S The first coil and the Nth coil of the slice C K-space data for each coil; It is vectorized K-space data of the same size as f. The spatial data of the i-th coil in the s-th slice is represented as follows: and Representing K-space data respectively The first coil of the first slice and the Nth C K-space data of each coil, and Representing K-space data respectively The Nth S The first coil and the Nth coil of the slice C K-space data of each coil, f 0 This represents the initial value of f. express initial value, This represents the multi-slice, multi-coil K-space undersampling operator. I represents the single-coil K-space undersampling matrix. S For N S ×N S The identity matrix, I C For N C ×N C The identity matrix, Indicates the Kronecker product. This represents undersampled multi-slice multi-coil K-space data, N = N x ×N y N represents the number of pixels in the single-coil image to be reconstructed. x N y 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 in the K-space data of a single slice. For multi-slice multi-coil Fourier transform operators, F x F yN x N y 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 represents extracting from a single coil image of a single slice. Image blocks are recombined into The matrix; P1(·) denotes the linear operator extracted from the first matrix, Indicates the Nth p Extracting linear operators from matrices, Let be the extraction operator for the s-th slice image, and let represent the matrix for extracting a single slice multi-coil image from a multi-slice multi-coil image; This indicates that the extracted image patches are combined into matrix A. s , A represents s The initial value of N p This indicates that a single slice or single coil image contains... The number of image patches; Represents the s-th slice Adaptive sparse transformation matrix of image patch express Point discrete cosine transform matrix, W s The initial value; Indicates with W s P1E s (F -1 f) is the corresponding auxiliary variable. Indicates and Corresponding auxiliary variables, Indicates with W s P j (F -1 f) is the auxiliary variable, j = 1...N p Indicates an index variable. B represents the value obtained in the k-th iteration. s The value, B represents the k-th iteration. s The j-th block; B s The initial value; Represents auxiliary variable B s The corresponding Lagrange multipliers, express The j-th block, express The first block, express The Nth p Each block, express Initial value; operator This indicates the construction of a K-space data from multiple slices and multiple coils. The third-order tensor; denoted as the joint hard thresholding function, z represents the input vector or matrix, θ represents the threshold, and α > 0, β > 0, μ1 > 0 are parameters; The time acceleration factor, t 0 t represents the initial value of t; k represents the number of iterations;

[0014] S1: Set s = 1;

[0015] S2: Yes Perform Singular Value Decomposition (SVD), that is... Σ s U is a singular value matrix. s and V s Let the left singular matrix and the right singular matrix represent the image patch in the (k+1)th iteration, respectively. Then, we can obtain the adaptive sparse transformation matrix of the image patch in that iteration. The calculation formula is as follows:

[0016]

[0017] in, Let represent the matrix composed of extracted image patches in the k-th iteration. It is the auxiliary variable in the k-th iteration. It is the kth iteration. The corresponding Lagrange multipliers;

[0018] S3: Calculate the auxiliary variables for the (k+1)th iteration. The calculation formula is as follows:

[0019]

[0020] Where f k For the multi-slice, multi-coil K-space data reconstructed in the k-th iteration, and They are respectively and The value of the kth iteration;

[0021] S4: Determine whether each slice has been completed. Update, when s = N S If the condition is met, proceed to step S5; otherwise, set s = s + 1 and return to S2.

[0022] S5: Calculation: right Expanding using mode 1 yields the matrix. It can be represented as:

[0023] S6: For the matrix Perform SVD to obtain By truncating the SVD using the rank parameter r1, the reconstructed matrix is ​​obtained.

[0024]

[0025] Where Σ (1) U is a singular value matrix. (1) and V (1) Let these represent the left singular matrix and the right singular matrix, respectively. and Not for U (1) , Σ (1) and V (1) The matrix obtained by rank truncation.

[0026] S7: Will Reduced to a third-order tensor It can be represented as: right Expanding in mode 2 yields the matrix. It can be represented as:

[0027] S8: For the matrix Perform SVD to obtain By truncating the SVD using the rank parameter r2, the reconstructed matrix is ​​obtained.

[0028]

[0029] Where Σ (2) U is a singular value matrix. (2) and V (2) Let these represent the left singular matrix and the right singular matrix, respectively; and They are respectively for U (2) , Σ (2) and V (2) The matrix obtained by rank truncation.

[0030] S9: Will Reduced to a third-order tensor It can be represented as:

[0031] S10: Calculate the K-space data for the (k+1)th iteration. The calculation formula is as follows:

[0032]

[0033] Where, in the formula It is the adjoint operator of R(·), R * R is a diagonal matrix; It is P j The adjoint operator of (·), Form an N C N×N C A diagonal matrix of N, I is a diagonal matrix of N. C N×N C The identity matrix of N, It is E s The adjoint operator of (·), It is a diagonal matrix, and α>0 and μ1>0 are parameters;

[0034] S11: Update the time acceleration factor t for the (k+1)th iteration. k+1 The calculation formula is as follows:

[0035]

[0036] Among them, t k This represents the time acceleration factor for the k-th iteration;

[0037] S12: Update the K-space data f of the multi-slice multi-coil image to be reconstructed in the (k+1)th iteration. k+1 The calculation formula is as follows:

[0038]

[0039] in, This represents the K-space data of the multi-slice, multi-coil image in the k-th iteration.

[0040] S13: Set s = 1;

[0041] S14: Calculate the matrix synthesized from the extracted image patches in the (k+1)th iteration. The calculation formula is as follows:

[0042]

[0043] S15: Calculate the (k+1)th iteration The corresponding Lagrange multipliers The calculation formula is as follows:

[0044]

[0045] S16: Determine whether each slice image has been completed. Update, when s = N S If the condition is met, proceed to step S17; otherwise, set s = s + 1 and return to S14.

[0046] S17: f k+1 Perform an inverse Fourier transform, followed by a square root of sum-of-squares (SOS) operation, to obtain the reconstructed multi-slice amplitude image. The calculation formula is as follows:

[0047]

[0048] in, This represents the multi-slice reconstructed image in the (k+1)th iteration;

[0049] S18: Calculation and The relative error (RE) between them is calculated using the following formula:

[0050]

[0051] S19: Determine whether the iteration stopping criterion is met. If the maximum number of iterations K is reached, i.e., k = K, or RE is less than tol, proceed to step S20; otherwise, let k = k + 1 and return to S1, where tol is the set tolerance.

[0052] S20: Output the final reconstructed multi-slice image.

[0053] The beneficial effects of this invention are as follows: In this invention, a JTL regularization term is introduced into the MS-HTC model, and ADMM is used to decompose the problem into several sub-problems for separate solving. Secondly, FISTA is used to accelerate the convergence speed of the proposed algorithm, and finally, GPU parallel computing is used to accelerate the reconstruction speed. This results in a fast multi-slice magnetic resonance imaging reconstruction algorithm, FMS-JTLHTC. Theoretical analysis and imaging experiments show that, compared with MS-HTC, the fast multi-slice magnetic resonance imaging reconstruction method FMS-JTLHTC proposed in this invention can reconstruct images of better quality and significantly reduce the reconstruction time of multi-slice magnetic resonance imaging. Attached Figure Description

[0054] Figure 1 This is a flowchart of the method of the present invention;

[0055] Figure 2 Four consecutive slice images (i.e., dataset1) were used to acquire in vivo human brain data of the subject by using 12-channel coils for full sampling.

[0056] Figure 3 Four one-dimensional variable-density undersampling masks with a 5× sampling rate but different sampling locations;

[0057] Figure 4 and Figure 5 The images reconstructed from one-dimensional variable-density undersampled data using the FMS-JTLHTC algorithm and the MS-HTC algorithm are shown in dataset1.

[0058] Figure 6 and Figure 7 The text shows the corresponding Figure 4 and Figure 5 Error image. Detailed Implementation

[0059] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings.

[0060] Example 1: This invention is an efficient reconfiguration method based on the MS-HTC framework.

[0061] Assuming, This represents the column-vectorized multi-coil K-space data to be reconstructed, i = 1...N C This represents the coil index variable, s = 1...N S N represents the slice index variable. C N represents the number of receiving coils used. S Indicates the number of slices. Superscript " H "" represents the conjugate transpose operation of a vector or matrix, and the spatial data of the i-th coil of the s-th slice of f is represented as: f 11 and Let f represent the first coil of the first slice and the Nth coil of the first slice of the column-vectorized multi-slice multi-coil K-space data f to be reconstructed. C K-space data of each coil, and Represent the Nth column of the multi-slice, multi-coil K-space data f to be reconstructed, which is vectorized into columns. S The first coil and the Nth coil of the slice C K-space data for each coil; It is vectorized K-space data of the same size as f. The spatial data of the i-th coil in the s-th slice is represented as follows: and Representing K-space data respectively The first coil of the first slice and the Nth C K-space data of each coil, and Representing K-space data respectively The Nth S The first coil and the Nth coil of the slice C K-space data of each coil, f 0 This represents the initial value of f. express initial value, This represents the multi-slice, multi-coil K-space undersampling operator. I represents the single-coil K-space undersampling matrix. S For N S ×N S The identity matrix, I C For N C ×N C The identity matrix, Indicates the Kronecker product. This represents undersampled multi-slice multi-coil K-space data, N = N x ×N y N represents the number of pixels in the single-coil image to be reconstructed. x N y These represent the number of pixels in the horizontal and vertical directions of the magnetic resonance image, respectively; M represents the actual number of sampling points in the spatial data of a single coil in a single slice; operator This indicates the construction of a K-space data from multiple slices and multiple coils. A third-order tensor. Introducing auxiliary variables. The MS-HTC reconstruction model can then be obtained as follows:

[0062]

[0063] Where ||·||2 represents the L2 norm of the vector, ||·|| F Let rank(Γ) represent the Frobenius norm of the tensor, where rank(Γ) = (r1, r2, N). S Let r1, r2, and N be the ranks of the three modal expansion matrices of the third-order tensor Γ. S .

[0064] To improve the reconstruction quality of 2D multi-slice MRI, this invention combines the JTL regularization term with the MS-HTC model to propose the FMS-JTLHTC model, resulting in the following optimization problem:

[0065]

[0066] Where α and β are regularization parameters, Extracting a linear operator for the j-th matrix represents extracting from a single coil image of a single slice. Image blocks are recombined into Matrix; Represents the s-th slice Adaptive sparse transformation matrix of image patch W s The conjugate transpose matrix, H " represents the conjugate transpose operation of a matrix; Let be the extraction operator for the s-th slice image, and let represent the matrix for extracting a single slice multi-coil image from a multi-slice multi-coil image; For multi-slice multi-coil Fourier transform operators, F x F y N x N y Point Fourier transform matrix, I S For N S ×N S The identity matrix, I C For N C ×N C The identity matrix, F represents the Kronecker product. -1 Indicates the inverse Fourier transform, -1 " indicates the inverse operation of a matrix.

[0067] Introducing auxiliary variables and the corresponding Lagrange multipliers Using ADMM technology, problem (2) can be transformed into the following sub-problems:

[0068]

[0069]

[0070]

[0071]

[0072]

[0073] In formulas (3)-(7), 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), Let μ1 represent the adaptive sparse transformation matrix of the image patch in the (k+1)th iteration within the s-th slice, where μ1 > 0 is a parameter. Indicates with W s P1E s (F -1 f) is the corresponding auxiliary variable. Indicates with W s P j E s (F-1 f) is the auxiliary variable, j = 1...N p Indicates an index variable. B represents the value obtained in the k-th iteration. s The value, B represents the k-th iteration. s The j-th block; Represents auxiliary variable B s The corresponding Lagrange multipliers, express The j-th block, f k Let K represent the K-space data of the multi-coil image to be reconstructed in the k-th iteration; in formula (4), This represents the (k+1)th iteration. The j-th block; in formula (5), Γ k+1 Let f denote the low-rank approximate tensor Γ for the (k+1)th iteration; in formula (6), f k+1 Let K represent the K-space data of the multi-slice multi-coil image to be reconstructed in the (k+1)th iteration; in formula (7), This represents the (k+1)th iteration. The j-th block.

[0074] Next, we will solve the above sub-problems.

[0075] Solving transformation learning W s :because as well as Problem (3) can then be transformed into the following form:

[0076]

[0077] in, Let represent the matrix composed of extracted image patches in the k-th iteration. Describes the auxiliary variable for the k-th iteration. Indicates the k-th iteration The corresponding Lagrange multipliers, ||·|| F This represents the Frobenius norm.

[0078] Therefore, problem (8) can be solved using the Singular Value Decomposition (SVD) method:

[0079]

[0080]

[0081] Where, Σ s U is a singular value matrix. s and Vs Let represent the left singular matrix and the right singular matrix, respectively.

[0082] Solving sparse coding B sj Problem (4) can be simplified to:

[0083]

[0084] Problem (11) can be solved using the joint hard thresholding method, and its solution can be expressed as:

[0085]

[0086] Where H J (z,θ) denotes the joint hard thresholding operator, defined as follows:

[0087]

[0088] Where z is the input vector or matrix, and θ is the threshold.

[0089] according to Then formula (12) can be expressed as:

[0090]

[0091] in, This represents the auxiliary variable for the (k+1)th iteration.

[0092] Solving for the low-rank approximate tensor Γ: The solution to equation (5) can be obtained by using Sequentially Truncated Higher-Order Singular Value Decomposition (ST-HOSVD) on R(f k The decomposition yields the following result, which can be represented in the following form:

[0093] Γ k+1 =ST-HOSVD(R(f k (15)

[0094] The specific operation is as follows: First, R(f) k Expanding using Mode 1 yields the matrix. Then on the matrix Perform SVD to obtain The SVD is truncated using the rank parameter r1 to obtain the matrix. Where Σ (1) U is a singular value matrix. (1) and V (1) Let these represent the left singular matrix and the right singular matrix, respectively; and Not for U(1) , Σ (1) and V (1) The matrix obtained by rank truncation. Then... Reduced to a third-order tensor right Expanding in mode 2 yields the matrix. For matrix Perform SVD to obtain The SVD is truncated using the rank parameter r2 to obtain the matrix. Where Σ (2) U is a singular value matrix. (2) and V (2) Let these represent the left singular matrix and the right singular matrix, respectively; and They are respectively for U (2) , Σ (2) and V (2) The matrix obtained by rank truncation. Finally, Reduced to a third-order tensor

[0095] Update image f: Set the derivative of the objective function of problem (6) with respect to f to zero, and we get the following form:

[0096]

[0097] Where α>0 and μ1>0 are parameters. It is the adjoint operator of R(·), R * R is a diagonal matrix; It is P j The adjoint operator of (·), Form an N C N×N C A diagonal matrix of N, I is a diagonal matrix of N. C N×N C The identity matrix of N, It is E s The adjoint operator of (·), It is a diagonal matrix, and

[0098] Therefore, problem (16) can be simplified to:

[0099]

[0100] Among them, f k+1 This represents the K-space data of the multi-slice multi-coil iteration in the (k+1)th iteration.

[0101] Using the FISTA algorithm to accelerate f, we have:

[0102]

[0103]

[0104]

[0105] in, and Let t represent the K-space data of the (k+1)th and (k)th iterations, respectively. k+1 and t k These represent the time acceleration factors for the (k+1)th and kth iterations, respectively.

[0106] Solve for B s Corresponding Lagrange multipliers Because there is as well as Formula (7) can then be expressed as:

[0107]

[0108] f k+1 Perform an inverse Fourier transform, followed by a SOS (Sum Of Squares) operation to obtain a single amplitude image. This operation can be represented as:

[0109]

[0110] In summary, this invention presents a novel fast multi-slice parallel magnetic resonance imaging reconstruction method, FMS-JTLHTC. Finally, the method is terminated when the relative error (RE) is less than the tolerance tol or the maximum number of iterations K is reached. The definition of RE is as follows:

[0111]

[0112] The specific process is as follows: Figure 1 As shown, the steps are as follows:

[0113] S0: Initialization, f 0 =D H y, t 0 =0, k=0,

[0114] superscript " 0 " represents the initial value before iteration, f 0 and They represent f and initial value, Indicate u B initial value, A representss initial value, W s initial value, B s The initial value of t 0 The initial value of the time acceleration factor t is represented, and k represents the number of iterations;

[0115] S1: Set s = 1;

[0116] S2: Yes Perform singular value decomposition, i.e. Σ s U is a singular value matrix. s and V s Let the left singular matrix and the right singular matrix represent the image patch in the (k+1)th iteration, respectively. Then, we can obtain the adaptive sparse transformation matrix of the image patch in that iteration. Formula (10) is shown below;

[0117] S3: Calculate the auxiliary variables for the (k+1)th iteration. Formula (14) is shown below;

[0118] S4: Determine whether each slice has been completed. Update, when s = N S If the condition is met, proceed to step S5; otherwise, set s = s + 1 and return to S2.

[0119] S5: Calculation: right Expanding using mode 1 yields the matrix. It can be represented as:

[0120] S6: For the matrix Perform SVD to obtain By truncating the SVD using the rank parameter r1, the reconstructed matrix is ​​obtained.

[0121]

[0122] Where Σ (1) U is a singular value matrix. (1) and V (1) Let these represent the left singular matrix and the right singular matrix, respectively; and Not for U (1) , Σ (1) and V (1) The matrix obtained by rank truncation.

[0123] S7: Will Reduced to a third-order tensor It can be represented as: right Expanding in mode 2 yields the matrix. It can be represented as:

[0124] S8: For the matrix Perform SVD to obtain By truncating the SVD using the rank parameter r2, the reconstructed matrix is ​​obtained.

[0125]

[0126] Where Σ (2) U is a singular value matrix. (2) and V (2) Let these represent the left singular matrix and the right singular matrix, respectively; and They are respectively for U (2) , Σ (2) and V (2) The matrix obtained by rank truncation.

[0127] S9: Will Reduced to a third-order tensor It can be represented as:

[0128] S10: Calculate the K-space data for the (k+1)th iteration. Formula (18) is shown below;

[0129] S11: Update the time acceleration factor t for the (k+1)th iteration. k+1 The formula is as shown in (19);

[0130] S12: Update the K-space data f of the multi-slice multi-coil image to be reconstructed in the (k+1)th iteration. k+1 The formula is as shown in (20);

[0131] S13: Set s = 1;

[0132] S14: Calculate the matrix composed of the extracted image patches in the (k+1)th iteration. The formula is as follows:

[0133]

[0134] S15: Calculate the (k+1)th iteration The corresponding Lagrange multipliers Formula (21) is shown below;

[0135] S16: Determine whether each slice image has been completed. Update, when s = NS If the condition is met, proceed to step S17; otherwise, set s = s + 1 and return to S14.

[0136] S17: f k+1 Perform an inverse Fourier transform, followed by a square root of sum-of-squares (SOS) operation, to obtain a multi-slice single-coil amplitude reconstruction image. Formula (22) is shown below;

[0137] S18: Calculation and The relative error (RE) between them is given by formula (23);

[0138] S19: Determine whether the iteration stopping criterion is met. If the maximum number of iterations K is reached, i.e., k = K, or RE is less than tol, proceed to step S20; otherwise, let k = k + 1 and return to S1, where tol is the set tolerance.

[0139] S20: Output the final reconstructed multi-slice single-coil amplitude image.

[0140] Experimental results:

[0141] To verify the reconstruction performance of the proposed algorithm, this invention compares it with the MS-HTC algorithm. All algorithms were implemented using MATLAB (MathWorks, Natick, MA), and all experiments were conducted on a desktop computer configured with an Intel(R) Core(TM) i9-13900K@3.00GHz processor, 64GB of memory, and an NVIDIA GeForce RTX 4090 GPU.

[0142] To compare the reconstruction performance of the algorithms, this invention used four consecutive slices of a human brain map for simulation experiments, named dataset1 (e.g., ...). Figure 2 This invention generates four one-dimensional variable-density undersampling masks with an acceleration factor of 5 but different sampling positions for four slices to artificially undersampling the fully sampled dataset, thereby generating a test dataset. The undersampling masks are as follows: Figure 3 As shown.

[0143] In the experiments, Peak Signal-to-Noise Ratio (PSNR) and Feature Similarity Index (FSIM) were used to evaluate the reconstruction performance of the method. Higher PSNR and FSIM values ​​indicate higher image reconstruction quality. The metrics were calculated within the region of interest. All experiments involved adjusting the method parameters to obtain the optimal PSNR value.

[0144] The formula for calculating PSNR is as follows:

[0145]

[0146] Where MSE is the reconstructed image The mean square error between the reference image x and the reference image x, where max is the maximum pixel value of x.

[0147] FSIM's following:

[0148]

[0149] Among them, PC x and It is the phase consistency between the reference image and the reconstructed image, S PC It is a PC x and Similarity measure, S G It is a measure of the similarity of gradient magnitude values ​​between the reference image and the reconstructed image.

[0150] This invention compares the visual representations reconstructed from a selected dataset using two methods when the acceleration factor is 5. Figure 4 and Figure 5 The images reconstructed using the MS-HTC and FMS-JTLHTC algorithms are shown in dataset1. Figure 4 and Figure 5 It can be seen that the FMS-JTLHTC algorithm more effectively removes blur artifacts and preserves image details, resulting in better reconstruction quality compared to the MS-HTC algorithm.

[0151] This invention is in Figure 6 and Figure 7 The image shows the error graphs of the MS-HTC and FMS-JTLHTC algorithms on dataset1 (white dots in the error graph represent errors; fewer white dots indicate better reconstruction quality). Figure 6 and Figure 7 It can be seen that the FMS-JTLHTC algorithm has fewer white dots in its error map compared to the MS-HTC algorithm, indicating that the FMS-JTLHTC algorithm reconstructs images of better quality.

[0152] Furthermore, this invention also analyzes the numerical PSNR of the reconstructed images from the two methods. The reconstructed image from the MS-HTC algorithm (e.g.) Figure 4 The PSNR of the image was 32.49 dB, and the FSIM was 0.9969. The reconstructed image using the FMS-JTLHTC algorithm (e.g.) Figure 5 The PSNR of the proposed FMS-JTLHTC algorithm is 34.93 dB, and the FSIM is 0.9986. It can be seen that, under the evaluation criteria, the proposed FMS-JTLHTC algorithm has a PSNR value 2.44 dB higher and an FSIM value 0.0017 higher than the MS-HTC algorithm, further validating the results of the visual comparative analysis.

[0153] Finally, although the proposed FMS-JTLHTC algorithm has good reconstruction quality, its reconstruction speed is slow. Therefore, this invention utilizes a GPU to accelerate the proposed FMS-JTLHTC algorithm, resulting in the accelerated version FMSJTLHTC-GPU, which significantly shortens the running time of the method. For dataset1, the reconstruction running time of the FMSJTLHTC-GPU algorithm and the MS-HTC algorithm are 17.4 seconds and 272.5 seconds, respectively. The reconstruction speed of the FMSJTLHTC-GPU algorithm is 15 times that of the MS-HTC algorithm, making it more valuable for practical applications.

[0154] This invention introduces JTL regularization terms into MS-HTC, proposing a fast 2D multi-slice MRI reconstruction algorithm, FMS-JTLHTC. Experimental results on multiple sets of experiments show that the proposed FMS-JTLHTC algorithm significantly improves reconstruction quality compared to the MSHTC algorithm, and also increases reconstruction speed by 15 times.

[0155] The specific embodiments of the present invention have been described in detail above with reference to the accompanying drawings. However, the present invention is not limited to the above embodiments. Within the scope of knowledge possessed by those skilled in the art, various changes can be made without departing from the spirit of the present invention.

Claims

1. A fast reconstruction method for multi-slice parallel magnetic resonance imaging based on transform learning, characterized in that: Includes the following steps: 1) Determine a multi-slice parallel magnetic resonance imaging reconstruction model based on transform learning; 2) The alternating direction multiplier method is used to solve the reconstruction model. Specifically, singular value decomposition is used to solve the transformation learning matrix, joint hard thresholding is used to solve the sparse coding matrix, sequential truncated high-order singular value decomposition is used to solve the approximate low-rank tensor, fast iterative thresholding is used to accelerate the solved multi-slice parallel magnetic resonance K-space data, and GPU is used to accelerate the entire reconstruction model solution process. 3) Obtain the multi-slice parallel magnetic resonance image by taking the square root of the sum of squares of the solved multi-slice parallel magnetic resonance K-space data; The specific steps are as follows: S0: Initialization, , , , , , , , ; Among them, the superscript " " represents the initial value before iteration, This represents the column-vectorized multi-coil K-space data to be reconstructed. Indicates the coil index variable. Represents the slice index variable. Indicates the number of receiving coils used. Indicates the number of slices, indicated by the superscript " " represents the conjugate transpose operation of a vector or matrix. The The first slice The spatial data of each coil K is represented as follows: , and These represent the column-vectorized, multi-slice, multi-coil K-space data to be reconstructed. The first slice, the first coil, and the first... K-space data of each coil, and These represent the column-vectorized, multi-slice, multi-coil K-space data to be reconstructed. The The first coil of the slice and the first... K-space data for each coil; Is and Vectorized K-space data of the same size, The The first slice The spatial data of each coil K is represented as follows: , and Representing K-space data respectively The first coil of the first slice and the first... K-space data of each coil, and Representing K-space data respectively The The first coil of the slice and the first... K-space data of each coil, express initial value, express initial value, This represents the multi-slice, multi-coil K-space undersampling operator. This represents the K-space undersampling matrix of a single coil. for The identity matrix, for The identity matrix, Indicates the Kronecker product. This represents undersampled multi-slice multi-coil K-space data. This represents the number of pixels in the single-coil image to be reconstructed. , These represent the number of pixels in the horizontal and vertical directions of the magnetic resonance image, respectively. This represents the actual number of sampling points for a single coil in the K-space data of a single slice; For multi-slice, multi-coil Fourier transform operators, , They represent , Point Fourier transform matrix, Indicates the inverse Fourier transform, " indicates the inverse operation of a matrix; For the first A matrix extraction linear operator is used to extract linear data from a single coil image of a single slice. Image blocks are recombined into Matrix; This indicates that the first matrix extracts the linear operator. Indicates the first Extracting linear operators from matrices, For the first The extraction operator for each slice image represents a matrix for extracting a single slice multi-coil image from a multi-slice multi-coil image; This indicates that the extracted image patches are combined into a matrix. , express initial value, This indicates that a single slice or single coil image contains... The number of image patches; Indicates the first A slice Adaptive sparse transformation matrix of image patch express Point discrete cosine transform matrix, express The initial value; , Indicates and Corresponding auxiliary variables, Indicates and Corresponding auxiliary variables, Indicates and Corresponding auxiliary variables, Indicates an index variable. Indicates the first The result of the next iteration The value, Indicates the first The next iteration The Each block; express The initial value; Representing auxiliary variables The corresponding Lagrange multipliers, express The Each block, express The first block, express The Each block, express Initial value; operator This indicates the construction of a K-space data from multiple slices and multiple coils. The third-order tensor; Represents the joint hard thresholding function. Represents the input vector or matrix. Indicates the threshold. , , For parameters; Represents the time acceleration factor. express The initial value; Indicates the number of iterations; S1: Settings ; S2: Yes Perform singular value decomposition, i.e. , It is a singular value matrix. and Let represent the left singular matrix and the right singular matrix respectively, then we can obtain the... The adaptive sparse transformation matrix of the image patch in the next iteration The calculation formula is as follows: ； in, Indicates the first The matrix formed by combining the extracted image patches in the next iteration. It is the first Auxiliary variables for the next iteration It is the first The next iteration The corresponding Lagrange multipliers; S3: Calculate the... Auxiliary variables for the next iteration The calculation formula is as follows: ； in For the first The multi-slice, multi-coil K-space data reconstructed in the next iteration and They are respectively and No. The value of the next iteration; S4: Determine whether each slice has been completed. Update, when Proceed to step S5 if the condition is met; otherwise, let... Then return to S2; S5: Calculation: ,right Expanding using mode 1 yields the matrix. , can be represented as: ; S6: For the matrix Perform SVD to obtain Using the rank parameter By truncating the SVD, the reconstructed matrix is ​​obtained. : ； in It is a singular value matrix. and Let these represent the left singular matrix and the right singular matrix, respectively; , and Don't be right , and The matrix obtained by rank truncation; S7: Will Reduced to a third-order tensor , can be represented as: ,right Expanding in mode 2 yields the matrix. , can be represented as: ; S8: For the matrix Perform SVD to obtain Using the rank parameter By truncating the SVD, the reconstructed matrix is ​​obtained. : ； in It is a singular value matrix. and Let these represent the left singular matrix and the right singular matrix, respectively; , and They are respectively for , and The matrix obtained by rank truncation; S9: Will Reduced to a third-order tensor , can be represented as: ; S10: Calculate the... K-space data of the next iteration The calculation formula is as follows: ； Where, in the formula yes The adjoint operator, It is a diagonal matrix; yes The adjoint operator, Form a diagonal matrix, for The identity matrix, yes The adjoint operator, It is a diagonal matrix, and , , For parameters; S11: Update Time acceleration factor for the next iteration The calculation formula is as follows: ； in, Indicates the first The time acceleration factor for each iteration; S12: Update K-space data of the multi-slice, multi-coil image to be reconstructed in the next iteration The calculation formula is as follows: ； in, Indicates the first K-space data of multi-slice, multi-coil images in the next iteration; S13: Settings ; S14: Calculate the... The matrix synthesized from the extracted image patches in the next iteration The calculation formula is as follows: ； S15: Calculate the... The next iteration The corresponding Lagrange multipliers The calculation formula is as follows: ； S16: Determine whether each slice image has been completed. Update, when Proceed to step S17 if the condition is met; otherwise, let Then return to S14; S17: Will Perform an inverse Fourier transform, followed by a square root SOS operation on the sum of squares to obtain the reconstructed multi-slice amplitude image. The calculation formula is as follows: ； in, Indicates the first Image reconstruction via multiple slices in subsequent iterations; S18: Calculation and The relative error RE between them is calculated using the following formula: ； S19: Determine if the iteration stopping criterion is met; if the maximum number of iterations has been reached... ,Right now or satisfy RE Less than Proceed to step S20; otherwise, let Return to S1, where, Tolerance set; S20: Output the final reconstructed multi-slice image. .

Images