A sensitivity coding reconstruction method based on dictionary learning of classified image patches

Abstract

The present invention relates to a sensitivity coding reconstruction method based on classified image block dictionary learning, and belongs to the field of magnetic resonance imaging. Introducing a fast dictionary learning on classified patches (FDLCP) method in magnetic resonance image reconstruction can improve sparsity. Based on a sensitivity coding (SENSE) reconstruction framework, the present invention proposes a magnetic resonance imaging reconstruction model based on classified image block dictionary learning, named FDLCP‑SENSE, which is solved using the alternating direction method of multipliers (ADMM). Experimental comparison with the TV‑SENSE algorithm shows that FDLCP‑SENSE can obtain clearer reconstructed images with smaller reconstruction errors, and is an effective magnetic resonance imaging reconstruction method.

Description

Technical Field

[0001] The invention relates to a sensitivity coding reconstruction method based on classified image block dictionary learning, and belongs to the technical field of magnetic resonance imaging. Background Art

[0002] Magnetic resonance imaging (MRI), an important clinical diagnostic imaging tool, provides clear images of human tissue without ionizing radiation damage. It acquires image data in K-space through frequency and phase encoding. However, its imaging speed is primarily limited by the lengthy relaxation time between phase encoding steps, resulting in a relatively slow imaging speed. Parallel imaging (PI) is an important technology for accelerating MRI.

[0003] Parallel imaging uses multi-channel phased array coils to simultaneously acquire magnetic resonance signals and encodes spatial information by exploiting the differences in the spatial sensitivity of each coil. Sensitivity encoding (SENSE) is a widely used parallel MRI technique that explicitly utilizes coil sensitivity information to recover image information.

[0004] Image reconstruction based on coil sensitivity estimation using the SENSE model requires regularization to suppress noise and reduce aliasing. Common SENSE model-based regularization methods include the TV regularization term combined with the SENSE model (TV-SENSE algorithm) and the Lp pseudo-norm total variation regularization term combined with the SENSE model (LpTV-SENSE algorithm). These algorithms enhance image sparsity and effectively preserve edge and contour information, but they still exhibit artifacts. Summary of the Invention

[0005] The purpose of the present invention is to overcome the shortcomings of the prior art and provide a sensitivity coding reconstruction method based on classified image block dictionary learning, which can further improve the quality of reconstructed images in magnetic resonance imaging.

[0006] The technical solution adopted by the present invention is: a sensitivity coding reconstruction method based on classified image block dictionary learning, comprising the following steps:

[0007] S0: Initialization, let x 0 =(RFS) H y,z 0 =0,α 0 =0,

[0008] Among them, the superscript " 0 " indicates the initial value before iteration, and the superscript " H" represents the conjugate transpose operation, Represents the image to be reconstructed, N=N v ×N h , N v and N h are the number of rows and columns of the image, respectively, x 0 represents the initial value of x; It is an undersampling matrix, M represents the number of points actually sampled by the single coil K-space data, and M<<N, I L represents the L×L unit matrix, L represents the number of receiving coils used for parallel imaging, represents the matrix for selecting sampling points from the K-space grid of a single coil, R T represents the transpose of R, represents undersampled multi-coil K-space data, represents the two-dimensional Fourier transform, U h 、U v Represents N h 、N v The Fourier transform matrix of the point, represents the Kronecker product, represents the sensitivity matrix, Indicates the sensitivity of the first coil, The sensitivity of the Lth coil is expressed as Indicates that S l is a diagonal matrix.

[0009] represents auxiliary variables, n represents the number of pixels in the image block, z 0 and α 0 represents the initial values of z and α, and denote the Lagrange multipliers corresponding to z and α, respectively, and Indicates u z and u α The initial value of

[0010] S1: reconstruct image x from the initial 0 The matrix for extracting the j-th image block containing n pixels is: The extracted j-th image block is expressed as: G j x 0 , there are N overlapping image patches, j = 1...N. According to the fast dictionary learning on classified patches (FDLCP) method, the N image patches are divided into geometric directions ω∈{ω1...ωQ}, and learn the adaptive analysis dictionary with geometric direction ω make For the adaptive analysis dictionary with the geometric direction ω of the j-th image block, a transformation matrix Φ containing multiple information can be constructed: Then there is

[0011] S2: Initialization, number of iterations k = 0;

[0012] S3: Calculate the auxiliary variable z for the k+1th iteration k+1 , the calculation formula is as follows:

[0013]

[0014] in, represents the auxiliary variable for the k+1th iteration, (·) -1 Indicates the inversion operator, superscript " T " represents the transpose operation of the matrix, represents the image to be reconstructed at the kth iteration, Indicates the number corresponding to z in the kth iteration k The Lagrange multiplier of I NL represents the NL×NL identity matrix, μ1>0 represents the parameter;

[0015] S4: Calculate the auxiliary variable α for the k+1th iteration using the hard threshold method k+1 , the calculation formula is as follows:

[0016]

[0017] in, Indicates the kth iteration corresponding to α k Lagrange multiplier, λ>0 and μ2>0 represent parameters, H(β,θ) represents the hard threshold function, and the calculation formula is: β represents the input matrix, θ represents the threshold;

[0018] S5: Calculate the image x to be reconstructed at the k+1th iteration k+1 , the calculation formula is as follows:

[0019]

[0020] Among them, S H represents the conjugate transpose of the sensitivity matrix S, F H represents the conjugate transpose of the Fourier transform matrix F, Φ H represents the conjugate transpose of the transformation matrix Φ, I N represents the N×N identity matrix;

[0021] S6: Update the Lagrange multiplier of the k+1th iteration Among them, δ represents the step size of 1;

[0022] S7: Update the Lagrange multiplier of the k+1th iteration

[0023] S8: Calculate x k+1 and x k The relative error (RE) between them is calculated as follows:

[0024]

[0025] S9: Determine whether the stopping condition is met. If RE is less than the tolerance tol, or greater than the maximum number of iterations K1, that is, k>K1, proceed to step S10; otherwise, set k=k+1 and return to S2;

[0026] S10: output x=x k+1 , and obtain the reconstructed image x.

[0027] The present invention has the following beneficial effects: Sensitivity coding is a method for image reconstruction that utilizes the sensitivity information of multiple coils. Based on the SENSE model, the present invention proposes a sensitivity coding reconstruction method based on dictionary learning of classified image blocks. The present invention uses the alternating direction multiplier method to transform the reconstruction problem into multiple subproblems for iterative solution. Experimental results show that compared with the TV-SENSE algorithm, the proposed method has significant advantages in signal-to-noise ratio, structural similarity index, and high-frequency error norm, resulting in better image reconstruction quality. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] Figure 1 Flow chart of the method of the present invention;

[0029] Figure 2 This is human brain data acquired using an 8-channel coil;

[0030] Figure 3 It is a 2D Poisson disk undersampling mask with a 4x acceleration factor and a 24×24 center fully sampled self-calibration area;

[0031] Figures 4-5 Images reconstructed from dataset 1 using the TV-SENSE and FDLCP-SENSE methods, respectively, which are undersampled from a two-dimensional Poisson disk undersampling mask with a 4x acceleration factor and a 24×24 central fully sampled self-calibration area;

[0032] Figures 6-7 They are Figures 4-5The error map corresponding to the reconstructed image is the error map reconstructed by the TV-SENSE algorithm and the FDLCP-SENSE method. DETAILED DESCRIPTION

[0033] The technical solution of the present invention is further described in detail below with reference to the accompanying drawings;

[0034] Example 1: The present invention is an efficient reconstruction method proposed based on the SENSE framework, comprising the following steps:

[0035] 1) SENSE framework:

[0036] SENSE is a parallel imaging method, and its reconstruction model can be expressed as:

[0037] y=RFSx (1)

[0038] Where: is the image to be reconstructed, N=N v ×N h , N v and N h are the number of rows and columns of the image, respectively; represents undersampled multi-coil K-space data, M represents the number of points actually sampled by single-coil K-space data, M<<N, and L represents the number of receiving coils used for parallel imaging. is the K-space undersampling matrix, I L represents the L×L identity matrix, represents the matrix for selecting sampling points from the K-space grid of a single coil, R T Represents the transpose of the matrix R. represents the two-dimensional Fourier transform, U h 、U v Represents N h 、N v The Fourier transform matrix of the point, represents the Kronecker product; represents the sensitivity matrix, Indicates the sensitivity of the first coil, The sensitivity of the Lth coil is expressed as Indicates that S l is a diagonal matrix, the superscript " H ” indicates the conjugate transpose operation.

[0039] Since the regularization term is effective in reducing aliasing artifacts and noise effects in reconstruction, it can effectively improve the sparsity of data. Therefore, the reconstruction method combining the regularization term with the SENSE model is expressed as:

[0040]

[0041] Where: Ψ(x) is the regularization term. Common regularization terms are: Regularization term, total variation (TV) regularization term, etc.

[0042] 2) Algorithm derivation:

[0043] To further improve the quality of MR image reconstruction, the fast dictionary learning on classified patches (FDLCP) method is combined with the SENSE model of parallel MRI reconstruction to obtain the following optimization problem:

[0044]

[0045] Where: λ is the regularization parameter, j=1...N, ω∈{ω1,ω2,...,ω Q} represents the geometric direction, Q is the number of categories of image blocks classified according to the geometric direction; is the adaptive analysis dictionary with the geometric direction ω of the j-th image block; is the matrix for extracting the jth image block containing n pixels from the reconstructed image x, where n represents the number of pixels in the image block.

[0046] According to the FDLCP method, N image blocks are divided into geometric directions ω∈{ω1...ω Q}, and learn the adaptive analysis dictionary with geometric direction ω make For the adaptive analysis dictionary with the geometric direction ω of the j-th image block, the transformation matrix Φ can be constructed:

[0047]

[0048] So we have: Φ H Φ=I.

[0049] Using the above formula, model (3) can be reformulated as:

[0050]

[0051] Introducing auxiliary variables and the corresponding Lagrange multipliers and Using ADMM technology, problem (5) can be transformed into iterative solutions to the following sub-problems:

[0052]

[0053]

[0054]

[0055]

[0056]

[0057] In equations (6)-(10), the superscripts “k+1” and “k” of the variables represent the variables of the k+1th iteration and the kth iteration, respectively; in equation (6), represents auxiliary variables, represents the auxiliary variable of the k+1th iteration, μ1>0, represents the image to be reconstructed at the kth iteration, Indicates the auxiliary variable z in the kth iteration k Lagrange multiplier of; In formula (7), represents auxiliary variables, α k+1 represents the auxiliary variable of the k+1th iteration, μ2>0, Indicates the auxiliary variable α in the kth iteration k Lagrange multiplier of; In formula (8), x k+1 represents the image to be reconstructed at the k+1th iteration; in equations (9) and (10), is the auxiliary variable z in the k+1th iteration k+1 The Lagrange multiplier of is the auxiliary variable α in the k+1th iteration k+1 The Lagrange multiplier of , δ represents the step size of 1.

[0058] Formulas (9)-(10) are for u z and u α The update method is then used to calculate the other formulas mentioned above.

[0059] Regarding the sub-problem of z, according to formula (6), the solution of z is:

[0060]

[0061] in,(·) -1 Represents the inversion operator, I NL Represents the NL×NL identity matrix.

[0062] Regarding the sub-problem of α, according to formula (7), using the hard threshold method, the solution of α can be obtained as:

[0063]

[0064] The hard threshold operator H(β,θ) is defined as follows:

[0065]

[0066] In formula (13), β represents the input vector and θ represents the threshold.

[0067] Regarding the subproblem of x, according to formula (8), by keeping other variables unchanged and setting the derivative of the objective function with respect to x to 0, we can obtain:

[0068]

[0069] Among them, Φ H represents the transpose of the transformation matrix Φ, S H F H FS=S H S is a diagonal matrix, μ1S H The inverse of S+μ2I is easy to obtain, so the analytical solution of formula (8) is:

[0070]

[0071] Among them, S H represents the conjugate transpose of the sensitivity matrix S, I N Represents the N×N identity matrix.

[0072] Calculate x k+1 and x k The relative error (RE) between them is calculated as follows:

[0073]

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

[0075] S0: Initialization, let x 0 =(RFS) H y,z 0 =0,α 0 =0,

[0076] Among them, the superscript " 0 " indicates the initial value before iteration, and the superscript " H " represents the conjugate transpose operation, Represents the image to be reconstructed, N=N v ×N h , N v and N h are the number of rows and columns of the image, respectively, x 0 represents the initial value of x; is the K-space undersampling matrix, M represents the number of points actually sampled by the single-coil K-space data, and M<<N, I L represents the L×L unit matrix, L represents the number of receiving coils used for parallel imaging, represents the matrix for selecting sampling points from the K-space grid of a single coil, R T represents the transpose of R, represents undersampled multi-coil K-space data, represents the two-dimensional Fourier transform, U h 、U v Represents N h 、N v The Fourier transform matrix of the point, represents the Kronecker product, represents the sensitivity matrix, Indicates the sensitivity of the first coil, The sensitivity of the Lth coil is expressed as Indicates that S l is a diagonal matrix. represents auxiliary variables, n represents the number of pixels in the image block, z 0 and α 0 represents the initial values of z and α, and denote the Lagrange multipliers corresponding to z and α, respectively, and Indicates u z and u α The initial value of .

[0077] S1: reconstruct image x from the initial 0 The matrix for extracting the j-th image block containing n pixels is: The extracted j-th image block is expressed as: G j x 0 , there are N overlapping image blocks in total, j = 1...N. According to the FDLCP method, the N image blocks are divided into geometric directions ω∈{ω1...ω Q}, and learn the adaptive analysis dictionary with geometric direction ω make For the adaptive analysis dictionary with the geometric direction ω of the j-th image block, the transformation Φ can be constructed: Then there is

[0078] S2: Initialization, number of iterations k = 0;

[0079] S3: Calculate the auxiliary variable z for the k+1th iteration k+1, the calculation formula is as follows (11);

[0080] S4: Calculate the auxiliary variable α for the k+1th iteration using hard thresholding k+1 , the calculation formula is as follows (12);

[0081] S5: Calculate the image x to be reconstructed at the k+1th iteration k+1 , the calculation formula is as follows (15);

[0082] S6: Update the Lagrange multiplier of the k+1th iteration Formula (9)

[0083] S7: Update the Lagrange multiplier of the k+1th iteration Formula is as follows (10);

[0084] S8: Calculate x k+1 and x k The relative error (RE) between them is calculated as follows (16);

[0085] S9: Determine whether the stopping condition is met. If RE < tol or k is greater than the maximum number of iterations K1, that is, k > K1, proceed to step S10; otherwise, set k = k + 1 and return to S1;

[0086] S10: output x=x k+1 , and obtain the reconstructed image x.

[0087] Experimental results:

[0088] In the following experiments, to validate the performance of our proposed method, we compared it with the TV-SENSE algorithm, which combines the TV regularization term with the SENSE model. All algorithms were implemented in MATLAB. All experiments were performed on a laptop equipped with an i7 7500U @ 2.70GHz CPU, 8GB of RAM, and the Windows 10 operating system.

[0089] In order to compare the performance of each algorithm, the present invention uses a human brain slice image to conduct simulation experiments, named dataset 1. Figure 2 As shown in Figure 2, a 2D Poisson disk undersampling mask with an acceleration factor of AF = 4 and a 24 × 24 center fully sampled self-calibration area is used, as shown in Figure 2. Figure 3 shown.

[0090] A visual comparison is performed on the two reconstruction algorithms under the dataset dataset1. The present invention selects the reconstructed images when the acceleration factor is 4 for comparison. Figures 4-5 The reconstructed images of the TV-SENSE algorithm and the FDLCP-SENSE algorithm are shown respectively. Figure 4This is the image reconstructed using the TV-SENSE algorithm. It can be seen that the edge contours of the image are not clear and some details are not reconstructed. Figure 5 This is the image reconstructed using FDLCP-SENSE. The reconstruction quality of this algorithm is better than that of the TV-SENSE algorithm, and the reconstructed image is highly consistent with the original image.

[0091] In order to further illustrate that the new algorithm FDLCP-SENSE proposed in this invention is superior to the TV-SENSE algorithm in terms of reconstruction performance, Figures 6-7 The error plots for dataset 1 using the two algorithms are shown (whiter error plots indicate greater error). The error plots show that TV-SENSE has more granular artifacts and a larger error, while the FDLCP-SENSE algorithm has a smaller error, indicating better image reconstruction quality.

[0092] In summary, experiments were conducted on different reconstruction algorithms for the selected dataset dataset1. The FDLCP-SENSE algorithm showed its advantages in both evaluation indicators and visual effects, and its visual effect was significantly better than that of the TV-SENSE algorithm.

[0093] The above is a detailed analysis and description of the specific embodiments of the present invention in conjunction with the accompanying drawings, but the present invention is not limited to the above embodiments. Various changes can be made within the knowledge scope of ordinary technicians in this field without departing from the purpose of the present invention.

Claims

1. A sensitivity coding reconstruction method based on dictionary learning of classified image blocks, characterized by: The following steps are involved: S0: Initialization, let x 0 =(RFS) H y,z 0 =0,α 0 =0, Among them, the superscript " 0 " indicates the initial value before iteration, and the superscript " H " represents the conjugate transpose operation, Represents the image to be reconstructed, N=N v ×N h , N v and N h are the number of rows and columns of the image, respectively, x 0 represents the initial value of x; is the K-space undersampling matrix, M represents the number of points actually sampled by the single-coil K-space data, and M<<N, I L represents the L×L unit matrix, L represents the number of receiving coils used for parallel imaging, represents the matrix for selecting sampling points from the K-space grid of a single coil, R T represents the transpose of R, represents undersampled multi-coil K-space data, represents the two-dimensional Fourier transform, U h 、U v Represents N h 、N v The Fourier transform matrix of the point, represents the Kronecker product, represents the sensitivity matrix, Indicates the sensitivity of the first coil, The sensitivity of the Lth coil is expressed as Indicates that S l is a diagonal matrix; represents auxiliary variables, n represents the number of pixels in the image block, z 0 and α 0 represents the initial values of z and α, and denote the Lagrange multipliers corresponding to z and α, respectively, and Indicates u z and u α The initial value of S1: reconstruct image x from the initial 0 The matrix for extracting the j-th image block containing n pixels is: The extracted j-th image block is expressed as: G j x 0 , there are N overlapping image blocks in total, j = 1...N, according to the fast dictionary learning FDLCP method based on classified image blocks, the N image blocks are divided into geometric directions ω∈{ω1...ω Q }, and learn the adaptive analysis dictionary with geometric direction ω make For the adaptive analysis dictionary with the geometric direction ω of the j-th image block, a transformation matrix Φ containing multiple information can be constructed: Then there is S2: Initialization, number of iterations k = 0; S3: Calculate the auxiliary variable z for the k+1th iteration k+1 , the calculation formula is as follows: in, represents the auxiliary variable for the k+1th iteration, (·) -1 Represents the inverse operator, superscript " T " represents the transpose operation of the matrix, represents the image to be reconstructed at the kth iteration, Indicates the number corresponding to z in the kth iteration k The Lagrange multiplier of I NL represents the NL×NL identity matrix, μ1>0 represents the parameter; S4: Calculate the auxiliary variable α for the k+1th iteration using the hard threshold method k+1 , the calculation formula is as follows: in, Indicates the k-th iteration corresponding to α k Lagrange multiplier, λ>0 and μ2>0 represent parameters, H(β,θ) represents the hard threshold function, and the calculation formula is: β represents the input matrix, θ represents the threshold; S5: Calculate the image x to be reconstructed at the k+1th iteration k+1 , the calculation formula is as follows: Among them, S H represents the conjugate transpose of the sensitivity matrix S, F H represents the conjugate transpose of the Fourier transform matrix F, Φ H represents the conjugate transpose of the transformation matrix Φ, I N represents the N×N identity matrix; S6: Update the Lagrange multiplier of the k+1th iteration Among them, δ represents the step size of 1; S7: Update the Lagrange multiplier of the k+1th iteration S8: Calculate x k+1 and x k The relative error (RE) between them is calculated as follows: S9: Determine whether the stopping condition is met. If RE is less than the tolerance tol, or greater than the maximum number of iterations K1, that is, k>K1, proceed to step S10; otherwise, set k=k+1 and return to S2; S10: output x=x k+1 , and obtain the reconstructed image x.

Images