A fast parallel imaging reconstruction method based on SIDWT and iterative self-consistency

Abstract

The present invention relates to a fast parallel imaging reconstruction method based on SIDWT and iterative self-consistency, and belongs to the field of magnetic resonance imaging technology. Parallel magnetic resonance imaging reconstruction has been a research hotspot in recent years. Based on the SIDWT transform, the present invention introduces the L1 norm regularization term of the transform coefficient into the iterative self-consistency parallel imaging reconstruction model in the K-space domain, and proposes an efficient reconstruction method named fSIDWT‑SPIRiT. This method aims at complex optimization problems containing data consistency terms, calibration consistency terms and L1 norm regularization terms. The data consistency terms and calibration consistency terms are first merged, and then the projected fast iterative thresholding method (projected Fast Iterative Shrinkage‑Thresholding Algorithm, pFISTA) is used to solve them. Experimental results show that the proposed fSIDWT‑SPIRiT method not only provides higher image reconstruction quality, but also effectively improves the convergence speed and reduces the image reconstruction time.

Description

Technical Field

[0001] The invention relates to a fast parallel imaging reconstruction method based on SIDWT and iterative self-consistency, and belongs to the technical field of magnetic resonance imaging. Background Art

[0002] Magnetic resonance imaging (MRI) is a widely used and indispensable tool in medical diagnosis due to its advantages, including multi-directional imaging, lack of ionizing radiation, and non-invasiveness and harmlessness to human tissue. However, its long scanning times, slow imaging speeds, and high costs, combined with the image artifacts and low temporal resolution of dynamic MRI caused by patient movement during long scans, have limited its application. Therefore, shortening scanning and imaging times and reconstructing higher-quality images have been hot topics in MRI research.

[0003] To achieve fast imaging, parallel imaging (PI) and compressed sensing (CS) are two effective methods. PI uses a coil array with multiple receiving channels to simultaneously acquire partial K-space data. By leveraging the correlation between the coils, high-quality images can be reconstructed, thereby accelerating MRI scans. CS theory, unlike traditional sampling theorems, exploits signal redundancy to break through the Nyquist sampling rate, allowing a signal very close to the original to be recovered from a small number of measured values. Therefore, applying CS to MRI can significantly shorten scan times. Since PI and CS imaging methods are based on different prior knowledge and both achieve fast imaging, they are often combined to achieve even faster imaging.

[0004] After in-depth research on pMRI imaging methods, Lustig et al. proposed the iterative self-consistent parallel imaging reconstruction (SPIRiT) model, which is a generalized reconstruction framework based on self-consistency. The reconstruction problem is formulated as a sparse optimization problem, which is solved by calibration consistency and data acquisition consistency. The image can be accurately reconstructed from undersampled data of any K-space sampling pattern.

[0005] To improve MRI reconstruction quality or shorten reconstruction time, many researchers have conducted a series of studies and improvements based on the SPIRiT model. Duan et al. proposed a fast reconstruction method based on SPIRiT by adding joint total variation (JTV) and joint L1 norm (JL1) regularization constraints to the SPIRiT model. This method combines the calibration consistency term and the data consistency term into one term, then uses operator splitting (OS) to decompose it into two easily computable subproblems, which are solved using the split Bregman denoising method. Finally, the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) is used for acceleration. This method not only ensures image reconstruction quality but also significantly reduces image reconstruction time. Zhang et al. introduced the L1 norm regularization term into the SPIRiT model and proposed a reconstruction method, namely the pFISTA-SPIRiT method, which uses the projected Fast Iterative Shrinkage-Thresholding Algorithm (pFISTA) to solve the SPIRiT problem under different sampling modes and different tight frames. This method has only one adjustable parameter, and the reconstruction error is insensitive to the method parameters. It also allows the widespread use of different tight frames in MRI image reconstruction. This method achieves better reconstruction and converges faster.

[0006] In summary, adding a regularization term to the reconstruction model can effectively improve the method's reconstruction quality, and merging the calibration consistency term and data consistency term in the reconstruction model can effectively speed up the method's convergence. However, existing SPIRiT-based improved methods only introduce a simple L1-norm regularization term and directly use pFISTA to solve in the image domain, leaving room for improvement in both reconstruction speed and quality. Summary of the Invention

[0007] The present invention provides a fast parallel imaging reconstruction method based on SIDWT and iterative self-consistency, aiming to overcome the shortcomings of the existing technology, further improve the reconstruction quality of magnetic resonance imaging and effectively shorten the reconstruction time.

[0008] The technical solution adopted by the present invention is: a fast parallel imaging reconstruction method based on SIDWT and iterative self-consistency, comprising the following steps:

[0009] S0: Initialization, let x 0 =0,z 0 =0,t0 =0, k=0;

[0010] Among them, the superscript " 0 " indicates the initial value, represents the multi-coil K-space data to be reconstructed, Indicates the K-space data of the first coil, Indicates the Qth coil K-space data, Q represents the number of receiving coils, represents the qth coil K-space data, q=1,...,Q represents the coil index variable, (·) H Represents the conjugate transpose operation of a vector or matrix, N = n y ×n x Indicates the number of pixels in a single image, n y and n x Represents the number of rows and columns of a single image, x 0 represents the initial value of x; represents the intermediate variable, z 0 represents the initial value of z; represents the factor related to acceleration, t 0 represents the initial value of t; k represents the number of iterations;

[0011] S1: Calculate the intermediate variable u related to the gradient problem k , the calculation formula is as follows:

[0012]

[0013] Among them, z k represents the intermediate variable obtained after the kth iteration, L represents The Lipschitz constant of the gradient, G represents the frequency domain self-consistent convolution operator obtained from the multi-coil K-space self-calibration region, and I represents a QN×QN identity matrix;

[0014] S2: intermediate variables u related to the gradient problem k Calibrate to get The calculation formula is as follows:

[0015]

[0016] in, and denote operators for selecting sampled and unsampled points from the multi-coil K-space, and denote operators for placing sampled and unsampled points back to their correct positions in the multi-coil K-space, respectively. T represents the transpose operation of the matrix, represents undersampled multi-coil K-space data, M represents the number of points actually sampled in single-coil K-space data, M<<N;

[0017] S3: Calculate the multi-coil K-space data x reconstructed after the k+1th iteration k+1 , the calculation formula is as follows:

[0018]

[0019] in, represents the coil-by-coil Fourier transform, I Q represents a Q×Q identity matrix, represents the Kronecker product, represents the two-dimensional Fourier transform, F x and F y Represents n x and n y Point Fourier transform matrix, Represents coil-by-coil sparse transform, which is used to thin out the image. It is a tight framework. In this invention, SIDWT is selected as the tight framework in the experiment. * represents the adjoint of Ψ, and specifically satisfies Ψ * Ψ=I,(·) * represents the adjoint operation of the matrix, represents the point-by-point soft threshold function, λ>0 represents the regularization parameter, represents the coil-by-coil inverse Fourier transform, (·) -1 Represents the inverse operation of the matrix;

[0020] S4: Update the intermediate variable t of the kth iteration k+1 and z k+1 , the calculation formulas are as follows:

[0021]

[0022]

[0023] Among them, z k+1 represents the intermediate variable obtained after the k+1th iteration, t k+1 and t k They represent the acceleration factors for the k+1th and kth iterations, respectively, and x k represents the multi-coil K-space data reconstructed after the k-th iteration;

[0024] S5: Determine whether the number of iterations k reaches the maximum number of iterations. If so, proceed to step S6; otherwise, set k = k + 1 and return to step S1;

[0025] S6: Reconstruct the multi-coil K-space data x k+1Each coil is subjected to inverse Fourier transform, and then the square root of sum of squares (SOS) method is used to obtain the final reconstructed image. The calculation formula is as follows:

[0026]

[0027] The present invention provides the following benefits: Based on the SIDWT, a fast parallel MRI reconstruction method, fSIDWT-SPIRiT, is proposed that combines an L1-norm regularization term with the SPIRiT model to accelerate reconstruction speed and improve reconstruction quality. This method addresses complex optimization problems involving data consistency, calibration consistency, and an L1-norm regularization term by first merging the data and calibration consistency terms into a single term, then solving the problem in the K-space domain using SPIRiT using the pFISTA technique. Simulation results demonstrate that the proposed method effectively improves the quality of reconstructed images and achieves faster convergence. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] Figure 1 is a flow chart of the method of the present invention;

[0029] Figure 2 The heart slices of the subjects were acquired using a 12-channel receiving coil for full sampling (i.e., dataset 1);

[0030] Figure 3 is a 2D Poisson disk undersampling mask with 3× sampling rate and 24×24 ACS lines;

[0031] Figure 4-Figure 5 The images reconstructed from undersampled data of a 2D Poisson disk with 3× acceleration and 24×24 center calibration are shown in dataset 1 using the pFISTA-SPIRiT method and the fSIDWT-SPIRiT method.

[0032] Figure 6-Figure 7 The corresponding Figure 4 and Figure 5 Error map of

[0033] Figure 8 Shows the corresponding Figure 4 and Figure 5 Comparison of reconstruction speed. DETAILED DESCRIPTION

[0034] The technical solution of the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0035] Example 1: The present invention proposes a fast parallel imaging reconstruction method based on SIDWT and iterative self-consistency based on the SPIRiT framework.

[0036] Assumptions represents the multi-coil K-space data to be reconstructed, Indicates the K-space data of the first coil, Indicates the Qth coil K-space data, Q represents the number of receiving coils, represents the qth coil K-space data, q=1,...,Q represents the coil index variable, (·) H Represents the conjugate transpose operation of a vector or matrix, N = n y ×n x Indicates the number of pixels in a single image, n y and n x Represent the number of rows and columns of a single image respectively, then the undersampled data obtained from all coils It is given by:

[0037] y=Dx (1)

[0038] in, Indicates the operator for selecting sampled points from the multi-coil K-space. M represents the number of points actually sampled in the single-coil K-space data, where M<<N.

[0039] The SPIRiT model performs consistency between each point on the grid and its entire neighborhood across all coils, so the calibration consistency for all channels can be given by the matrix form:

[0040] x=Gx (2)

[0041] where G represents the frequency-domain self-consistent convolution operator obtained from the multi-coil K-space self-calibration region.

[0042] The image reconstruction quality can be improved by introducing the L1 norm regularization term in the SPIRiT model. Therefore, the SPIRiT parallel MRI reconstruction model with L1 norm regularization term can be expressed as the following optimization problem:

[0043]

[0044] Among them, ||·||1 represents the L1 norm, Represents coil-by-coil sparse transform, which is used to thin out the image. It is a tight framework. In this paper, SIDWT is selected as the tight framework in the experiment. represents the coil-by-coil Fourier transform, I Q represents a Q×Q identity matrix, represents the Kronecker product, represents the two-dimensional Fourier transform, Fx and F y Represents n x and n y Point Fourier transform matrix, represents the coil-by-coil inverse Fourier transform, (·) -1 Represents the inverse operation of a matrix.

[0045] Using penalty function techniques, we can obtain the unconstrained version of problem (3):

[0046]

[0047] Among them, ||·||2 represents the L2 norm, γ and λ are the parameters of the data consistency term and the L1 norm regularization term, respectively. Lustig et al. pointed out that it is important to keep the sampled data unchanged, so x is expressed as Then problem (4) can be transformed into the following form:

[0048]

[0049] in, represents the unsampled multi-coil K-space data, represents the operator for selecting unsampled points from multi-coil K-space, and denote operators for placing sampled and unsampled points back to their correct positions in the multi-coil K-space, respectively. T Represents the transpose operation of a matrix.

[0050] Using the pFISTA technique and considering the consistency of data acquisition, the solution to problem (5) can be obtained by iteratively solving the following problem:

[0051]

[0052]

[0053]

[0054]

[0055]

[0056] In formula (6) to formula (10), the superscripts of variables are k+1 "and" k " denote the variables obtained after the k+1th and kth iterations respectively. In formula (6), u k represents the intermediate variable related to the gradient problem obtained after the kth iteration, z k represents the intermediate variable obtained after the kth iteration, L represents The Lipschitz constant of the gradient, I represents a QN×QN identity matrix; in formula (7), Indicates u k The intermediate variable obtained after calibration; in formula (8), x k+1 represents the multi-coil K-space data reconstructed after the k+1th iteration, Ψ * represents the adjoint of Ψ, and specifically satisfies Ψ * Ψ=I,(·) * represents the adjoint operation of the matrix, represents the point-by-point soft threshold function; in formula (9), t k+1 and t k denote the acceleration factors of the k+1th and kth iterations respectively; in formula (10), z k+1 represents the intermediate variable obtained after the k+1th iteration, x k represents the multi-coil K-space data reconstructed after the k-th iteration.

[0057] In summary, the reconstructed multi-coil K-space data x k+1 Each coil is subjected to inverse Fourier transform, and then the square root of sum of squares (SOS) method is used to obtain the final reconstructed image. The calculation formula is as follows:

[0058]

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

[0060] S0: Initialization, let x 0 =0,z 0 =0,t 0 =0, k=0;

[0061] Among them, the superscript " 0 " represents the initial value, x 0 represents the initial value of the multi-coil K-space data x to be reconstructed; z 0 represents the initial value of the intermediate variable z; t 0 represents the initial value of the acceleration factor t; k represents the number of iterations;

[0062] S1: Calculate the intermediate variable u related to the gradient problem k , the calculation formula is as follows (6);

[0063] S2: intermediate variables u related to the gradient problem k Calibrate to get Calculation formula (7);

[0064] S3: Calculate the multi-coil K-space data x reconstructed after the k+1th iteration k+1 , the calculation formula is as follows (8);

[0065] S4: Update the intermediate variable t of the kth iteration k+1 and z k+1 , the calculation formulas are as follows (9) and (10) respectively;

[0066] S5: Determine whether the number of iterations k reaches the maximum number of iterations. If so, proceed to step S6; otherwise, set k = k + 1 and return to step S1;

[0067] S6: Reconstruct the multi-coil K-space data x k+1 Each coil is subjected to inverse Fourier transform, and then the square root of sum of squares (SOS) method is used to obtain the final reconstructed image. The calculation formula is as shown in (11).

[0068] The present invention will be further described below with reference to specific experiments.

[0069] To validate the effectiveness of the proposed fSIDWT-SPIRiT method, we compared its reconstruction performance with that of the pFISTA-SPIRiT method in the following experiments. All experiments were performed on a laptop equipped with an Intel(R) Core(TM) i5-7200U CPU @ 2.50GHz processor, 16GB of RAM, and a Windows 10 operating system (64-bit). All methods were implemented using MATLAB.

[0070] In order to compare the reconstruction performance of the two methods, the present invention selected a heart image for simulation experiments, named dataset 1 (e.g. Figure 2 ), dataset 1 is a cardiac dataset acquired using a 28-channel coil, which is then compressed into 12 virtual coils with a size of 192 × 192 using coil compression technology. To generate the test dataset, a two-dimensional Poisson disk sampling pattern with an acceleration factor of R× (excluding ACS) is used for undersampling, and all methods use a calibration area of 24 × 24 and a 5 × 5 SPIRiT kernel. Figure 3 The 2D Poisson disk subsampling mode with acceleration factors of 3× and 24×24ACS is demonstrated.

[0071] First, the present invention compares the visual presentation of dataset 1 reconstructed using the pFISTA-SPIRiT method and the fSIDWT-SPIRiT method when the acceleration factor is 3. Figure 4 and Figure 5 The reconstructed images obtained by using the pFISTA-SPIRiT method and the fSIDWT-SPIRiT method for dataset 1 are shown respectively. Figure 4 and Figure 5 It can be seen that the reconstructed image obtained using the pFISTA-SPIRiT method has more artifacts, while the reconstructed image obtained using the fSIDWT-SPIRiT method is clearer and closer to the original image.

[0072] Secondly, in order to more intuitively compare the reconstruction quality of the two methods, the present invention Figure 6 and Figure 7 The error graphs obtained by reconstructing dataset 1 using the pFISTA-SPIRiT method and the fSIDWT-SPIRiT method at 3 times acceleration are shown in Figure 1 (white dots represent errors, the more white dots, the greater the error). Figure 6 and Figure 7 It is obvious that the error between the reconstructed image obtained by the pFISTA-SPIRiT method and the original image is large, while the error between the reconstructed image obtained by the fSIDWT-SPIRiT method and the original image is smaller, indicating that the proposed method can reconstruct more details of the image and has better reconstruction quality.

[0073] Finally, this paper compares the reconstruction performance of various methods in terms of speed. This paper uses the signal-to-noise ratio (SNR) to objectively measure the quality of image reconstruction and uses time to measure the speed of image reconstruction. The higher the SNR value, the better the image reconstruction quality. The definition of SNR is as follows:

[0074]

[0075] Among them, Var represents the variance of the reference image x, and MSE represents the reconstructed image The mean squared error between the image and the reference image x.

[0076] Figure 8 The speed comparison of the reconstruction of dataset 1 using the pFISTA-SPIRiT method and the fSIDWT-SPIRiT method is shown at 3 times acceleration. Figure 8It can be seen that the SNR value of the reconstructed image obtained by the pFISTA-SPIRiT method is 22.87dB, while the SNR value of the reconstructed image obtained by the fSIDWT-SPIRiT method is 23.83dB. Therefore, the reconstructed image obtained by the fSIDWT-SPIRiT method has a higher SNR value, indicating that the reconstructed image quality obtained by the fSIDWT-SPIRiT method is better. It can also be found that it only takes 33 seconds to reconstruct dataset 1 using the fSIDWT-SPIRiT method, while it takes 133 seconds to reconstruct dataset 1 using the pFISTA-SPIRiT method, indicating that the convergence speed of the fSIDWT-SPIRiT method for reconstructing dataset 1 is significantly faster than that of the pFISTA-SPIRiT method. This shows that the method fSIDWT-SPIRiT proposed in the present invention not only achieves faster reconstruction, but also reconstructs higher quality images, further verifying the effectiveness of the proposed method.

[0077] In summary, this paper combines the L1-norm regularization term with the SPIRiT model to propose a new parallel magnetic resonance imaging method, which we call the fast parallel imaging reconstruction method based on SIDWT and iterative self-consistency (fSIDWT-SPIRiT). Experiments were then conducted on selected datasets using both the fSIDWT-SPIRiT method and the pFISTA-SPIRiT method. The experimental results show that the proposed method, fSIDWT-SPIRiT, outperforms the pFISTA-SPIRiT method in both reconstruction quality and reconstruction speed.

[0078] The above describes the specific embodiments of the present invention in detail with reference to the accompanying drawings. However, the present invention is not limited to the above embodiments. Various changes can be made within the knowledge of ordinary technicians in this field without departing from the scope of the present invention.

Claims

1. A fast parallel imaging reconstruction method based on SIDWT and iterative self-consistency, comprising the following steps: S0: Initialization, let x 0 =0,z 0 =0,t 0 =0, k=0; The superscript "0" indicates the initial value. represents the multi-coil K-space data to be reconstructed, Indicates the K-space data of the first coil, Indicates the Qth coil K-space data, Q represents the number of receiving coils, represents the qth coil K-space data, q=1,...,Q represents the coil index variable, (·) H Represents the conjugate transpose operation of a vector or matrix, N = n y ×n x Indicates the number of pixels in a single image, n y and n x Represents the number of rows and columns of a single image, x 0 represents the initial value of x; represents the intermediate variable, z 0 represents the initial value of z; represents the factor related to acceleration, t 0 represents the initial value of t; k represents the number of iterations; S1: Calculate the intermediate variable u related to the gradient problem k , the calculation formula is as follows: Among them, z k represents the intermediate variable obtained after the kth iteration, L represents The Lipschitz constant of the gradient, G represents the frequency domain self-consistent convolution operator obtained from the multi-coil K-space self-calibration region, and I represents a QN×QN identity matrix; S2: intermediate variables u related to the gradient problem k Calibrate to get The calculation formula is as follows: in, and denote operators for selecting sampled and unsampled points from the multi-coil K-space, and denote operators for placing sampled and unsampled points back to their correct positions in the multi-coil K-space, respectively. T represents the transpose operation of the matrix, represents undersampled multi-coil K-space data, M represents the number of points actually sampled in single-coil K-space data, M<<N; S3: Calculate the multi-coil K-space data x reconstructed after the k+1th iteration k+1 , the calculation formula is as follows: in, represents the coil-by-coil Fourier transform, I Q represents a Q×Q identity matrix, represents the Kronecker product, represents the two-dimensional Fourier transform, F x and F y Represents n x and n y Point Fourier transform matrix, Represents coil-by-coil sparse transform, which is used to thin out the image. It is a tight framework. SIDWT is selected as the tight framework in the experiment. * represents the adjoint of Ψ, and specifically satisfies Ψ * Ψ=I,(·) * represents the adjoint operation of the matrix, represents the point-by-point soft threshold function, λ>0 represents the regularization parameter, represents the coil-by-coil inverse Fourier transform, (·) -1 Represents the inverse operation of the matrix; S4: Update the intermediate variable t of the kth iteration k+1 and z k+1 , the calculation formulas are as follows: Among them, z k+1 represents the intermediate variable obtained after the k+1th iteration, t k+1 and t k They represent the acceleration factors for the k+1th and kth iterations, respectively, and x k represents the multi-coil K-space data reconstructed after the k-th iteration; S5: Determine whether the number of iterations k reaches the maximum number of iterations. If so, proceed to step S6; otherwise, set k = k + 1 and return to step S1; S6: Reconstruct the multi-coil K-space data x k+1 Each coil is inverse Fourier transformed, and then the square root of the sum of squares (SOS) method is used to obtain the final reconstructed image. The calculation formula is as follows:

Images