CN113892938A - Improved sensitivity encoding reconstruction method based on non-local low-rank constraint - Google Patents

Abstract

The invention relates to an improved sensitivity encoding reconstruction method based on non-local low-rank constraint, and belongs to the technical field of magnetic resonance imaging. Parallel Magnetic Resonance Imaging (MRI) reconstruction has been a research focus in recent years, and the sensitivity encoding (SENSE) model is a widely used parallel MRI reconstruction model that reconstructs images from multi-coil undersampled k-space data by directly using coil sensitivity information. The invention provides a parallel magnetic resonance imaging reconstruction algorithm based on image non-local low-rank (NLR) constraint based on a SENSE reconstruction model, which is named as NLR-SENSE, uses a weighted nuclear norm as a rank proxy function, and uses a multiplier Alternating Direction (ADMM) method to effectively solve an objective optimization function. Compared with an LpTV-SENSE algorithm combined with Lp pseudo-norm Total Variation (TV) constraint, simulation experiment results show that the NLR-SENSE algorithm can obtain a clearer reconstructed image and is an effective SENSE improved algorithm.

Description

Improved sensitivity encoding reconstruction method based on non-local low-rank constraint

Technical Field

The invention relates to an improved sensitivity encoding reconstruction method based on non-local low-rank constraint, and belongs to the technical field of magnetic resonance imaging.

Background

Magnetic Resonance Imaging (MRI) is a few radiation-free medical diagnostic imaging methods at present, and has wide application in clinical medicine. However, magnetic resonance imaging has the disadvantages of long time and more imaging artifacts, and cannot meet the requirements of efficient and accurate medical diagnosis. Therefore, the research of the fast magnetic resonance imaging technology and the high-quality magnetic resonance imaging reconstruction technology has high practical significance.

Parallel imaging is a common method for accelerating MRI, which uses a multi-channel phased array to simultaneously acquire magnetic resonance signals, and uses the spatial sensitivity difference of each coil to perform sensitivity encoding of spatial information. Sensitivity Encoding (SENSE) is a widely used parallel imaging technique that recovers image information by directly using information based on coil sensitivity.

The process of reconstructing an image from undersampled k-space data is often expressed as a goal optimization problem. But data reconstruction is usually a morbid problem, and direct solution may have more aliasing artifacts. In order to obtain a better reconstruction result, the past research proposes that a regularization term with data prior information is added into a target optimization problem, so that the performance of a reconstruction algorithm can be effectively improved. At present, the parallel magnetic resonance imaging reconstruction algorithm based on the SENSE model has the following steps: the L1-SENSE algorithm combined with the L1 norm regularization term uses sparsity of the image in projection subspace (e.g., wavelet domain) as prior information; a TV-SENSE algorithm combined with a Total Variation (TV) regularization term, using local similarity characteristics of images as prior information; and an Lp-SENSE algorithm incorporating an Lp pseudo-norm total variation (LpTV) regularization term, the use of an Lp pseudo-norm can further improve algorithm performance. However, images reconstructed by using the algorithms still have more artifacts, and the reconstruction performance of the algorithms still needs to be improved.

Disclosure of Invention

The invention aims to overcome the defects of the prior art and provide an improved sensitivity encoding reconstruction method based on non-local low-rank constraint, which can further improve the quality of a parallel magnetic resonance imaging reconstruction image.

The purpose of the invention is realized by the following technical scheme: an improved sensitivity encoding reconstruction method based on non-local low-rank constraint. It is characterized by comprising the following steps:

s0: initialization, order

Z⁰＝0，z⁰＝0，D⁰＝0，d⁰＝0，k＝0；

Where the superscript "0" represents the initial value before iteration and the superscript "H" represents the conjugate transpose operator.

Representing an image to be reconstructed with column vectorization, N ═ N_h×N_vNumber of pixel points, N, representing a single-coil amplitude image_hAnd N_vNumber of rows and columns, X, respectively, of the single-coil amplitude image⁰Is the initial value of X.

The method comprises the steps of representing column vectorized undersampled multi-coil k-space data, wherein a superscript T represents a transposition operator, M represents the number of undersampled data points of a single coil k-space, and C represents the number of receiving coils used for parallel imaging. For c-th coil data of Y

Denotes, C1, C denotes a coil index, Y denotes₁And Y_CThe 1 st coil data and the C th coil data of the column-vectorized under-sampled multi-coil k-space data Y are respectively represented.

The representation of the fourier operator is shown as,

representing an undersampling operator. I is_CIs a C x C identity matrix and,

and

are each N_hAnd N_vA two-dimensional discrete fourier transform matrix of points,

a matrix representing the selection of sample point locations from a k-space grid,

representing the kronecker product. Wherein

The sensitivity operator is represented. For sensitivity information of the c-th coil

Denotes S_cIs a diagonal matrix. S₁Indicating sensitivity information of the 1 st coil, S_CRepresenting the C-th coil sensitivity information.

Denotes an auxiliary variable, Z⁰Is the initial value of Z.

Representing the corresponding Lagrangian multiplier, Z⁰Is the initial value of z.

Defining Block matching operator V_i:

The operator performs block matching on X as follows: will be provided with

Divided into overlapping N_pEach size is

Reference image block of, N_pFor the number of reference image blocks, N represents the number of pixels of the image block, m is the number of similar blocks, and i is 1_pIs a block index. Assume that the reference image block is

In a neighborhood search window of size sxs (e.g. 40 × 40), m image blocks with the minimum euclidean distance to the reference block are selected, each image block column is vectorized and a similar block group V is constructed_iAnd (X) column.

Similar block group V for representing ith reference image block_i(X) corresponding auxiliary variables, all { D }_iHorizontally splicing into a matrix in sequence

D⁰Is the initial value of D.

Is corresponding to D_iLagrange multiplier of (d) will all of d_iHorizontally splicing into a matrix in sequence

d⁰Is the initial value of d.

S1: initializing, i is 1;

s2: judging whether mod (k, T) is 0 or not, and if so, carrying out reconstruction on the k iteration X^kSimilar block matching is carried out to obtain similar block groups V_i(X^k) Otherwise go to S3;

where mod (·) represents the remainder calculation.

S3: to pair

Performing singular value decomposition to obtain

Wherein, Σ is a singular value diagonal matrix, and the jth diagonal element thereof is the jth singular value:

j＝1,...,min(m, n) is singular value index, and all singular values are combined into a row vector

And

respectively represent

Y is a unitary matrix whose columns are left singular vectors of singular values, Γ is a unitary matrix whose columns are right singular vectors of singular values,

indicating correspondence to auxiliary variables in the k-th iteration

Lagrange multiplier.

S4: calculating the weight w ═ w₁,...,w_min(m,n)]Wherein w is₁And w_min(m,n)Representing the first and last values of w, respectively. w is calculated as:

wherein ε is 10^-16To ensure that the denominator is not 0.

S5: computing an auxiliary variable D for the (k + 1) th iteration^k+1The calculation formula is as follows:

D^k+1＝Υdiag(soft(σ,τw))Γ^H；

wherein: soft (σ, τ w) max (σ - τ w,0) is a soft threshold operator, τ μ/γ₁Mu > 0 and gamma₁> 0 is a parameter.

S6: judging whether N is completed_pAn

When i is equal to N_pThe routine proceeds to step S7; otherwise, returning to S2 after i is equal to i + 1;

s7: calculating an auxiliary variable Z for the (k + 1) th iteration^k+1The calculation formula is as follows:

wherein the superscript "-1" represents the operator for the matrix inversion,

is a diagonal matrix of the grid,

represents an identity matrix, z^kIs the k-th iteration corresponding to the auxiliary variable Z^kOf lagrange multiplier, gamma₂> 0 is a parameter.

S8: calculating the image X to be reconstructed of the (k + 1) th iteration^k+1The calculation formula is as follows:

wherein

Representation operator V_iIs associated with an operator, represents

The corresponding position in the image is put back.

Is a diagonal matrix whose each diagonal element is equal to the corresponding pixel point appearing in all similar groups of blocks { V }_i(X^k) The number of times in (c).

Is a diagonal matrix.

S9: updating Lagrange multiplier z for the (k + 1) th iteration^k+1Calculated by the following formula:

z^k+1＝z^k+η₂(SX^k+1-Z^k+1)；

wherein eta is₂The parameter of the Lagrange multiplier z is updated when the value is more than 0;

s10: initializing, i is 1;

s11: updating Lagrange multiplier for the (k + 1) th iteration

The calculation is made by the following formula:

wherein eta is₁> 0 denotes updating the Lagrangian multiplier d_iParameters of (1) };

s12: judging whether N is completed_pAn

When i is equal to N_pThe routine proceeds to step S13; otherwise, returning to S11 after i is equal to i + 1;

s13: judging whether the maximum iteration number is reached, and if K is satisfied, entering step S14; otherwise, k is set to k +1, and the process returns to S2;

s14: outputting the final reconstructed image X^k+1。

The invention has the beneficial effects that: the SENSE model is a parallel magnetic resonance imaging reconstruction algorithm directly using coil sensitivity information, and converts the process of reconstructing an image from multi-coil undersampled k-space data into the solution of a least square optimization problem. The invention combines the regularization constraint of non-local low-rank (NLR), provides a high-quality SENSE improved algorithm which is called as NLR-SENSE algorithm, and can further improve the performance of the parallel magnetic resonance imaging reconstruction algorithm. For the objective optimization problem of the proposed algorithm, the present invention uses a multiplier alternating direction method (ADMM) and solves it using the weighted kernel norm as the proxy function of the rank. Simulation experiments show that compared with the LpTV-SENSE algorithm, the algorithm provided by the invention has higher parallel magnetic resonance image reconstruction performance.

Drawings

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

FIG. 2 is an image of a brain (i.e., dataset 1) obtained using a full sampling of eight channel receive coils;

FIG. 3 is a two-dimensional Poisson disk undersampling mask diagram with a 5 times acceleration factor for a 24 × 24 center fully sampled self-correcting region;

FIGS. 4 and 5 show images reconstructed from undersampled dataset 1 data from a two-dimensional Poisson disk undersampled mask with a 5 times acceleration factor for a 24 × 24 center fully sampled self-correcting region using the LpTV-SENSE algorithm and the NLR-SENSE algorithm, respectively;

fig. 6 and 7 show the error images corresponding to fig. 4 and 5, respectively, i.e. the error images resulting from the LpTV-SENSE algorithm and the NLR-SENSE reconstruction.

Detailed Description

The technical scheme of the invention is further described in detail in the following with reference to the attached drawings.

Example 1: the invention provides an improved sensitivity encoding reconstruction method based on non-local low-rank constraint based on a SENSE framework, which comprises the following steps:

1) SENSE framework:

suppose that

Representing an image to be reconstructed with column vectorization, N ═ N_h×N_vNumber of pixel points, N, representing a single-coil amplitude image_hAnd N_vThe number of rows and columns, respectively, of the single-coil amplitude image.

Under-sampled multi-coil k-space data representing column vectorization, the superscript "T" representing the transpose operator, and M representing single-coil k-space undersamplingThe number of sample data points, C, represents the number of receive coils used for parallel imaging. For c-th coil data of Y

Denotes, C1, C denotes a coil index, e.g., Y₁And Y_CThe 1 st coil data and the C th coil data of the column-vectorized under-sampled multi-coil k-space data Y are respectively represented.

And

are each N_hAnd N_vA two-dimensional discrete fourier transform matrix of points,

a matrix representing the selection of sample point locations from a k-space grid.

And

the fourier operator and the undersampling operator are indicated separately.

Denotes the kronecker product, I_CIs a C × C identity matrix.

In the SENSE reconstruction model, the sensitivity operator is

Indicating sensitivity information of the c-th coil, S_cIs a diagonal matrix. S₁Indicating sensitivity information of the 1 st coil, S_CRepresenting the C-th coil sensitivity information. The relationship between the image X to be reconstructed and the k-space undersampled image Y can then be found as:

the solution of the SENSE model is often written as the following least squares problem:

2) and (3) algorithm derivation:

although the SENSE algorithm combined with the regularization methods of L1, TV, and LpTV can improve the reconstruction quality to some extent, there are still many artifacts in the image reconstructed using these algorithms, and there is still a large room for improvement in the image reconstruction performance of the SENSE technique. In order to further improve the reconstruction quality, the invention combines non-local low rank regularization (NLR) with a SENSE model to provide an NLR-SENSE model, and obtains the following optimization problem:

wherein

In the case of a data fidelity item,

is a regularization term, rank (V)_i(X)) represents V_iRank of (X), μ > 0, represents a regularization parameter.

Is a block matching operator, the step of block matching X is as follows: dividing X into overlapping N_pEach size is

Reference image block of, N_pRepresenting the number of reference image blocks, N representing the number of pixel points of the image blocks, i ═ 1_pIs a block index. Assume that the reference image block is

In a neighborhood search window (for example, setting the size of a local window to be 40 × 40), m image similar blocks with the minimum Euclidean distance from a reference block are found, m is the number of the similar blocks, each image block column is vectorized and forms a similar block group V_iAnd (X) column.

Introducing auxiliary variables

And

and corresponding lagrange multiplier

And

using the ADMM technique, problem (3) can be transformed into the following sub-problem to be solved with alternate updates:

z^k+1＝z^k+η₂(SX^k+1-Z^k+1) (7)

in the formulas (4) to (8), the superscripts "k + 1" and "k" of the variables denote the k +1 th of the algorithmThe sub-iteration and the k-th iteration. In formula (4) { D_iContains N_pSub-matrix, all of { D_iHorizontally splicing into a matrix in sequence

Is an auxiliary variable for the (k + 1) th iteration,

is in the kth iteration corresponds to

Of lagrange multiplier, gamma₁> 0 is a parameter. In formula (5), Z^k+1Is an auxiliary variable for the (k + 1) th iteration, z^kIs the k-th iteration corresponding to Z^kOf lagrange multiplier, gamma₂> 0 is a parameter. In formula (6), X^k+1Is the image to be reconstructed for the (k + 1) th iteration. In formula (7), z^k+1Is the k +1 th iteration corresponding to Z^k+1Of lagrange multiplier, η₁> 0 is a parameter that updates the lagrange multiplier z. In formula (8) { d_iContains N_pSub-matrix, all of d_iSplicing into a matrix in sequence

Is in the k +1 th iteration corresponds to

Of lagrange multiplier, η₂> 0 is the updated Lagrangian multiplier d_iThe parameters of (c).

1) Solving for problems with { D_iThe low rank approximate minimization problem of }: the low rank approximation problem (4) is a nondeterministic polynomial difficulty (NP) -hard) problem. Using the weighted kernel norm as a proxy function for the rank, the problem can be rewritten as:

wherein τ ═ μ γ₁，||.||_w,*Represents a weighted nuclear norm, defined as:

wherein sigma_jIs the j-th singular value of the matrix, w_jThe representation corresponds to σ_jThe weight of (c).

Suppose that

Is that

Singular Value Decomposition (SVD) of (1), wherein

Is a singular value diagonal matrix, and the superscript "H" represents the conjugate transpose operator. Forming a row vector from all singular values

Is expressed as a j-th singular value of

j 1.. min (m, n) denotes a singular value index,

and

respectively represent

The 1 st and last singular values of y and Γ are unitary matrices, the columns of which are the left and right singular vectors of the singular values, respectively. Then, a solution to the weighted nuclear norm minimization problem (9) is obtained as:

D^k+1＝Υsoft(Σ,τw)Γ^H (11)

where soft (Σ, τ w) is max (Σ -diag (τ w),0) denotes a soft threshold operator, and w is [ w ═ w [ [ w ]₁,w₂,...,w_min(m,n)]To represent

Weight corresponding to the singular value of, w₁And w_min(m,n)Respectively the 1 st and last value of w. w is calculated as:

wherein ε is 10^-16To ensure that the denominator is not 0.

2) Solving the sub-problem with respect to Z: according to equation (5), the solution to Z is:

where the superscript "-1" denotes the operator that inverts the matrix,

is a diagonal matrix of the grid,

representing an identity matrix.

3) Solving the sub-problem with respect to X. According to equation (6), the solution for X is given by:

wherein

Is V_iThe suffix ". sup." is the companion operator, indicating that

The corresponding position in the image is put back.

Is a diagonal matrix whose each diagonal element is equal to the corresponding pixel point occurring at all { V }_i(x^k) The number of times in (c).

Is a diagonal matrix.

By solving (4) - (8) circularly and alternately, the final reconstructed image can be obtained, and a new SENSE parallel MRI reconstruction algorithm with NLR regularization, namely an NLR-SENSE algorithm, is obtained. Because the block matching algorithm has higher computational complexity, in the NLR-SENSE algorithm, the invention proposes to update the block matching operator every T times, that is, block matching is performed when mod (k, T) is judged to be 0, where mod (·) represents remainder calculation.

The specific flow is shown in fig. 1, wherein the steps are as follows:

s0: initialization, order

Z⁰＝0，z⁰＝0,D⁰＝0,d⁰＝0，k＝0；

Where the superscript "0" indicates the initial value before the iteration. X⁰Representing an initial value, Z, of an image X to be reconstructed⁰And D⁰Initial values of auxiliary variables Z and D, respectively, Z⁰And d⁰Are initial values of Z and d, respectively, i.e. corresponding to Z, respectively⁰And D⁰Lagrange multiplier.

S1: initializing, i is 1;

s2: judging whether mod (k, T) is 0 or not, and if so, carrying out reconstruction on the k iteration X^kSimilar block matching is carried out to obtain similar block groups V_i(X^k) Otherwise go to S3;

s3: to pair

Performing singular value decomposition to obtain

S4: calculating a weight w by the equation (12);

s5: computing an auxiliary variable D for the (k + 1) th iteration^k+1Calculating by the formula (11);

s6: judging whether N is completed_pAn

When i is equal to N_pThe routine proceeds to step S7; otherwise, returning to S2 after i is equal to i + 1;

s7: calculating an auxiliary variable Z for the (k + 1) th iteration^k+1Calculating by the formula (13);

s8: calculating the image X to be reconstructed of the (k + 1) th iteration^k+1Calculated by equation (14);

s9: updating Lagrange multiplier z for the (k + 1) th iteration^k+1Calculating by formula (7);

s10: initializing, i is 1;

s11: updating Lagrange multiplier for the (k + 1) th iteration

Calculating by formula (8);

s12: judging whether N is completed_pAn

When i is equal to N_pThe routine proceeds to step S13; otherwise, returning to S11 after i is equal to i + 1;

s13: judging whether the maximum iteration number is reached, and if K is satisfied, entering step S14; otherwise, k is set to k +1, and the process returns to S2;

s14: outputting the final reconstructed image X^k+1。

The experimental results are as follows:

the NLR-SENSE method was compared with the LpTV-SENSE method in experiments to verify the reconstruction performance of the algorithm. All algorithms were implemented with MATLAB. All experiments were performed on a laptop configured to: intel core i7-5500U @2.4GHz dual core, 12GB memory and Windows 7 operating system (64 bits).

The reconstruction of the SENSE model needs to rely on coil sensitivity encoding information, but accurate sensitivity information is not readily available in practical applications. The coil sensitivity is estimated by using a sensitivity estimation method provided by an ESPIRiT model, the method depends on an auto-calibration signal (ACS) and estimates the sensitivity by using a characteristic vector decomposition method, and the method is a more accurate sensitivity estimation method. It should be noted that, any sensitivity estimation method may be used in the algorithm provided by the present invention to obtain the sensitivity information.

A simulation experiment was performed using 1 brain magnetic resonance imaging dataset acquired with 8 channels, named dataset 1 (as shown in figure 2). To generate the test dataset, the fully sampled dataset was artificially undersampled using a two-dimensional poisson disk undersampled mask with a 5 x acceleration factor of a 24 x 24 center fully sampled self-correcting region to obtain undersampled k-space data, fig. 3 illustrates the undersampled mask used.

In the following experiment, the quality of the reconstructed image was evaluated using a Signal-to-Noise Ratio (SNR) as a numerical index. A higher value of SNR indicates a better image reconstruction quality. In the experiment, the SNR index is calculated only in the region of interest (ROI) of the image. The parameters of all comparison algorithms are adjusted to SNR optimum. SNR (dB) is defined as:

where Var denotes the reference image X and the reconstructed image

Mean square error between, MSE is the variance of the reference image X.

And acquiring undersampled data by using an undersampled mask of 5 times of acceleration factor for the fully-sampled data dataset 1 acquired by using multiple coils, and reconstructing the undersampled data by using all comparative algorithms. In order to compare the reconstruction effects of different algorithms, experiments show the reconstructed images of different algorithms and error graphs between the reconstructed images and an original image (white dots indicate errors, and the more the white dots, the larger the errors). Specifically, fig. 4 and 5 show images reconstructed using the LpTV-SENSE algorithm and the NLR-SENSE algorithm, respectively, and fig. 6 and 7 show reconstructed error images corresponding to fig. 4 and 5. From the reconstructed image and the error image, it can be seen that the image reconstructed by the LpTV-SENSE method has more artifacts and more errors, and the reconstructed image by the NLR-SENSE has less reconstruction artifacts and reconstruction errors, so that the reconstructed image is clearer and is closer to the original image.

In addition, the experiment also analyzed the values of SNR of the reconstructed images of the two methods. The SNR of the LpTV-SENSE method (as in FIG. 4) was 21.27 dB; the SNR of the reconstructed image (as in fig. 5) of the NLR-SENSE method is 22.71 dB. It can be seen that the NLR-SENSE method proposed by the present invention has a high SNR, which is consistent with the visual image and the error image of the reconstructed result.

In conclusion, different algorithms are used for carrying out experiments on the selected test set, and the NLR-SENSE method based on non-local low-rank regularization is obviously superior to the LpTV-SENSE algorithm in SNR evaluation index and visual comparison.

While the present invention has been described in detail with reference to the embodiments shown in the drawings, the present invention is not limited to the embodiments, and various changes can be made without departing from the spirit and scope of the present invention.

Claims (1)

1. A method for improving sensitivity coding reconstruction based on non-local low rank constraint is characterized by comprising the following steps:

s0: initialization, order

Z⁰＝0，z⁰＝0，D⁰＝0，d⁰＝0，k＝0；

Wherein the superscript "0" represents the initial value before iteration, the superscript "H" represents the conjugate transpose operator,

representing an image to be reconstructed with column vectorization, N ═ N_h×N_vNumber of pixel points, N, representing a single-coil amplitude image_hAnd N_vNumber of rows and columns, X, respectively, of the single-coil amplitude image⁰Is the initial value of X and is,

the method comprises the steps of representing column vectorization undersampled multi-coil k-space data, using superscript T to represent a transposition operator, using M to represent the number of undersampled data points of a single-coil k-space, using C to represent the number of receiving coils used for parallel imaging, and using the data of the C-th coil of Y

Denotes, C1, C denotes a coil index, Y denotes₁And Y_CThe 1 st coil data and the C th coil data respectively representing column-vectorized under-sampled multi-coil k-space data Y,

the representation of the fourier operator is shown as,

representing undersampling operators, I_CIs a C x C identity matrix and,

and

are each N_hAnd N_vPoint two-dimensional discrete fourier transformThe matrix is a matrix of a plurality of matrices,

a matrix representing the selection of sample point locations from a k-space grid,

represents the kronecker product, wherein

Indicating sensitivity operator, for sensitivity information of the c-th coil

Denotes S_cIs a diagonal matrix, S₁Indicating sensitivity information of the 1 st coil, S_CThe C-th coil sensitivity information is represented,

denotes an auxiliary variable, Z⁰Is the initial value of Z and is,

representing the corresponding Lagrangian multiplier, Z⁰Is the initial value of z;

defining Block matching operator V_i:

The operator performs block matching on X as follows: will be provided with

Divided into overlapping N_pEach size is

Reference image block of, N_pFor the number of reference image blocks, N represents the number of pixels of the image block, m is the number of similar blocks, and i is 1_pFor the block index, a reference image block is set to

In a neighborhood search window with the size of s multiplied by s, m image blocks with the minimum Euclidean distance from a reference block are selected, each image block column is vectorized, and a similar block group V is formed_iThe column of (X) is a column,

similar block group V for representing ith reference image block_i(X) corresponding auxiliary variables, all { D }_iHorizontally splicing into a matrix in sequence

D⁰Is the initial value of D and is,

is corresponding to D_iLagrange multiplier of (d) will all of d_iHorizontally splicing into a matrix in sequence

d⁰Is the initial value of d and is,

s1: initializing, i is 1;

s2: judging whether mod (k, T) is 0 or not, and if so, carrying out reconstruction on the k iteration X^kSimilar block matching is carried out to obtain similar block groups V_i(X^k) Otherwise go to S3;

where mod (-) represents the remainder calculation;

s3: to pair

Performing singular value decomposition to obtain

Wherein, Σ is a singular value diagonal matrix,the jth diagonal element is the jth singular value:

for singular value indexing, all singular values are combined into a row vector

And

respectively represent

Y is a unitary matrix whose columns are left singular vectors of singular values, Γ is a unitary matrix whose columns are right singular vectors of singular values,

indicating correspondence to auxiliary variables in the k-th iteration

Lagrange multiplier of (a);

s4: calculating the weight w ═ w₁,...,w_min(m,n)]Wherein w is₁And w_min(m,n)The first and last values of w are represented, respectively, and w is calculated as:

wherein ε is 10^-16For ensuring that the denominator is not 0;

s5: computing an auxiliary variable D for the (k + 1) th iteration^k+1The calculation formula is as follows:

D^k+1＝Υdiag(soft(σ,τw))Γ^H；

wherein: soft (σ, τ w) max (σ - τ w,0) is a soft threshold operator, τ μ/γ₁Mu > 0 and gamma₁More than 0 is a parameter;

s6: judging whether N is completed_pAn

When i is equal to N_pThe routine proceeds to step S7; otherwise, returning to S2 after i is equal to i + 1;

s7: calculating an auxiliary variable Z for the (k + 1) th iteration^k+1The calculation formula is as follows:

wherein the superscript "-1" represents the operator for the matrix inversion,

is a diagonal matrix of the grid,

represents an identity matrix, z^kIs the k-th iteration corresponding to the auxiliary variable Z^kOf lagrange multiplier, gamma₂More than 0 is a parameter;

s8: calculating the image X to be reconstructed of the (k + 1) th iteration^k+1The calculation formula is as follows:

wherein V_i ^*:

Representation operator V_iIs associated with an operator, represents

The corresponding position in the image is put back,

is a diagonal matrix whose each diagonal element is equal to the corresponding pixel point appearing in all similar groups of blocks { V }_i(X^k) The number of times in (c) },

is a diagonal matrix of the grid,

s9: updating Lagrange multiplier z for the (k + 1) th iteration^k+1Calculated by the following formula:

z^k+1＝z^k+η₂(SX^k+1-Z^k+1)；

wherein eta is₂The parameter of the Lagrange multiplier z is updated when the value is more than 0;

s10: initializing, i is 1;

s11: updating Lagrange multiplier for the (k + 1) th iteration

The calculation is made by the following formula:

wherein eta is₁> 0 denotes updating the Lagrangian multiplier d_iParameters of (1) };

s12: judging whether N is completed_pAn

When i is equal to N_pThe routine proceeds to step S13; otherwise, returning to S11 after i is equal to i + 1;

s13: judging whether the maximum iteration number is reached, and if K is satisfied, entering step S14; otherwise, k is set to k +1, and the process returns to S2;

s14: outputting the final reconstructed image X^k+1。

Images