Sparse reconstruction of magnetic resonance image combined with two-step iteration and adaptive shrinkage factor

1.   Introduction

Nyquist sampling theorem is the criterion that must be followed in traditional image acquisition system. The traditional acquisition mode often acquires raw image data at a very high sampling frequency, it may take a long time and require high equipment performance. The same, traditional magnetic resonance (MR) imaging also takes too long time, which is easy to cause motion artifacts, and patient compliance is poor. In addition, for imaging that requires injection of contrast agents, prolonged sampling can result in reduced contrast, which can result in poor quality or undiagnosable images.

In order to solve the above problems, Compressed Sensing (CS) theory was proposed ^[1], based on the theory of sparse decomposition and signal approximation. In this theory, the raw data which are undersampled can be reconstructed by performing nonlinear optimisation, and without reducing the image quality. In 2007, CS was applied to the field of MR imaging and achieved amazing results ^[2]. After that, efficient reconstruction algorithms had always been a research hotspot in undersampling MR image reconstruction field. Among them, Iterative shrinkage thresholding algorithm (ISTA) ^[3], Split Bregman Iteration ^[4] and conjugate gradient (CG) ^[5] were the mostly used algorithms. Compared with Bregman and CG, ISTA had lower algorithm complexity and directly solved the L1 minimization problem to get faster iterations ^[3,6,7,8,9]. In addition, with the deepening of deep learning in medical image processing research, it also provides new solutions for fast magnetic resonance reconstruction ^[10,11].

Therefore, a large number of people have studied the ISTA algorithm to improve the reconstruction accuracy and improve the iterative efficiency, and a series of improved ISTA were brought up. In 2007, Bioucas-Dias and Figueiredo ^[12] proposed Two-step Iterative Shrinkage Thresholding Algorithm (TwIST), the first two iteration estimates and threshold functions are used to estimate the current value at each iteration. In 2009, Beck et al. ^[13] proposed Iterative shrinkage thresholding algorithm (FISTA) based on the results of the first two iteration to improve global rate of convergence. In 2015, Zhang et al. ^[14] proposed an algorithm used exponential wavelet as sparse transformation, it can improve MR image reconstruction effect by increasing signal sparsity. In 2015, Li et al. ^[15] proposed Hessian Schatten norm-based regularization instead of the total variation (TV), it can keep the favorable properties of TV and reduce the staircase effect appearing in TV-based reconstructions. In 2016, Wu et al. ^[16] used regularization parameter to be taken as a threshold in a fixed-point ISTA. In 2019, Shang et al. ^[17] proposed an easy-to-implement algorithm based on the framework of alternative direction method, named iterative p-shrinkage thresholding algorithm (PISTA), for solving the low Tucker rank tensor recovery problem. And it has been successfully applied in MRI image recovery.

However, in order to improve the efficiency of the iteration, early in the iterations, the weight of the fidelity term should be increased to speed up the convergence. Late in the iterations, the weight of the regular term should be increased to improve the convergence accuracy. But in above algorithms, the weight of fidelity item and regular item are fixed. In this paper, two-step iterative Thresholding Algorithm with dynamic shrinkage factor (DTwIST) was proposed, DTwIST can dynamic adjust shrinkage factor to rebalance the weight of fidelity item and regular item to improve the convergence efficiency, and can improve the final convergence accuracy.

2.   Materials and methods

2.1.   Magnetic resonance imaging model under compressed sensing

To apply the theory of compressed sensing to rapid imaging of magnetic resonance, the first is to use part of k-space data obtained by undersampling as a known condition, and then use the a priori knowledge of the sparseness of the magnetic resonance image in the specific transform domain, and then construct the convex optimization problem of L1 norm in combination with the compressed sensing theory, at last an optimal solution for the magnetic resonance image can be got. The basic mathematical model is shown as follows:

where $u$ is the magnetic resonance image to be reconstructed, $\Psi$ stands for linear sparse transform, it can be total variation (TV), wavelet transform, curvelet transform or other sparse transform, $F_{u}$ stands for k-space sampling mode, $y$ is the k-space data obtained by undersampling measurement, $\mathcal{E}$ represents the noise level of the measured data. Because k-space data is a two-dimensional Fourier transform of magnetic resonance images, the measurement matrix in the model is the undersampled Fourier transform matrix. Convert Eq (1) into unconstrained optimization problem form by lagrange multiplier method:

The regularization parameter $\alpha$ in Eq (2) is used to balance the proportion of the fidelity term $\left\|F_{u} u-y\right\|_{2}^{2}$ with the regular term $\left\|\Psi_{u}\right\|_{1}$. The sparse transform domain coefficient of MR image $\boldsymbol{u} u$ is $m$, which defines as follows:

$\Psi$ is an orthogonal matrix, and thus its transpose is equal to its inverse, Eq (2) can be rewritten as follows:

Reconstructed magnetic resonance images can be obtained by solving Eq (4).

2.2.   Reweighted iterative shrinkage thresholding algorithm

RISTA mainly uses the first two reconstructed results to adjust the weights $\alpha$ of the fidelity and regular item. The specific process of RISTA for solving Eq (2) is as blow.

Therefore solving Eq (4) is equivalent to solving the following Eq:

Suppose $f(m)$ satisfies the condition of being Lipschitz continuous. Then, it can be approximated as follows near previous iteration $m_{k-1}$.

where $L$ is the lower bound of the derivative of $f(m)$, called Lipschitz constant, and $k$ is the current number of iterations. Combining the gradient descent method, a more general quadratic approximation model is used to solve the optimization problem of Eq (5).

After ignoring the constant term $m_{k-1}$ and $\nabla f({m_{k - 1}})$, the above Eq (7) can be changed into the following:

Based on separable characteristics of L1 norm and the square of the L2 norm, Eq (8) can be solved as following 3:

where $C = F_{\mathrm{u}} \Psi^{T}$, the function $G_{k}(x)$ is defined as Eq (10):

where $\operatorname{sgn}(x)$ is a symbolic function, $\alpha$ is a constant and determines the basic shrinkage step size for each iteration thresholding.

2.3.   Two-step iterative thresholding algorithm with dynamic shrinkage factor

Here we introduced a dynamic shrinkage factor $\mu_{{k}}$ in Eq (10) to dynamically adjust $\alpha$. Then Eq (10) can be rewritten as:

The update function of the regularization parameter $\mu_{k}$ is updated as follows:

where $m_{k-1}$ and $m_{k-2}$ are the results of the previous two iterations. As the number of iterations increases, $\left\|m_{k-1}-m_{k-2}\right\|_{2}$ tends to zero and $\mu_{k}$ tends to 1. By way of adjusting the regularization parameters $\alpha$, the weight of fidelity items is larger in the early stage of iterations, and it is conducive to rapid convergence. But in the late iterations, the weight of regular items is larger, so that it can ensure the image reconstruction accuracy. Parameter $\sigma$ is used to adjust the change rate of the regularization parameter.

Then combined with method of iterative reweighted shrinkage ^[13], we used the first two iteration values to adjust the current iteration value:

where $m_{k-2}$, $m_{k-1}$ and $m_{{k}}$ are the iterative values of last three times, we use Eq.(13) to adjust the current iteration value. Where $\gamma$ and $\beta$ are calculated as follows:

where $\lambda_{1}$ and $\lambda_{2}$ are chosen as reference ^[10], here we set them as $\lambda_{1} = 0.001, \lambda_{2} = 1$.

Next the iterative step size of the algorithm needs to be considered. In paper ^[18,19,20], the convergence of ISTA algorithm was proved and the initial value $\alpha$ could be chosen as below 21:

The termination condition of the image iteration is given by the following:

where $\varepsilon$ is the maximum allowable error set in advance.

2.4.   Algorithm implementation

The sparse transform maps the original image into the sparse domain, making the non-zero coefficients of the sparse domain more sparse. Generally different sparse transforms are available for different types of images. The requirement of the construction of the sparse transform is to make the transform coefficients of the image on the set of bases more sparse by finding a set of orthogonal bases. In the ISTA series of algorithms, two-dimensional discrete wavelet transform (DWT2) were widely used, this sparse transform was also used in this paper.

In CS-MRI, since the sparse transform is fixed, the incoherent undersampling can only be achieved by designing the measurement matrix to ensure that the RIP condition is satisfied ^[22,23]. In two-dimensional magnetic resonance imaging, the measurement matrix is embodied in the acquisition mode of the K space. The variable density spiral trajectories 24 was used in this experiment.

The specific algorithm implementation steps are as follows:

Algorithm 1. Two-step iterative Thresholding Algorithm with dynamic shrinkage factor

Inputs:

$y$ : undersampled observation value in k-space

$\psi$ : discrete wavelet transform

Initialization:

Repeat

Update $\mu_{k}$ by Eq (12)

Update $G_{k}(x)$ by Eq (11)

Update $m_{{k}}$ by Eq (9)

Update $\gamma_{1} \beta$ by Eq (14)

Renew $m_{{k}}$ by Eq (13)

$k = k + 1;$

Until $E < t o l$

Output: reconstructed MR image $u_{k} = \Psi^{F} m_{k}$

3.   Simulation setting

As shown in Figure 1, MR images of the abdomen and lumbar spine were selected as experimental images to test the reconstruction performance of CS-MR images. The MRI raw data came from the Department of Radiology, Jiangsu Provincial People's Hospital, and they were acquired on a 1.5T Philips Achieva scanner with a 8-channel receiver coil.

The abdomen image matrix size = 448 × 448 and the lumbar spine image matrix size = 528 × 528. To control the effect of different matrix sizes, the images are resampled before algorithm processing, and the matrix size after the two images are resampled = 512 × 512.

The sampling pattern of variable density spiral trajectories were set to 20, 30 and 40%, and they represented 5×, 3.3× and 2.5× accelerations. Spiral is a circular spiral, radial curve. It starts from the center and does not collect corner data. The spatial resolution of the image is isotropic 20% spiral sampling trajectory are shown in Figure 2.

Two general stopping criterias were used: maximum number of iterations (50) and tolerance (10^-5), and when either criterion is satisfied, the reconstruction stops.

ISTA, PISTA and TwIST were chosen as comparison algorithms for proposed method. The reconstruction effects were compared. 4 evaluation indicators were selected: the Peak Signal-to-Noise Ratio (PSNR in dB), Structural Similarity Index (SSIM) 25, Transferred Edge Information (TEI) 26 and normalized mutual information (NMI).

here, $Q_g^{\tilde x\hat x}$ and $Q_\alpha ^{\tilde x\hat x}$ represent the edge strength and direction preservation values respectively. The larger the value of TEI is, the better the edge information of the reconstructed image is kept.

where $H(A)$ and $H(B)$ are the entropies of images A and B, respectively, and $H(A, B)$ is the joint entropy of images A and B. The larger the NMI value, the higher the similarity between the two images.

4.   Results

In order to fully evaluate the performance of the proposed algorithm (DTwIST) in this paper, we have designed several sets of experiments, and ISTA, PISTA, TwIST were used as comparison algorithms.MATLAB (R2016a) is used to implement the algorithm.

First, we compare the reconstruction performances of different images at the same sampling rate (20%) by 4 algorithms. Figure 3 shows the images were reconstructed by 4 different algorithms at 20% sampling rate after 50 iterations. 4 columns of images represent ISTA, PISTA, TwIST and DTwIST. The results show that the reconstructed images with DTwIST can get sharper edges and better local detail than other 3 at 20% sampling rate.

Figures 4 and 5 show the convergence curves of the reconstruction performance. PSNR, SSIM, TEI and NMI were used to quantitatively evaluate the performance of the 4 algorithms. The abscissa is the number of iteration. The results show that the image reconstructed based on DTwIST have a faster convergence rate and better reconstruction performance than the other three algorithms during the iterative process.

Table 1 shows the numerical comparison of evaluation parameters- iterations, PSNR, SSIM, TEI and NMI reconstructed by the ISTA, PISTA, TwIST and DTwIST algorithms when the convergence iteration is stopped. The results show the superiority of DTwIST, especially excellent convergence speed.

In order to verify the applicability of the algorithm at other sampling rates. We also did experiments at 30 and 40% sampling rate.

Figures 6 and 7 show the reconstruction performance at different sampling rates after 50 iterations. PSNR, SSIM, TEI and NMI were used to quantitatively evaluate the performance of the 4 algorithms. The abscissa represents the sampling rate. The results show that the proposed algorithm has similar performance at different sampling rates.

In order to further verify the performance of the proposed algorithm, we selected a total of 60 MR images from different slices of the brain, ankle and knee joints for further experiments.

Figures 8-10 show the box plots of the reconstruction performance of 60 MR images. The sampling rate is 20, 30, 40% respectively. The results show that the proposed algorithm has better reconstruction performance for MR images of other body parts.

5.   Discussion and conclusions

The proposed algorithm (DTwIST) overcomes the disadvantage that the TwIST algorithm is difficult to balance speed and efficiency. By dynamically adjusting the shrinkage factor, it can greatly accelerate the iterative convergence while ensuring the convergence accuracy.

By compared with ISTA, PISTA, and TwIST in terms of PSNR, SSIM, TEI, NMI and number of iterations, a conclusion can be drawn than the proposed methods can retain more global information and local information of reconstructed images, especially with excellent convergence speed. At the same time, it provides a reference for the improvement of similar algorithms.

Acknowledgments

This work was supported in part by the Jiangsu Provincial Key Research and Development Program under Grant BE2020714, the Postgraduate Research and Practice Innovation Program of Jiangsu Province under Grant KYCX21_1556. Ethical Approval: The Research Ethics Committee of the Nanjing Medical University granted approval for this study and the collection and use of the data analyzed in this study. Approval number: NMUE2021301.

Conflict of interest

The authors declare there is no conflict of interest.