CN113628298B - Feature vector based self-consistency and non-local low-rank parallel MRI reconstruction method - Google Patents

Abstract

The invention relates to a self-consistency and non-local low-rank parallel MRI reconstruction method based on feature vectors, and belongs to the technical field of magnetic resonance imaging. The parallel imaging reconstruction model based on the iterative self-consistency of the feature vectors estimates a plurality of groups of sensitivity maps through the feature value decomposition of the k-space automatic calibration data, and performs reconstruction by using a method similar to sensitivity coding. The invention provides a parallel magnetic resonance imaging reconstruction model based on non-local low-rank constraint of images, which is named NLR-ESPIRiT, based on an ESPIRiT reconstruction framework, and is effectively solved by using a rapid iteration shrinkage threshold technology. Experiments compare the L1-ESPIRiT model with the L1-norm regularization term and the LpJTV-ESPIRiT model with the Lp pseudo-norm joint total variation regularization term. Experimental results show that the NLR-ESPIRiT algorithm can better remove undersampling artifacts.

Description

Feature vector based self-consistency and non-local low-rank parallel MRI reconstruction method

Technical Field

The invention relates to a self-consistency and non-local low-rank parallel MRI reconstruction method based on feature vectors, and belongs to the technical field of magnetic resonance imaging.

Background

Magnetic resonance imaging (Magnetic Resonance Imaging, MRI) is a technique for medical imaging using the principles of nuclear magnetic resonance, and is currently a few safe and accurate clinical diagnostic method without any harm to the human body. Although MRI has been widely used in clinical medicine, MRI technology has long scanning time and slow imaging, and not only can image artifacts be generated due to organ motion, but also the requirement of dynamic real-time imaging cannot be met, which greatly limits the development of MRI technology.

Parallel imaging is a common method for accelerating MRI, which uses a multi-channel phased array to acquire magnetic resonance signals simultaneously, and uses the spatial sensitivity difference of each coil to perform sensitivity encoding of spatial information. One type of parallel imaging technique knows sensitivity information, such as sensitivity encoding (SENSitivity Encoding, SENSE) based on coil sensitivity information, when accurate spatial sensitivity information is known, SENSE can achieve optimal reconstruction results. However, in practical applications, it is difficult to measure the sensitivity information with high accuracy, and a small sensitivity error may cause a deviation in the reconstruction result. Another type of automatic calibration method based on k-space kernel, such as a general self-calibration partial parallel sampling (GeneRalized Autocalibarating Partially Parallel Acquisition, GRAPPA) and Iterative self-consistency parallel imaging reconstruction (spiit) method, can complete accurate reconstruction only by using a small number of central full sampling signals (ACS) for calibration.

SENSE, however, is constrained by the difficulty of obtaining accurate sensitivity information; the SPIRiT method requires coil-by-coil reconstruction, which takes a longer time. In combination with the sensitivity information-based and k-space-based method of automatically calibrating kernels, uecker et al propose an iTerative Self-consistency parallel imaging reconstruction (ESPIRiT) model based on feature vectors. First, the k-space auto-calibration data is applied for eigenvalue decomposition, and the principal eigenvector may represent the coil sensitivity information. Data reconstruction is then performed using a SENSE-like method based on the estimated sensitivity information. Unlike SENSE, the ESPIRiT method can estimate coil sensitivity through k-space calibration information without directly using the sensitivity information; compared with SPIRIT, the ESPIRiT method only needs to reconstruct images according to two groups of sensitivity maps, and can reconstruct parallel magnetic resonance images more effectively.

The regularization technology based on the image priori information is a common method for recovering image data, and the performance of a reconstruction algorithm can be effectively improved by adding regularization terms into an algorithm reconstruction model. However, the existing improved method based on the ESPIRiT only comprises an L1-ESPIRiT algorithm combined with a simple L1 norm regularization term and an LPJTV-ESPIRiT algorithm combined with an Lp pseudo-norm joint total variation (LpJTV) regularization term, and the quality of an algorithm reconstructed image still needs to be improved.

Disclosure of Invention

The invention aims to overcome the defects of the prior art and provide a self-consistency and non-local low-rank parallel MRI reconstruction method based on feature vectors, which can further improve the quality of a magnetic resonance imaging reconstructed image.

The aim of the invention is realized by the following technical scheme: a parallel MRI reconstruction method based on self-consistency of eigenvectors and non-local low rank. The method is characterized by comprising the following steps of:

s0: initializing, lettingz ⁰ ＝0,D ⁰ ＝0,d ⁰ ＝0,t ⁰ ＝0，k＝0；

Wherein the superscript is " ⁰ "represents the initial value before iteration, superscript" ^T "means the conjugate transpose operator.Representing a column-vectorized multicomponent image to be reconstructed, n=n _h ×N _v Pixel count, N, representing a single coil amplitude image _v and N_h The number of lines and columns of the single coil amplitude image, J represents the number of sensitivity map sets, and the number of image components of x is equal to the number of sensitivity map sets. x for the jth component image of x _j Expressed, j=1..j, expressed an index of the sensitivity map group, for example, x ₁ and x_J Respectively representing the 1 st and the J-th component images of the column-vectorized multicomponent image x to be reconstructed. X is x ⁰ The initial value of x is indicated. />Representing column vectorized undersampled multi-coil k-space data, M representing the number of single coil k-space undersampled data points, and C representing the number of receive coils used for parallel imaging. y for the c-th coil data of y _c And c=1,.. representing the index of the coil, e.g. y ₁ and y_C Respectively representing the 1 st and C-th coil data of the column-vectorized undersampled multi-coil k-space data y. />Andrepresenting Fourier operator and under-run, respectivelySampling operators. I _C Is a unit matrix of C x C, and />Respectively represent N _h and N_v Point Fourier transform matrix, < >>Matrix representing the selection of sampling point positions from a k-space grid,/for>Representing the kronecker product. />For the sensitivity map group matrix, the j-th sensitivity map of the c-th coil of S is obtained by using the SPIRIT model based on the feature vector _cj Representing, for example: s is S ₁₁ Group 1 sensitivity maps representing the 1 st coil, S _CJ A group J sensitivity map of the C-th coil is shown. />And (3) withRepresenting intermediate variables, z ⁰ ，t ⁰ And->Are their corresponding initial values. Defining block matching operatorsThe operator performs the block matching on x as follows: will->Divided into overlapping N _p The size is +.>N represents the number of pixels of the image block, m is the number of similar blocks, i=1,.. _p Is a block index. Let the reference picture block be +.>In a search window with the size of s x s, selecting m image blocks with the smallest Euclidean distance with a reference block, vectorizing each image block column and forming a similar block group->Is a column of (c).Similar block group V representing the ith reference image block of the jth component image _ji (x) Corresponding auxiliary variables, all { D } _ji Sequentially and horizontally spliced into a matrix +.>D ⁰ Is the initial value of D. />Is corresponding to D _ji Lagrangian multipliers of (2), all { d }, will be _ji Sequentially spliced into matrix ++>d ⁰ Is the initial value of d.

S1, calculating the intermediate variable z of the (k+1) th iteration ^k+1 The calculation formula is as follows:

wherein ,intermediate variable representing the kth iteration +.>L represents Lipschitz constant.

S2, initializing, j=1;

s3, initializing, i=1;

s4, reconstructing an image x of the kth iteration ^k Performing similar block matching to obtain a similar block group V _ji (x ^k )；

S5, pairingSingular value decomposition to obtain ++>

Wherein, sigma is a singular value diagonal matrix,both Γ and gamma are unitary matrices, the columns of which are left and right singular vectors of singular values, respectively,/->Is expressed as the first singular value of (2) and />Respectively indicate->I=1,..min (m, n), which is a singular value index, +.>Representing the corresponding auxiliary variable in the kth iterationLagrangian multiplier of (c).

S6, initializing estimation by the following formula

S7, calculating weight w= [ w ] ₁ ,...,w _min(m,n) ]，w ₁ And w is equal to _min(m,n) The 1 st and last value of w, respectively. The first value of w is w _l The calculation formula is as follows:

wherein b₀ Epsilon=10 as a parameter ^-16 For ensuring that the denominator is not 0.

S8, calculating auxiliary variables of the (k+1) th iterationThe calculation formula is as follows:

wherein: soft (Σ, w) =max (Σ -w, 0) is a soft threshold operator.

S9, judging whether N in each coil is completed _p Personal (S)When i=n _p Step S10 is entered at the moment; otherwise, let i=i+1 and return to S4;

s10, judging whether all image components are finished, and if j=J, proceeding to step S11; otherwise, let j=j+1 and return to S3;

s11, calculating an image x to be reconstructed of the (k+1) th iteration ^k+1 The calculation formula is as follows:

wherein the superscript is " ^* "is an accompanying operator, superscript" ^-1 "means that the inversion operator is used to calculate the inversion,is V _ji The companion matrix representing ∈>And replacing the corresponding position in the original multi-component image.

Is a diagonal matrix, each diagonal element of which is equal to the corresponding pixel point appearing in all the groups of similar blocks V _ji (x ^k ) Is a number of times in the past.

S12, updating d of the (k+1) th iteration ^k+1 ，d ^k+1 Is a sub-matrix of (2)The calculation can be performed by the following formula:

wherein η > 0 represents a parameter for updating the Lagrangian multiplier;

s13, reconstructed multicomponent image x ^k+1 The reconstructed single coil amplitude image is calculated using the square root of the sum of squares (square root of sum of squares, SOS) as follows:

wherein The single coil amplitude reconstruction image representing the k+1th iteration.

S14: updating intermediate variables for the k+1st iteration

S15: calculating intermediate variables for the k+1th iteration

S16, calculating X ^k+1 and X^k Relative Error (RE) between (a) and (b) is calculated as:

wherein The single coil amplitude reconstruction image representing the kth iteration.

S17: judging whether an algorithm stopping standard is met, if so, entering step S18; otherwise, setting k=k+1, and returning to S1;

s18: outputting the final reconstructed single-coil amplitude image X ^k+1 。

The beneficial effects of the invention are as follows: the ESPIRiT model can estimate sensitivity information by eigenvalue decomposition of k-space auto-calibration data and reconstruct the data using a method similar to sensitivity encoding. Based on an ESPIRiT model, the invention provides a high-quality parallel imaging reconstruction algorithm adopting an image non-local low rank as a constraint term. First, the present invention uses a fast iterative shrinkage threshold algorithm (Fast iterative shrinkage/thresholding algorithm, FISTA) to translate the reconstruction problem into a gradient problem, a denoising sub-problem with respect to rank and two acceleration steps. Then, for the non-local low-rank denoising problem therein, a multiplier-alternate direction method (alternating direction method of multipliers, ADMM) is used, and a weighted kernel norm is used as a rank agent for solving. Theoretical analysis and imaging experiments show that compared with other ESPIRiT-based algorithms (such as L1-ESPIRiT and LpJTV-ESPIRiT), the novel parallel imaging algorithm provided by the invention can effectively improve the quality of the reconstructed image.

Drawings

FIG. 1 is a flow chart of the method of the present invention;

FIG. 2 is a diagram of a brain image (i.e., dataset 1) acquired with an eight-channel receive coil for full sampling with a Field of view (FOV) smaller than the subject;

FIG. 3 is a two-dimensional poisson disk undersampling mask map with a 5-fold acceleration factor and a 24 x 24 center full-sampling self-correcting region;

FIGS. 4-6 show images reconstructed from a dataset1 dataset undersampled using NLR-ESPIRIT algorithm, L1-ESPIRIT algorithm, and LpJTV-ESPIRIT algorithm, respectively, from a two-dimensional poisson disk undersampling mask containing 5-fold acceleration factors and 24X 24 center full-sampling self-correction regions;

fig. 7-9 show the error images corresponding to the reconstructed images of fig. 4-6, i.e., the error images reconstructed by the NLR-ESPIRiT algorithm, the L1-ESPIRiT algorithm, and the LPJTV-ESPIRiT algorithm, respectively.

Detailed Description

The technical scheme of the invention is further described in detail below with reference to the attached drawings and specific embodiments.

Example 1: the invention provides a high-efficiency reconstruction method based on an ESPIRiT framework.

1) ESPIRiT framework:

assume thatRepresenting a column-vectorized multicomponent image to be reconstructed, n=n _h ×N _v Pixel count, N, representing a single coil amplitude image _v and N_h The number of rows and columns of single coil amplitude images, J represents the number of coilsThe number of sensitivity map sets, x, is equal to the number of sensitivity map sets. x for the jth component image of x _j Expressed, j=1..j, expressed an index of the sensitivity map group, for example, x ₁ and x_J Respectively representing the 1 st and the J-th component images of the column-vectorized multicomponent image x to be reconstructed. X is x ⁰ The initial value of x is indicated. />Representing column vectorized undersampled multi-coil k-space data, M representing the number of single coil k-space undersampled data points, and C representing the number of receive coils used for parallel imaging. y for the c-th coil data of y _c And c=1,.. representing the index of the coil, e.g. y ₁ and y_C Respectively representing the 1 st and C-th coil data of the column-vectorized undersampled multi-coil k-space data y. Then, the ESPIRiT reconstruction model is obtained as follows:

wherein , and />Respectively representing fourier operators and undersampled operators. />Representing fourier operator +_> and />Respectively represent N _h and N_v Point two-dimensional Fourier transform matrix,>represents the Cronecker product, I _C Is a matrix of units of CxC, < >>Representing a matrix of sampling point locations selected from a k-space grid. />Representing the sensitivity operator in ESPIRiT:

wherein the j-th group sensitivity map of the c-th coil of S is S _cj Representing, that is: s is S ₁₁ Group 1 sensitivity maps representing the 1 st coil, S _CJ A group J sensitivity map of the C-th coil is shown.

In ESPIRiT, the calculation formula for the sensitivity operator is:

where r denotes the k-space position index,the sensitivity vector at the r-position is represented, representing a convolution operator with a matrix value kernel, defined as:

wherein R_r Representing selection of k-space r-position neighborsOperators of k-space blocks, superscripts' ^-1” and “^H "represents an inversion operation and a conjugate transpose operation, respectively. V (V) _|| From the k-space calibration matrix a: singular value decomposition of a is performed to obtain a=uΛv ^H Where Λ is a singular value diagonal matrix of a, U and V are unitary matrices, and their columns are left and right singular vectors of singular values, respectively. Decomposing V to obtain zero space span V of A _⊥ And row space span V of A _|| . Since the columns of V are the basis of the rows of matrix A, all self-calibration information is located at V _|| Is a kind of medium.

From (3), it can be found that the sensitivity map s (r) is equivalent toThe feature vector corresponding to the feature value 1 of (a) can be obtained by the method of (a) for +.>And decomposing the characteristic value to obtain the product. When not finding +>When the eigenvalue is equal to 1, s (r) is set to zero. When only one set of sensitivity maps is used, the reconstruction problem can be expressed as standard SENSE. The ESPIRiT model uses multiple sets of sensitivity maps for parallel MRI reconstruction, as small FOVs can cause aliasing artifacts. Therefore, multiple sets of sensitivity maps need to be considered:

wherein λ_j (r) is a parameter at the r position, typically close to 1.s is(s) _cj (r) is s _j The j-th component of (r). Diagonal matrix S _cj S through all positions _j (r) the rearrangement is carried out to obtain,

2) Algorithm derivation:

although the ESPIRiT algorithm incorporating the regularization method of LpJTV can effectively improve the reconstruction quality, there is still much room for improvement. In order to further improve reconstruction quality, the invention combines non-local low-rank regularization (NLR) with an ESPIRIT model, and proposes an NLR-ESPIRIT algorithm to obtain the following optimization problems:

wherein rank (V) _ji (x) For V) _ji (x) Regular terms of low rank constraint are performed. Wherein the method comprises the steps ofIs a block matching operator, and the step of using the operator to perform block matching on x is as follows: will beDivided into overlapping N _p The size is +.>N represents the number of pixels of an image block, m is the number of similar blocks, i=1,.. _p Is a block index. Let the reference picture block be +.>In a search window with the size of s x s, selecting m image blocks with the smallest Euclidean distance with a reference block, vectorizing each image block column and forming a similar block group->Is a column of (c).

Using the FISTA technique, problem (6) can be translated into solving a gradient problem, a rank minimization problem and two acceleration steps:

in the formulae (7) - (10), the superscript "for variables" ^k+1” and “^k "represents the k+1st iteration and the kth iteration of the algorithm.Representing intermediate variables related to gradient problem, +.>And->Representing an intermediate variable related to the acceleration step, L is +.>Is a parameter, μ=α/L. z ^k+1 ，t ^k+1 and />Is the intermediate variable of the k+1st iteration, z ^k ，t ^k and />Is the intermediate variable of the kth iteration. X is x ^k Representing the image to be reconstructed for the kth iteration, x ^k+1 Representing the image to be reconstructed for the k+1th iteration.

In (8), the similarity of the ith reference image block of the jth component image is introducedBlock group V _ji (x) Corresponding auxiliary variableCorresponding Lagrangian multiplier +.>Using the ADMM method, problem (8) can be solved by converting into the following sub-problems:

in the formula (11) { D _ji Comprises J×N _p Sub-matrix, all { D } _ji Sequentially and horizontally spliced into a matrix Is an auxiliary variable for the k+1st iteration, ">Is the k-th iteration corresponding to +.>Lagrangian multiplier of (c). In formula (12), β > 0 is a parameter. In the formula (13) { d _ji Comprises J×N _p Sub-matrix, all { d } _ji Sequentially spliced into matrix ++> Is the k+1 iteration corresponding to +.>Is a parameter that updates the lagrangian multiplier.

The low rank approximation problem (11) is a non-deterministic polynomial difficulty (nondeterministic polynomial time (NP) -hard) problem. Using the weighted kernel norms as a proxy function of rank, the problem can be rewritten as a weighted kernel norms minimization (weighted nuclear norm minimization, WNNM) problem:

wherein δ² ＝2μ/β，Representation D _ji Is defined as:

wherein σ_l (D _ji ) Representation D _ji Is the first singular value of (c). Weight w _l The calculation formula of (2) is as follows:

wherein b₀ Epsilon=10 as a parameter ^-16 For ensuring that the denominator is not 0.The estimation can be obtained by the following formula:

wherein Representation->I=1,..min (m, n) represents the singular value index.

Assume thatIs->Is a singular value decomposition (singular value decomposition, SVD), wherein Σ is a singular value diagonal matrix, andrespectively indicate->Both the 1 st and last singular values, Γ and y, are unitary matrices, with columns being left and right singular vectors of singular values, respectively. Then, the solving process of the WNNM problem (14) is as follows:

wherein soft (·) is the soft threshold operator:

soft(Σ,w)＝max(Σ-w,0) (19)

wherein w＝[w₁ ,...,w _min(m,n) ]Represents the weight, w ₁ And w is equal to _min(m,n) The 1 st and last value of w, respectively, the first value of w being w _l 。

The sub-problem with x is then solved. According to equation (12), the closed form solution for x is given by:

wherein Is V _ji The superscript "×" is the companion operator, indicating that +.>And replacing the corresponding position in the original multi-component image. />Is a diagonal matrix, each diagonal element of which is equal to the corresponding pixel point appears in all { V } _ji (x ^k ) Number of times in }.

After the reconstruction problem (6) about NLR-ESPIRiT is solved, a solution image of each group of sensitivity mapping images is obtained. To obtain a visual representation of the images, a square root of sum of squares (square root of sum of squares, SOS) method may be used on groups of images to form computationally reconstructed single coil amplitude imagesThe calculation formula is as follows:

a new ESPIRiT parallel MRI reconstruction algorithm with NLR regularization term, called the NLR-ESPIRiT algorithm, was then obtained. When the Relative Error (RE) is below the tolerance tol, the algorithm stops. X is X ^k+1 and X^k RE is defined as:

wherein X^k Representing the kth reconstructed single coil amplitude image.

The specific flow of the invention is shown in figure 1, wherein the steps are as follows:

s0: initializing, lettingz ⁰ ＝0,D ⁰ ＝0,d ⁰ ＝0,t ⁰ =0, k=0; wherein the superscript is " ⁰ "represents the initial value before iteration, x ⁰ Representing the initial value, z, of the reconstructed image x ⁰ ，t ⁰ And->Is corresponding to the intermediate variables z, t and +.>Initial value of D ⁰ Is the initial value of the auxiliary variable D, D ⁰ Is the initial value of the lagrange multiplier d.

S1, calculating z ^k+1 Formulas such as (7);

s2, initializing, j=1;

s3, initializing, i=1;

s4, for x ^k Performing similar block matching to obtain a similar block group V _ji (x ^k )；

S5, pairingSingular value decomposition to obtain ++>

S6, initializing estimationFormulas such as (17);

s7, calculating weight w= [ w ] ₁ ,...,w _min(m,n)], wherein w_l By calculation (16);

s8, calculating auxiliary variablesFormulas such as (18):

s9, judging whether N in each coil is completed _p D (D) _ji When i=n _p Step S10 is entered at the moment; otherwise, let i=i+1 and return to S4;

s10, judging whether all image components are finished, and if j=J, proceeding to step S11; otherwise, let j=j+1 and return to S3;

s11, calculating x ^k+1 Formulas such as (20);

s12, updating d ^k+1 Formulas such as (13);

s13, calculating X ^k+1 Formula (21):

s14: calculating t ^k+1 Formulas such as (9);

s15: calculation ofFormulas such as (10);

s16, calculating RE, wherein the formula is (22);

s17: judging whether an algorithm stopping standard is met, if so, entering step S18; otherwise, setting k=k+1, and returning to S1;

s18: outputting the final reconstructed single-coil amplitude image X ^k+1 。

Experimental results:

in the following experiments, NLR-ESPIRiT was compared with L1-ESPIRiT and LPJTV-ESPIRiT methods in order to verify the performance of the proposed method NLR-ESPIRiT. All algorithms are implemented with MATLAB. All experiments were performed on a server configured to: intel Core i7-5500U@2.4GHz dual Core, 12GB memory, windows 7 system (64 bit).

A simulation experiment was performed using 1 brain magnetic resonance imaging dataset with FOV smaller than the subject, named dataset1 (as shown in fig. 2). To generate the test dataset, the fully sampled dataset was manually undersampled using a two-dimensional poisson disk undersampling mask with a 5-fold Acceleration Factor (AF) and a 24 x 24 center fully sampled self-correcting region to obtain undersampled k-space data, fig. 3 illustrates the undersampling mask used.

Parameters of the compared algorithm are adjusted for SNR optima. In the following experiments, the signal-to-noise ratio (Signal Noise Ratio, SNR) was used to quantitatively evaluate the quality of the reconstructed image. SNR is calculated in a region of interest (region of interest, ROI) of an image, the higher the value of SNR, the better the image reconstruction quality. For reference image X and reconstructed imageSNR (dB) is defined as:

wherein Var represents the reference image X and the reconstructed imageMean square error between, MSE is the variance of reference image X.

To compare the reconstruction performance of each algorithm, dataset1 was undersampled using an undersampled mask of 5-fold acceleration factor, and the reconstruction was performed using all the compared algorithms. Figures 4-6 show reconstructed images using the NLR-ESPIRiT algorithm, the L1-ESPIRiT algorithm, and the LPJTV-ESPIRiT algorithm, respectively. As can be seen from the reconstructed images in fig. 4-6, the reconstructed image using the L1-ESPIRiT algorithm has significant blurring artifacts, the method using the LPJTV-ESPIRiT can further remove the artifacts, and the reconstructed image using the NLR-ESPIRiT has minimal reconstruction artifacts, which are clearer and closer to the original image.

To more intuitively compare the reconstruction effects of several algorithms, the experiments also showed corresponding reconstructed error images (white point indicates error, the more white points, the greater the error). Figures 7-9 show reconstructed error images of the NLR-ESPIRiT algorithm, the L1-ESPIRiT algorithm, and the LPJTV-ESPIRiT algorithm, respectively. As apparent from the error image, the reconstructed image of the L1-ESPIRiT has larger error, and then the reconstructed image of the L1-ESPIRiT has smaller reconstruction error and better reconstruction quality.

In addition, experiments also analyzed the values of SNR of reconstructed images of several methods. The SNR of the reconstructed image of the L1-ESPIRiT method (as in fig. 5) is 17.74dB; the SNR of the LPJTV-ESPIRiT method (e.g., fig. 6) is 19.37dB; the SNR of the reconstructed image (as in fig. 4) of the NLR-ESPIRiT method is 20.27dB. It can be seen that the NLR-ESPIRIT method provided by the invention can obviously improve the SNR of the reconstructed image.

In summary, the selected test set is used to experiment different algorithms, and the NLR-ESPIRiT method based on non-local low-rank regularization is obviously superior to the L1-ESPIRiT and LPJTV-ESPIRiT algorithms in SNR evaluation index and vision.

While the present invention has been described in detail with reference to the drawings, the present invention is not limited to the above embodiments, and various changes can be made without departing from the spirit of the present invention within the knowledge of those skilled in the art.

Claims (1)

1. A feature vector-based self-consistency and non-local low-rank parallel MRI reconstruction method, the method comprising the steps of:

s0: initializing, lettingz ⁰ ＝0,D ⁰ ＝0,d ⁰ ＝0,t ⁰ ＝0，k＝0；

Wherein the superscript "0" represents the initial value before iteration, and the superscript "T" represents the conjugate transpose operator;representing a column-vectorized multicomponent image to be reconstructed, n=n _h ×N _v Pixel count, N, representing a single coil amplitude image _v and N_h The number of lines and columns of the single coil amplitude image respectively, J represents the number of sensitivity groups, the number of image components of x is equal to the number of sensitivity groups, and the J-th component image of x is x _j The expression j=1,..j, J, the index of the sensitivity map group, x ₁ and x_J Respectively representing the 1 st and the J-th component images of a column-vectorized multicomponent image x to be reconstructed, x ⁰ The initial value of x is indicated,representing column-vectorized undersampled multi-coil k-space data, M representing the number of single-coil k-space undersampled data points, C representing the number of receive coils used for parallel imaging, y for the C-th coil data of y _c And c=1,.. representing the coil index, y ₁ and y_C 1 st coil data and C-th coil data representing column-vectorized undersampled multi-coil k-space data y, respectively, +.> and />Respectively representing a Fourier operator and an undersampled operator, I _C Is a matrix of units of CxC, < >> and />Respectively represent N _h and N_v Point Fourier transform matrix, < >>Matrix representing the selection of sampling point positions from a k-space grid,/for>Representation of CroneckerAccumulation of pathogenic qi>For the sensitivity map group matrix, the j-th sensitivity map of the c-th coil of S is obtained by using the SPIRIT model based on the feature vector _cj Representation, S ₁₁ Group 1 sensitivity maps representing the 1 st coil, S _CJ Group J sensitivity map representing coil C,>and->Representing intermediate variables, z ⁰ ，t ⁰ And->Is their corresponding initial value, defines the block matching operator V _ji :/>The operator performs the block matching on x as follows: will->Divided into overlapping N _p The size is +.>N represents the number of pixels of the image block, m is the number of similar blocks, i=1,.. _p For block index, reference picture block is +.>In a search window with the size of s x s, selecting m image blocks with the smallest Euclidean distance with a reference block, vectorizing each image block column and forming a similar block group->Column of->Similar block group V representing the ith reference image block of the jth component image _ji (x) Corresponding auxiliary variables, all { D } _ji Sequentially and horizontally spliced into a matrix +.>D ⁰ Is the initial value of D, +.>Is corresponding to D _ji Lagrangian multipliers of (2), all { d }, will be _ji Sequentially spliced into a matrixd ⁰ Is the initial value of d;

s1, calculating the intermediate variable z of the (k+1) th iteration ^k+1 The calculation formula is as follows:

wherein ,intermediate variable representing the kth iteration +.>L represents Lipschitz constant;

s2, initializing, j=1;

s3, initializing, i=1;

s4, reconstructing an image x of the kth iteration ^k Performing similar block matching to obtain a similar block group V _ji (x ^k )；

S5, pairingSingular value decomposition to obtain ++>

Wherein, sigma is a singular value diagonal matrix,both Γ and gamma are unitary matrices, the columns of which are left and right singular vectors of singular values, respectively,/->Is expressed as the first singular value of (2) and />Respectively indicate->I=1,..min (m, n), which is a singular value index, +.>Representing the corresponding auxiliary variable in the kth iterationLagrangian multipliers of (2);

s6, initializing estimation by the following formula

S7, calculating weight w= [ w ] ₁ ,...,w _min(m,n) ]，w ₁ And w is equal to _min(m,n) The 1 st and last value of w, respectively, the first value of w being w _l The calculation formula is as follows:

wherein b₀ Epsilon=10 as a parameter ^-16 For ensuring that the denominator is not 0;

s8, calculating auxiliary variables of the (k+1) th iterationThe calculation formula is as follows:

wherein: soft (Σ, w) =max (Σ -w, 0) is the soft threshold operator;

s9, judging whether N in each coil is completed _p Personal (S)When i=n _p Step S10 is entered at the moment; otherwise, let i=i+1 and return to S4;

s10, judging whether all image components are finished, and if j=J, proceeding to step S11; otherwise, let j=j+1 and return to S3;

s11, calculating an image x to be reconstructed of the (k+1) th iteration ^k+1 The calculation formula is as follows:

wherein the superscript "×" is the accompanying operator, the superscript "-1" represents the inverting operator, is V _ji The companion matrix representing ∈>Put back the corresponding position in the original multicomponent image, < > in->Is a diagonal matrix, each diagonal element of which is equal to the corresponding pixel point appearing in all the groups of similar blocks V _ji (x ^k ) The number of times in (a);

s12, updating d of the (k+1) th iteration ^k+1 ，d ^k+1 Is a sub-matrix of (2)The calculation can be performed by the following formula:

wherein η > 0 represents a parameter for updating the Lagrangian multiplier;

s13, reconstructed multicomponent image x ^k+1 The reconstructed single coil amplitude image is calculated using the square root of the sum of squares (square root of sum of squares, SOS) as follows:

wherein A single coil amplitude reconstruction image representing the k+1th iteration;

s14: updating intermediate variables for the k+1st iteration

S15: calculating intermediate variables for the k+1th iteration

S16, calculating X ^k+1 and X^k Relative Error (RE) between (a) and (b) is calculated as:

wherein Reconstructing an image of single coil amplitudes representing the kth iteration;

s17: judging whether an algorithm stopping standard is met, if so, entering step S18; otherwise, setting k=k+1, and returning to S1;

s18: outputting the final reconstructed single-coil amplitude image X ^k+1 。