CN109345473B - Image processing method based on self-adaptive fast iterative shrinkage threshold algorithm - Google Patents

Abstract

The invention discloses an image processing method based on a self-adaptive fast iterative shrinkage threshold algorithm, which comprises eight steps and is used for reconstructing an original image. Compared with the existing three algorithms, the image acquired by the method can show more local details and clearer outlines; fewer error points are generated, providing more accurate reconstruction; the convergence rate is high, and the iteration efficiency is higher than that of other methods; the image reconstruction of different parts has stability.

Description

Image processing method based on self-adaptive fast iterative shrinkage threshold algorithm

Technical Field

The invention belongs to the field of image processing, and particularly relates to an image processing method based on a self-adaptive fast iterative shrinkage threshold algorithm.

Background

Magnetic Resonance (MR) imaging is a safe, fast and accurate image acquisition technique. It has the advantages of multiple directions, parameters and modes, and is harmless to human body. It can display anatomical and functional information of human tissue. MR imaging has a wide range of applications. However, MR imaging has long scanning time, slow scanning speed, and may not provide dynamic real-time images and navigation due to image blurring caused by organ motion. Therefore, the disadvantages of MR imaging limit the spread of functional imaging and cause additional pain to the user.

Therefore, the skilled person in the field makes intensive studies and proposes the theory of compressed sensing: at very low nyquist sampling rates, discrete signals are obtained by random sampling. In some known transform domains, the original signal is reconstructed by a non-linear reconstruction algorithm based on the sparsity of the signal.

At present, it is very important to construct a stable and efficient reconstruction algorithm based on compressed sensing. The reconstruction algorithm mainly comprises a greedy algorithm and a convex relaxation algorithm. For low-dimensional small-scale signals, the greedy tracking algorithm is fast and has good quality, such as matching tracking, orthogonal matching tracking and regularized orthogonal matching tracking. However, it is difficult to satisfy the requirement of reconstruction accuracy for such a large-scale signal with high dimension. The convex relaxation algorithm takes less time in the reconstruction. The classical convex function optimization algorithm mainly comprises conjugate gradient, bre gman iteration and iteration reweighted least squares.

In order to improve the speed of magnetic resonance image reconstruction, an iterative shrinkage threshold algorithm and a series of improved algorithms are provided. These methods directly solve the L1 minimization problem. Daubechies et al propose an iterative algorithm (ISTA) that applies a Landweber iteration of a threshold (or non-linear contraction) in each iteration step to solve the linear inverse problem. BeckA et al propose a novel iterative shrinkage threshold algorithm (FISTA) that preserves the computational simplicity of ISTA and improves the global convergence speed. Xiaoobo Qu et al sparsely represent curves and edges using contourlet transforms to solve the L1 norm optimization problem for CS-MRI. Bayram I et al studied a subband-adaptive version of the popular iterative contraction/thresholding algorithm, which employs different update steps and thresholds for each subband. They also present an algorithm to select the appropriate update steps and thresholds for the linearity and time invariance of the distortion operator. Wang proposes a novel sparse dubbing index wavelet transform (EWT) that provides sparser coefficients than the conventional wavelet transform. Elahi et al. An improved iterative algorithm (GTIA) for CS-MRI image reconstruction based on P-threshold is introduced, using which sparsity in the image is promoted, which is a key factor for CS-based image reconstruction.

However, the above algorithm is based on the fact that the shrinking factor has been reduced by a fixed step size in the iteration, which is not applicable to the whole iterative process. In fact, we want to converge quickly in earlier iterations and reduce the contraction speed to preserve reconstruction accuracy in later iterations.

Therefore, the present application proposes an adaptive fast iterative shrink threshold algorithm (SAFISTA) to adaptively adjust the iterative shrink operator. It uses feedback to dynamically adjust the iteration step size to increase the convergence speed.

The invention content is as follows:

aiming at the defects of the prior art, the invention provides an image processing method based on an adaptive fast iterative shrinkage threshold algorithm.

The technical scheme of the invention is as follows: an image processing method based on an adaptive fast iterative shrinkage threshold algorithm comprises the following steps.

The method comprises the following steps: acquiring K-space undersampled MR image data, and reconstructing a model as follows;

F_u＝kF (2)

wherein m is C^NIs a reconstructed image, y ∈ C^MIs the K-space undersampled signal of MR, F_uIs the observation matrix, k i is the sampling pattern, and F represents the two-dimensional fourier transform.

Step two: introducing a projection coefficient theta of an original image in a sparse domain, and sparsely representing m as phi theta, phi as phi [ psi ]₁,ψ₂,...,ψ_N]∈R^N×NConverting the L0 norm minimization problem of the model into a norm L1 minimization problem;

m＝Ψθ (4)

step three: combining the noise in the imaging process, the problem of minimizing the norm L1 is converted into the problem of solving the projection coefficient theta:

wherein A ═ F_uΨ^Tε is the allowable error

Step four: combining sparsity and K-space data consistency to convert into Lagrange constraint conditions

Wherein A ═ F_uΨ^TThe first term is a data fidelity term and the second term is a regularization term; λ is a regularization parameter that balances the ratio of the two terms;

step five: order to

g(θ)＝||θ||₁Converting the formula (6) into

Step six: modifying theta by adopting a gradient descent method so that the transformation of the formula (7) and the formula (6) is respectively obtained

And

taking into account the continuity of theta to obtain

Considering that the squares of the L1 norm and the L2 norm are separable, problem (8) is transformed into a plurality of minimization problems, and solved by a threshold shrinkage method

Wherein the content of the first and second substances,

shrink(x,β)＝sign(x)·max{|x|-β,0})，

λ＝max(A^Ty)，λ_k+1＝ρλ_k(ii) a ρ is a shrinkage parameter, which is a constant;

step seven: defining an adaptive contraction operator rho A-R rho, and iterating

Up to

Stopping the iteration when a preset value is reached, wherein

Step eight: using θ obtained in steps five to eight, a reconstructed image is obtained using equation (1) in step one.

Has the advantages that: compared with the existing three algorithms, the image acquired by the method can show more local details and clearer outlines; fewer error points are generated, providing more accurate reconstruction; the convergence rate is high, and the iteration efficiency is higher than that of other methods; the image reconstruction of different parts has stability.

Description of the drawings:

fig. 1 is an original image.

Fig. 2 shows the image, detail image and error image after head reconstruction, and the sampling rate is 20%.

Fig. 3 shows a reconstructed image of a blood vessel, a detail image and an error image, and the sampling rate is 20%.

Fig. 4 shows a knee reconstruction image, a detail image and an error image, with a sampling rate of 20%.

FIG. 5 comparison of performance of different methods for human knee and head

FIG. 6 shows the performance comparison of different methods of blood vessels

FIG. 7 is a comparison of performance of different methods of knee joint.

Fig. 8 shows a reconstructed image, detail image and error image of the head, with a sampling rate of 50%.

Fig. 9 is a comparison of performance at different sampling rates.

The specific implementation mode is as follows:

the present invention is further illustrated by the following figures and specific examples, which are to be understood as illustrative only and not as limiting the scope of the invention, which is to be given the full breadth of the appended claims and any and all equivalent modifications thereof which may occur to those skilled in the art upon reading the present specification.

The common model of incomplete MR image data acquisition for solving the sparsest solutions is according to the following formula:

where m ∈ C^NIs a reconstructed image y ∈ C^MIs the K-space pre-sampled signal of the MR_uIs that the observation matrix is expressed by:

F_u＝kF (2)

where k is the sampling pattern and F represents the two-dimensional Fourier transform.

Reconstruction of the original signal can be seen as an L0 norm minimization problem for non-deterministic polynomial time (NP). To address this problem, all possible non-zero values in the sparse coefficients are listed. However, if the observation matrix satisfies the constraint equidistant condition, the MR image can be sparsely represented. The equation transforms the L0 norm minimization problem into the L1 norm minimization problem according to equation (3) below:

where m is sparsely represented by the sparse transform Ψ:

m＝Ψθ (4)

where psi ═ psi₁,ψ₂,...,ψ_N]∈R^N×NHaar wavelet sparse transforms may be employed herein. θ is the projection coefficient of the original image in the sparse domain.

Equation (3) is rewritten as:

where F_uΨ^TIn order to encode the matrix, the encoding matrix,

A＝F_uΨ^T (6)

in combination with the noise contribution from the imaging process, equations (5) and (6) are equivalent to:

where ε is the allowable error, the MRI reconstruction problem includes sparsity and k-space data consistency, given by the following Lagrangian constraints:

in equation (8), the first term is the data fidelity term and the second term is the regularization term. λ is a regularization parameter that balances the ratio of the two terms.

The optimization algorithm is used to solve for θ by equation (8). Then, the reconstructed image may be obtained by equation (1). The threshold shrinking algorithm (SAFISTA) using an adaptive shrink threshold operator is used below for CS-MRI reconstruction. The method for solving equation (8) by SAFISTA is as follows:

order to

g(θ)＝||θ||₁Then, equation (8) can be written as:

consider first the solution of the following equation:

to solve equation (9), during each iteration, a gradient descent method is used to modify θ:

where a is greater than 0, let f (theta) satisfy the Lipschitz-continuous condition,

the same procedure is used to solve equation (8), the iterative equation is as follows:

because θ is continuous, equation (13) is equivalent to:

since the squares of the L1 norm and the L2 norm are separable, the solution of equation 13 is also separated into minimization problems for each element, which can be solved by threshold shrinkage:

here, the shrink soft threshold shrink operator is defined as:

shrink(x,β)＝sign(x)·max{|x|-β,0}) (16)

where sign is a sign function. To increase the convergence speed, parameters t and z are introduced. The specific iteration steps are as follows:

where theta is_k-1And theta_kIs the last two iteration value of θ, equation (15) is rewritten as:

λ is a parameter used to balance the fidelity term f (θ) and the regularization term g (θ). Here, the larger the value of λ, the smaller the specific gravity occupied by g (θ) and the larger the specific gravity occupied by f (θ). Since θ is sparse, the initial value is 0. That is, θ is small in the early stage, and λ should be selected as a large parameter. Typically λ is chosen according to the following method:

λ＝max(A^Ty) (20)

where a is defined by equation (6). And, the nd λ is updated by:

λ_k+1＝ρλ_k (21)

where ρ is a shrinkage parameter, usually taken as a constant, and 0 < ρ < 1. The updated lambda determines each iteration step of the threshold shrink function in equation (19). In order to further increase the convergence speed, an adaptive contraction operator is introduced to increase the initial convergence speed and maintain the final convergence precision. The adaptive shrink operator is defined as follows:

ρ_A＝Rρ (22)

where R is a parameter that we want to take less in early iterations and more in later iterations. According to this rule, the ratio of the values of the total variation of the image of the last two iterations is chosen to define R:

the TV total variation of an image is defined here as:

during the iterative process of MRI image reconstruction, R is initially less than 1, typically around 0.6, and gradually converges to 1 during the iterative process. λ is iterated by equation (25):

the stopping criterion for the iteration is defined as follows:

when ε decreases to a preset value, the iteration stops.

Example 1:

different MR images of various parts of the human body, including MR images of the head, blood vessels and knee joints of the human body, are used as experimental images. The MR images used for the test were acquired on 1.5T Philips Achieva magnetic resonance at national hospital of jiangsu province using 8-channel receiver coils, the gradient echo sequence having the following parameters in table 1.

In this embodiment, q is 2 in equation (25), and in this application, three parameters are used to evaluate the reconstruction quality: mean Square Error (MSE), peak signal-to-noise ratio (PSNR) Structured Similarity (SSIM). The MSE reflects the degree of difference between the estimated value and the original value. PSNR represents the ratio of the maximum possible power of the signal to the noise power. SSIM is a measure of the similarity of two images from brightness, contrast, and structural information. It is more consistent with human visual characteristics than PSNR.

Fig. 1 shows an MR image reconstructed from fully sampled K-space with human head, blood vessel and knee joint data as references.

Figures 2-4 have been constructed using four different methods for reconstructing three parts of a human body. FIGS. 2-4(a) - (d) show reconstructed images using an iterative contraction algorithm with 20% sampling rates of ISTA, FISTA, GTIA, and SAFISTA. FIGS. 2-4(E) - (H) show the same partial details of the reconstructed image by different algorithms. FIGS. 2-4(I) - (L) illustrate reconstruction errors generated by different methods. Because the reconstruction error is small, the contrast of each error image is increased for easy observation.

Fig. 5-7 show three evaluation parameters for each iteration in different ways, with a sampling rate of 20%, after 30 forced iterations. Fig. 5-7(a) - (c) show the performance of MSE, PSNR and SSIM, respectively. Table 2 lists the MSE, PSNR and SSIM of actual MRI data using ISTA, FISTA, GTIA and SAFISTA, with a sampling rate of 20%.

And (4) testing the reconstruction performance of the human head image at other sampling rates. Figure 8 shows a reconstructed human head image at a sampling rate of 50% in a different way. Furthermore, fig. 9 shows the reconstruction performance of MSE, PSNR and SSIM at sampling rates of 30%, 40% and 50%, respectively.

Example the invention was analyzed quantitatively and qualitatively for three evaluation parameters. To verify the effectiveness of the present invention, Iterative Shrinkage Threshold Algorithm (ISTA), fast iterative shrinkage threshold algorithm (SUTA) [ and generalized iterative threshold algorithm (GTIA) were compared.

FIGS. 2-4(a) - (d). Four algorithms are used to reconstruct three parts of the human body. Fig. 2-4(e) - (h) show details of the reconstructed image. The present invention is shown with more detail and a clearer outline than the other three algorithms. Fig. 2-4(i) - (l) show a comparison of reconstruction errors. By contrast, the present invention produces fewer error points and can provide a more accurate reconstruction. Fig. 8 shows a higher sampling rate, and the present invention provides a reconstructed image with richer details.

In an embodiment, the present invention is compared to three other algorithms for iterative performance. Fig. 5-7 show comparative plots of MSE, PSNR and SSIM for each iteration. It can be seen that as the number of iterations increases, MSE gradually decreases, PSNR and SSIM gradually increase. The convergence rate of the invention is obviously superior to that of the other three methods, and the convergence rates are obviously different. These results show that the present invention has a high iteration efficiency.

Table 2 shows that the method has better reconstruction performance in the aspects of MSE, PSNR and SSIM of different parts of a human body. This indicates that the proposed method has good stability for different types of MRI images.

TABLE 1 magnetic resonance imaging parameters

Brief description:

echo time TE

TR repeption time repetition time

FOV field of view moment

Slice of SL: slice

Table 2 reconstructed image evaluation parameter table

Brief description:

MSE mean square error

PSNR peak signal-to-noise ratio

Structural similarity of SSIM

Iterative shrinkage threshold algorithm for ISTA

FISTA (fast iterative shrinkage threshold) algorithm

GTIA generalized threshold iterative algorithm

SAFISTA (self-adaptive fast iterative shrinkage threshold) algorithm

Claims (1)

1. An image processing method based on an adaptive fast iterative shrinkage threshold algorithm is characterized in that: the method comprises the following steps:

the method comprises the following steps: acquiring K-space undersampled MR image data, and reconstructing a model as follows:

F_u＝kF (2)

wherein m is C^NIs a reconstructed image, y ∈ C^MIs the K-space undersampled signal of MR, F_uIs an observation matrix, k is a sampling pattern, and F represents a two-dimensional fourier transform;

step two: introducing a projection coefficient theta of an original image in a sparse domain, and carrying out sparse representation on m to obtain m-psi theta, psi-phi₁,ψ₂,...,ψ_N]∈R^N×NConverting the L0 norm minimization problem of the model into a norm L1 minimization problem;

m＝Ψθ (4)

step three: combining the noise in the imaging process, the problem of minimizing the norm L1 is converted into the problem of solving the projection coefficient theta:

wherein A ═ F_uΨ^TWhere ε is the allowable error;

step four: combining sparsity and K-space data consistency to convert into Lagrange constraint conditions

Wherein A ═ F_uΨ^TThe first term is a data fidelity term and the second term is a regularization term; λ is a regularization parameter that balances the ratio of the two terms;

step five: order to

Converting formula (6) to

Step six: modifying theta by adopting a gradient descent method so that the transformation of the formula (7) and the formula (6) is respectively obtained

And

considering the continuity of θ, the following equation is obtained:

considering that the squares of the L1 norm and the L2 norm are separable, problem (8) is transformed into a plurality of minimization problems, and solved by a threshold shrinkage method

Wherein the content of the first and second substances,

shrink(x,β)＝sign(x)·max{|x|-β,0})，

λ＝max(A^Ty)，λ_k+1＝ρλ_k(ii) a ρ is a shrinkage parameter, which is a constant;

step seven: defining an adaptive shrink operator ρ_AIterating for R ρ

Up to

Stopping the iteration when a preset value is reached, wherein

Step eight: using θ obtained in steps five to eight, a reconstructed image is obtained using equation (1) in step one.

Images