Iterative CT reconstruction based on ADMM using shearlet sparse regularization

1.   Introduction

To improve spatial resolution and reduce image artifacts and radiation dose in CT scanning ^[1,2,3,4], algorithms for iterative reconstruction have been proposed and developed. Regularization is an important type of iterative reconstruction algorithm that can create high-quality reconstructed images ^[5]. With the emergence of the compressed sensing theory ^[6], a regularization technique based on Total Variation (TV) minimization has demonstrated excellent performance for CT imaging with sparse-view scanning and interior tomography ^[7].

TV-based methods favor solutions with the hypothesis of the piece-wise constant of the expected image from few-view samplings, such as simple geometric images with flat intensity. For example, an innovative algorithm was formulated by Sidky et al. ^[8] using Projection Onto Convex Sets (POCS) with TV minimization. It has been reported that an algorithm with TV minimization can reconstruct normal-quality CT images from few-view projections without importing obvious artifacts ^[9]. A Penalized Weighted Least-Squares (PWLS) algorithm has been shown to maintain high-quality images with the regularization of Total Generalized Variation (TGV) ^[10]. The PWLS-TGV algorithm yields a more obvious increase in the spatial resolution and accuracy characteristics of the reconstructed images than algorithms based on TV regularization. Although TV minimization can reduce streaking artifacts, the patchy feature can be introduced into the images, potentially affecting medical diagnosis. Some experimental results indicate that at least 100 views are required for the realistic reconstruction of the human anatomy; therefore, TV minimization is unsuited to very low dose levels in CT reconstruction. Consequently, the key is to produce regularizers that offer an optimal sparse representation of medical images as compared with TV; wavelet transform is one of these sparse transforms. In comparison with the loss of soft edges following TV minimization, wavelet transformation permits the regularization to change the image using resolution scales. Haar wavelets are very similar to TV for piece-wise constant images. Garduno et al. ^[11] investigated the simplest possible wavelet, the Haar wavelet, for regularized CT reconstruction. Reconstruction with regularization using Haar transform is not more successful for CT reconstruction than that using the small TV value. Steidl et al. ^[12] found that for one-dimensional signals, Haar wavelet transform and TV regularization with a one-step spatial discrete are comparable for two pixel-pairs when a single scale is put into effect; therefore, regularization based on Haar transform does not exceed that with TV. Some regularization of ridgelets, shearlets ^[13] and curvelets ^[14] based on sparse wavelets and dictionaries of edge-preservation have been developed to overcome the limitation that images reconstructed by TV regularization may lose their contrast and fine structures and staircase artifacts may be produced. Curvelet transform demonstrates good performance in discontinuities along curves with bounded curvatures, while ridgelet transform demonstrates good performance in discontinuities along straight lines.

There are a lot of benefits to using shearlets in image reconstruction ^[15]. Firstly, a lower redundant sparse tight frame representation can be allowed by shearlets as compared with that allowed by other related multiresolution representations; moreover, shearlets also offer shift invariance and analysis of direction. The basic functions of shearlet transform are matched with the streak artifacts of X-ray noise via directional analysis in the reconstructed images. The streak artifacts can then be reduced, which permits the stripes to be approached by unimportant coefficients. Secondly, a series of multi-scale spaces can be generated, into which Shearlet transform can decompose the $L^{2}\left(\mathbb{R}^{2}\right)$ of images. While the image edges are expressed quantitatively within the smaller subclass of functions of bounded variation (BV), other noticeable features from homogeneous regions to oscillatory patterns of noise or texture are included in the intermediate ranges from the small $B V\left(\mathbb{R}^{2}\right)$ to the large $L^{2}\left(\mathbb{R}^{2}\right)$.

The Alternating Direction Method of Multipliers (ADMM) was proposed for the design of a CT reconstruction algorithm to solve the problem of PWLS ^[16]. The split Bregman and ADMM are two important methods based on primal-dual splitting. When applied to X-ray CT reconstruction, the split Bregman method achieves the optimization solution by minimizing the cost function involving the fidelity of raw statistics. To solve a generic PWLS model for Low-dose CT image reconstruction, He et al. ^[17] parameterize a plug-and-play ADMM method. Then, the fundamental principle of Deep Learning is applied to optimize the parameters of the parameterized plug-and-play ADMM (3pADMM). The ADMM is used to create the suggested algorithm, sparsity-constrained angiography (SCAN), with a focus on the sparse character of the angiography imaging ^[18].

In the present paper, we designed a compressed sensing CT iterative reconstruction algorithm with Shearlet Sparse Regularization (SSR) based on ADMM minimization and verified the algorithm using simulation experiments and real data. The advantages and disadvantages of this method were expressed through numerical phantom and clinical datasets. Lung CT reconstruction experiments show that for lung CT images with rich texture details, Shearlet sparse regularization reconstruction images show an excellent display of texture details. For comparison, reconstruction images based on FBP and TV regularization are shown. In the experiments using real data, projections of a spiral CT were used to evaluate whether this method applies to this application.

2.   Methods

The model of CT scanning and imaging is represented using a discrete linear system:

where the matrix A is the system matrix representing X-ray transformation; the vector y represents the projection data; and the vector x denotes the image data; and $\mathrm{\eta }$ denotes the noise measuring the projection data. Generally, the number of pixels in a reconstructed image is greater than that of projecting samples in CT imaging; thus, Eq (1) is ill-posed. The goal of regularization is to introduce further regularity to the problem to make the solution unique and stable. Subsequently, for noiseless and under-determined systems, we applied the method of least-squares to minimize the cost function to solve x:

Usually, we can solve the optimization problem (1) using the gradient descent algorithm. In practical CT imaging, the noise in the domain of projection makes minimizing (1) for x difficult, and even more difficult in sparse projection systems. Thus, we built the optimization problem based on the theory of Compress Sensing (CS). In CS, if the sparse transform vector $\varPsi$ of an image x is given, image x can be accurately reconstructed with a higher possibility using fewer samples than demanded by the Nyquist criterion, since the elements of the vector $\varPsi_x$ are mainly zeros. The l₁ norm of the vector $\varPsi_x$ can express this sparsity. As a result, under the guidance of CS theory, we generated the optimization constrained minimization problem for CT reconstruction:

where ε represents the error between the measurement y and the projection of the image x, and $\varPsi$ is the shearlet transform in our study, which is a linear sparse transform.

For shearlet regularization, we gave $\varPsi = SH$, using SH to imply the shearlet transform. For a two-dimensional space, the discrete shearlet system was expressed as follows:

where A is the anisotropic expanding matrix and S is the shear matrix. The discrete shearlet transform was defined as follows:

In the present study, the shearlet transform was implemented using the Matlab-coded shearlet toolkit.

To solve Eq (3), the augmented Lagrangian function corresponding to the unconstrained optimization problem is:

where the parameter $\rho > 0$ is the penalty parameter and $\mu$ denotes a Lagrange multiplier. We can obtain the convenient formula of (7) by combining the last two terms inside the quadratic term as follows:

Three subsequent steps finish one ADMM iteration:

in this particular case, it produces:

By solving the variables x, z, μ alternately, an exact solution was obtained after several iterations, and the x update was solved using 10 iterations of the steepest descent method. ADMM still converges even for this x.

For comparison, we simultaneously developed the TV regularization of anisotropy based on ADMM. For anisotropic TV, we set the sparse transform T = TV, where TV remains the discrete differential operator. For the other comparison method, we implemented the algorithm of the simultaneous algebraic reconstruction techniques (SART). In addition, this method uses a relaxing factor λ ∈ (0, 2).

3.   Experiments

3.1.   Data simulation

Firstly, we studied the FORBILD ^[19] high contrast resolution numerical phantom developed by Siemens Healthcare. The high contrast numerical phantom is a proposed algorithm for measuring performance and nature. The true image size of the phantom is 256 × 256 and the pixel value range is [0, 1]. The projection data of the phantom were simulated using Matlab's fanbeam routine with 360 projections over 360º. We added white Gaussian noise to the projection data using the standard deviation of 1% after forwarding projection.

3.2.   Experimental data acquisition

Due to the advantages of the algorithm SSR, we designed three experiments using a Neusoft NeuViz 64 X-ray CT scanning system to measure its performance and features. The CTP528 high contrast phantom contains 21 line pairs along the radius, with diameters ranging from 0.24 to 5 mm. The last dataset was projections of a human lung containing textural features. Data were collected using a diagnostic scanner (NeuViz 64 X-ray CT, Neusoft Chinese Medical Systems, Shenyang, Liaoning, China). The scanning radius was 570 mm; the detector had 64 rows; the row height was 0.625 mm; each row had 672 channels, and the fan angle was 52º. There were 1160 uniform scanning views over 360º. The size of the reconstructed image was 256 × 256. For the three scanning experiments, the tube voltages were 120 and 140 kV, and the tube currents were 300 and 150 mA.

3.3.   Parameter setting

The ADMM algorithm includes two parameters, the regularization parameter λ and the coupling parameter ρ. In the case of the numerical phantom, we evaluated the reconstructed images using the mean squared error (MSE) since we had the ground truth. For the real datasets, there were no true images; consequently, the metric Q could be used. A non-reference metric is proposed by the metric Q, which provides a quantitative measure of the image quality and content using single value decomposition for local image gradient. Zhu and Milanfar ^[20] developed a tool for the metric Q to select parameters in image processing and showed that it displays excellent performance in image denoising. In Zhu and Milanfar's method, the image $I(N\times N)$ is divided into disjoint patches of $M$ $\{{P}_{i}{\}}_{i\in \{0, ..., M\}}$. The gradient matrix ${G}_{i}$ for each patch is computed. The image content metric ${Q}_{i}$ is defined as follows:

where ${s}_{i, 1}$ and ${s}_{i, 2}$ represent single values of the gradient matrix ${G}_{i}$. In addition, formula (14) satisfies ${s}_{i, 1} > {s}_{i, 2}$, wherein the coherence ${R}_{i}$ is also defined as:

In the end, the final metric Q is computed as follows:

where $K\subset [1, M]$ is the index set of patches that are more than the thresholding. We used $8\times 8$ pixel patch given in the reference. Additionally, we also computed a contrast to noise ratio (CNR) and compared the reconstruction images with those based on shearlets and TV. The CNR formula that we adapted is defined as follows ^[21]:

where ${\overline{f}}_{S}$ and ${\sigma }_{S}$ denote the mean and standard deviation within a region of interest (ROI), and ${\overline{f}}_{b}$ and ${\sigma }_{b}$ denote these values of the ROI over the background region.

3.4.   Implementation details

All datasets were reconstructed by SART, ADMM using isotropic TV regularization (ADMM-TV), and ADMM using sparse shearlet regularization (SSR) in MATLAB. All methods were executed on a workstation in which the CPU was an Intel Xeon E5-2623. The ADMM was initialized with a zero matrix.

4.   Results and discussion

4.1.   Data simulation

The first experiment for evaluation was designed to reconstruct the Forbild high contrast simulation phantom, for which ground truth data were available as shown in Figure 1. Figure 1(a)−(c) shows the full image reconstruction for SSR, ADMM-TV and SART. Zoomed images of the region of interest (ROI), marked with a red square, are shown in Figure 1(e) for the ground truth; Figure 1(f) for SSR; and Figure 1(g) for ADMM-TV and Figure 1(h) for SART. Figure 1 also shows the CNR of the three reconstructed images. Figure 2(a), (b) shows the profile of marked red lines on the images reconstructed using the three methods and the true image. We selected the coupling parameter ρ and the regularization parameter λ according to the parameter of MSE because the ground truth image was available. For the Forbild simulation phantom dataset, we computed the MSE in the range of $\rho \in\left\{10^{0}, 10^{1}, \ldots, 10^{5}\right\}$ and $\lambda \in\{0, 10, 20, \ldots, 200\}$. Figure 2(b) shows the MSE profiles for a series of values for the regularization parameter λ based on SSR and TV, and we chose coupling parameter ρ in terms of the MSE minimum value. We selected the optimal parameters ρ = 10, λ = 70 for SSR and ρ = 10, λ = 40 for TV. The plots of metric Q over the regularization parameter λ for SSR and ADMM-TV are also shown in Figure 2(b). In comparison with the plots of MSE, the optimal regularization parameter λ varied by 20 for SSR, while it varied by 30 for ADMM-TV. For the selection of the iteration i, the metric Q search was performed in the range $i \in\{1, 2, 3, \ldots, 50\}$. Figure 2(b) shows the profiles of metric Q for iteration i for SSR, ADMM-TV, and SART. The chosen optimal iterations according to metric Q for the Forbild simulation phantom were i = 9 for SSR and ADMM-TV, and i = 500 for SART.

We used the FORBILD high contrast simulation phantom dataset to show the advantages of the shearlet regularization method. Figure 1(e)−(h) indicates that in the same window, for SSR method, the gray value of the holes is small and the other parts of the image are smoother than the results of the ADMM-TV method. For the SART method, although the gray values of the holes are lower than those for the SSR method, the noise of the reconstruction image is too high. We selected holes of 1.3 mm, which is the location of the line profile (Figure 2(a)), to compare the results of the three reconstruction methods. Figure 2(a) indicates that SSR reconstructed the Forbild high contrast features more precisely than ADMM-TV, while the results of SART are better than those of SSR. For the piece-wise constant region (the middle circle), the line profile of which is shown in Figure 2(b), the variation in the values on the line profile for the SSR method was lower than that for ADMM-TV and SART, and the variation of the SART method was large. The value of CNR for SSR was minimum, followed by ADMM-TV and SART. Figure 2(b) shows that the minimum MSE for SSR was less than that for ADMM-TV, and the maximum metric Q for SSR was greater than that for ADMM-TV.

4.2.   The high contrast phantom dataset

For the high contrast phantom dataset, the images reconstructed using the three algorithms are shown in Figure 3(a)−(c). The red rectangle range is an ROI, which is zoomed and indicated in Figure 3(d)−(f) with a red line marking the position of the line profiles. In the reconstructed images, the high contrast resolution was 16, 18 and 20 cm/lp followed by the line marked in red, as shown in Figure 3(a)−(c) from top to bottom. Additionally, the CNRs of all images are shown. Figure 3(h) shows the curve of the pixel value followed by the line marked in red using different reconstruction algorithms. This time, since there was no ground truth, we selected ADMM parameters using the metric Q for the range $\rho \in\left\{10^{0}, 10^{1}, \ldots, 10^{5}\right\}$ and $\lambda \in\{0, 10, 20, \ldots, 400\}$. Figure 3(g) shows the results of plotting the metric Q curve according to the $\lambda$ with SSR and ADMM-TV, while the coupling parameter $\rho$ was chosen according to the maximal value over all reconstructions. For the SSR method, we selected the optimal parameters with ρ = 1, λ = 100; for the ADMM-TV method, we selected the optimal parameters with ρ = 10, λ = 160. The plots of the metric Q over the iteration i varying in the range $i \in\{0, 1, 2, \ldots, 50\}$ are shown in Figure 3(c). The chosen optimal iterations were i = 9 in the case of SSR, i = 30 in the case of ADMM-TV, and i = 500 in the case of SART.

The high contrast phantom is usually used to evaluate the spatial resolution for spiral CT. The zoomed ROI in Figure 3(d)−(f) containing three spatial resolutions, 16, 18 and 20 cm/lp, shows that the three-line pairs for the SSR method were more apparent than those for the ADMM-TV and SART methods, while the other parts of the ROI for the ADMM-TV method were smoother than those for the SSR and SART methods. In addition, the CNR value for the TV method was larger than that for the SSR and SART methods.

4.3.   The lung dataset

For the lung CT experiments, Figure 4(a)−(c) shows the images reconstructed using the three methods. We marked the three ROIs with red rectangles (right, middle, and bottom) and used a red line for the line profiles. The labels in each figure indicate the CNR, and the zoomed ROIs are shown in Figure 4(d)–(f) for the three methods. Additionally, Figure 4(h) shows the line profile of the red marks. We selected parameters using the metric Q as a metric. ρ was varied in the range $\rho \in\{0, 10, \ldots, 500\}$, while λ was varied in the range $\lambda \in\{0, 1, 2, \ldots, 50\}$. Figure 4(f) shows the metric Q plots. The optimal parameters for SSR were ρ = 140, λ = 26, while those for ADMM-TV were ρ = 220, λ = 24 according to the maximum metric Q value. The iteration parameters were varied in the range $i \in\{0, 1, 2, \ldots, 100\}$. The optimal iterations were i = 100 for SSR, i = 29 for ADMM-TV, and i = 500 for SART according to the metric Q plot shown in Figure 4(g).

In the lung CT dataset (see Figure 4), the reconstructed images for SSR and ADMM-TV have similar edge characteristics; however, when the ADMM parameters were selected by the matrix Q for ADMM-TV, the noise increased significantly. This is because the metric Q selected the low iterations. The CNR values for the three reconstruction methods also indicate the characteristics. For the zoomed ROI images (see Figure 4(a)−(c) and (d)−(e)), it is again obvious that SSR has fewer block images as compared with ADMM-TV. Simultaneously, around the bone in Figure 5(d)−(e), SSR better preserves the features of the edges, while ADMM-TV introduces the familiar staircase artifacts. However, we found that SSR can strengthen the linear artifacts in all directions, as indicated by the arrows in Figure 4(h)−(j), while ADMM-TV can reduce these artifacts. This is due to the characteristics of shearlet transform, which provide optimally sparse representations for all directionality controlled by shear matrices. Normally, the value of the sparse image is set to zero through each iteration of the soft thresholding step of ADMM, making it below the threshold. The ADMM-TV method eliminates noise if the gradients of that noise are below the threshold. As a result, ADMM-TV sets the low gradients to zero. On the contrary, the coefficients of the noise characteristic edges are few because of the directionality of the shear frame. Therefore, SSR can strongly eliminate noise while preserving elongated characteristics and texture. The reconstructed image used the SSR method to better retain the edges than the ADMM-TV method, especially in the case of high noise (e.g., low-dose scanning), when we scanned phantoms that contained some information regarding edges or elongated features. The regularization parameter λ and the coupling parameter ρ involved in the ADMM algorithm can significantly impact the reconstruction properties; thus, parameter selection plays an equally important role. For real datasets, we employed the metric Q because there was no true image. To verify the accuracy of the metric Q, we drew the curves of the MSE and the metric Q in the same figure to compare the results of parameter selection (see Figure 4(g)). Although the value of the optimal parameter using MSE was not the same as that using the metric Q, these two quality measures were very similar.

In the present paper, the shearlet is limited to two dimensions. The expansion of three-dimensional iterative CT reconstruction is required in future work. Although the shearlet method can be directly extended to three dimensions, its calculation is relatively large. Future work should include the acceleration of the SSR method using GPU and the application of these developed techniques to different datasets to compare different regularization strategies.

5.   Conclusions

We propose an ADMM-based algorithm to perform iterative CT reconstruction using shearlet regularization. We investigated SSR and compared its results with those of ADMM-TV and SART on three phantoms: the FORBILD simulation phantom, the real high contrast phantom, and lung datasets. The results demonstrate that the SSR method exhibits the advantages of high directionality and contrast as compared with ADMM-TV. On the contrary, the ADMM-TV method had a better reconstruction effect in the piece-wise constant region.

Acknowledgements

This work was supported by "the Fundamental Research Funds for the Central Universities" (Project number: 2021MS089) and "Natural Science Foundation of Liaoning Province of China" (Project number: N2119005).

Conflict of interest

The authors declare there is no conflict of interest.