Convolutional Sparse Coding for Compressed Sensing CT Reconstruction

TL;DR

This paper introduces a convolutional sparse coding approach for sparse-view CT reconstruction that processes entire images directly, leading to improved detail preservation and artifact reduction compared to patch-based methods.

Contribution

The paper proposes a novel CSC-based method for CT reconstruction that avoids patch division, maintaining image consistency and detail, and demonstrates superior performance over existing methods.

Findings

Outperforms state-of-the-art methods in simulated data

Effective on real CT data with improved image quality

Reduces artifacts and preserves details in reconstructed images

Abstract

Over the past few years, dictionary learning (DL)-based methods have been successfully used in various image reconstruction problems. However, traditional DL-based computed tomography (CT) reconstruction methods are patch-based and ignore the consistency of pixels in overlapped patches. In addition, the features learned by these methods always contain shifted versions of the same features. In recent years, convolutional sparse coding (CSC) has been developed to address these problems. In this paper, inspired by several successful applications of CSC in the field of signal processing, we explore the potential of CSC in sparse-view CT reconstruction. By directly working on the whole image, without the necessity of dividing the image into overlapped patches in DL-based methods, the proposed methods can maintain more details and avoid artifacts caused by patch aggregation. With predetermined filters, an alternating scheme is developed to optimize the objective function. Extensive experiments with simulated and real CT data were performed to validate the effectiveness of the proposed methods. Qualitative and quantitative results demonstrate that the proposed methods achieve better performance than several existing state-of-the-art methods.

Tables

Table 1. TABLE I: Quantitative results obtained by different algorithms for the abdominal image.

Num. of Views		48			64			80
	PSNR	RMSE	SSIM	PSNR	RMSE	SSIM	PSNR	RMSE	SSIM
FBP	24.15	0.06203	0.41375	25.35	0.05417	0.45858	26.40	0.04785	0.54287
PWLS-TGV	41.06	0.00885	0.96068	44.05	0.00627	0.97697	46.91	0.00451	0.98736
PWLS-DL	41.75	0.00817	0.96863	44.84	0.00573	0.97975	47.65	0.00414	0.98915
PWLS-ULTRA	41.69	0.00824	0.96672	44.46	0.00599	0.97791	46.02	0.00500	0.98355
LEARN	40.43	0.00953	0.95421	43.28	0.00692	0.97254	44.74	0.00580	0.97945
PWLS-CSC	41.36	0.00855	0.96375	44.92	0.00568	0.98064	47.75	0.00410	0.98911
PWLS-CSCGR	42.52	0.00748	0.96919	45.73	0.00517	0.98329	48.38	0.00381	0.99034

Table 2. TABLE II: Quantitative results obtained by different algorithms for the thoracic image.

Num. of Views		48			64			80
	PSNR	RMSE	SSIM	PSNR	RMSE	SSIM	PSNR	RMSE	SSIM
FBP	21.85	0.08079	0.36232	22.99	0.07087	0.42335	23.9	0.0638	0.47914
PWLS-TGV	40.74	0.00919	0.96781	44.43	0.00600	0.98310	48.44	0.00379	0.99231
PWLS-DL	39.80	0.01023	0.96203	43.68	0.00655	0.98145	47.21	0.00436	0.99069
PWLS-ULTRA	42.55	0.00745	0.97535	45.13	0.00554	0.98399	47.04	0.00445	0.98884
LEARN	39.43	0.01068	0.95866	42.52	0.00743	0.97716	44.29	0.00613	0.98340
PWLS-CSC	41.55	0.00837	0.97177	45.07	0.00558	0.98517	48.69	0.00368	0.99269
PWLS-CSCGR	42.78	0.00726	0.97687	46.07	0.00497	0.98753	49.54	0.00333	0.99372

Table 3. TABLE III: Quantitative results associated with the different number of filters.

Image		Abdominal				Thoracic
Num. of Filters	PSNR	RMSE	SSIM	Time	PSNR	RMSE	SSIM	Time
4	44.77	0.00577	0.97990	722 s	45.00	0.00562	0.98450	788 s
8	45.40	0.00537	0.98205	776 s	45.74	0.00517	0.98654	695 s
16	45.54	0.00528	0.98262	745 s	46.01	0.00500	0.98737	763 s
24	45.61	0.00524	0.98288	787 s	46.07	0.00497	0.98753	917 s
32	45.73	0.00517	0.98329	970 s	46.07	0.00497	0.98753	985 s
64	45.70	0.00519	0.98334	1069 s	46.18	0.00491	0.98779	1371 s
150	45.65	0.00522	0.98308	2087 s	46.03	0.00500	0.98752	2418 s

Table 4. TABLE IV: Quantitative results associated with the different numbers of training samples.

Image		Abdominal			Thoracic
Num. of Samples	PSNR	RMSE	SSIM	PSNR	RMSE	SSIM
1	45.43	0.00535	0.98236	45.48	0.00532	0.98595
3	45.69	0.00519	0.98313	45.94	0.00504	0.98717
5	45.68	0.00520	0.98317	45.94	0.00504	0.98727
10	45.73	0.00517	0.98329	46.07	0.00497	0.98753
50	45.68	0.00520	0.98313	46.12	0.00494	0.98764

Table 5. TABLE V: Quantitative results associated with the different training sets.

Image		Abdominal			Thoracic
Training Set	PSNR	RMSE	SSIM	PSNR	RMSE	SSIM
Undersampled CT Images	44.12	0.00622	0.97759	43.73	0.00651	0.98006
Natural Images	44.71	0.00581	0.97974	44.65	0.00586	0.98329
Intermediate Images	45.48	0.00532	0.98237	45.65	0.00522	0.98655
Full-Sampled CT Images	45.73	0.00517	0.98329	46.07	0.00497	0.98753

Equations

$${\sigma }_p^2=r_p\times exp({\bar{\mathbf{y}}}_p/\epsilon )$$

$${argmin}_{\mathbf{u}}(\mathbf{y}-\mathbf{Au}{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{y}-\mathbf{Au})+\beta R(\mathbf{u})$$

$${argmin}_{\{{\mathbf{M}}_i\},\{{\mathbf{f}}_i\}}\frac{1}{2}{\Vert \sum _{i=1}^N{\mathbf{f}}_i\ast {\mathbf{M}}_i-\mathbf{x}\Vert }_2^2+\lambda \sum _{i=1}^N{\Vert {\mathbf{M}}_i\Vert }_1,$$

$$\begin{array}{rl}& {argmin}_{\mathbf{u},\{{\mathbf{M}}_i\},\{{\mathbf{f}}_i\}}\frac{1}{2}(\mathbf{y}-\mathbf{Au}{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{y}-\mathbf{Au})+\\ & \beta \left(\frac{1}{2}{\Vert \sum _{i=1}^N{\mathbf{f}}_i\ast {\mathbf{M}}_i-\mathbf{u}\Vert }_2^2+\lambda \sum _{i=1}^N{\Vert {\mathbf{M}}_i\Vert }_1\right).\end{array}$$

$$\begin{array}{rl}& {argmin}_{\mathbf{u},\{{\mathbf{M}}_i\}}\frac{1}{2}(\mathbf{y}-\mathbf{Au}{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{y}-\mathbf{Au})+\\ & \beta \left(\frac{1}{2}{\Vert \sum _{i=1}^N{\mathbf{f}}_i\ast {\mathbf{M}}_i-\mathbf{u}\Vert }_2^2+\lambda \sum _{i=1}^N{\Vert {\mathbf{M}}_i\Vert }_1\right).\end{array}$$

$$\begin{array}{rl}{argmin}_{\mathbf{u}}& (\mathbf{y}-\mathbf{Au}{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{y}-\mathbf{Au})\\ & +\beta {\Vert \sum _{i=1}^N{\mathbf{f}}_i\ast {\tilde{\mathbf{M}}}_i-\mathbf{u}\Vert }_2^2.\end{array}$$

$$\begin{array}{rl}{\mathbf{u}}_q^{t+1}={\mathbf{u}}_q^t& -\frac{{\sum }_{p=1}^P\left(\left(1/{\sigma }_p^2\right)a_{pq}\left({\left[{\mathbf{Au}}^t\right]}_p-{\mathbf{y}}_p\right)\right)}{{\sum }_{p=1}^P\left(\left(1/{\sigma }_p^2\right)a_{pq}{\sum }_{k=1}^Q{a}_{pk}\right)+\beta }\\ & -\frac{\beta \left({\left[{\sum }_{i=1}^N{\mathbf{f}}_i\ast {\tilde{\mathbf{M}}}_i\right]}_q^t-{\mathbf{u}}_q^t\right)}{{\sum }_{p=1}^P\left(\left(1/{\sigma }_p^2\right)a_{pq}{\sum }_{k=1}^Q{a}_{pk}\right)+\beta },\\ & q=1,2,\dots ,Q,\end{array}$$

$${argmin}_{\{{\mathbf{M}}_i\}}\frac{1}{2}{\Vert \sum _{i=1}^N{\mathbf{f}}_i\ast {\mathbf{M}}_i-\mathbf{u}\Vert }_2^2+\lambda \sum _{i=1}^N{\Vert {\mathbf{M}}_i\Vert }_1.$$

$${argmin}_{\mathbf{M}}\frac{1}{2}{\Vert \mathbf{FM}-\mathbf{u}\Vert }_2^2+\lambda {\Vert \mathbf{M}\Vert }_1.$$

$${argmin}_{\mathbf{M},\mathbf{B}}\frac{1}{2}{\Vert \mathbf{FM}-\mathbf{u}\Vert }_2^2+\lambda {\Vert \mathbf{B}\Vert }_{1}s.t.\mathbf{M}=\mathbf{B}.$$

$${\mathbf{M}}^{j+1}={argmin}_{\mathbf{M}}\frac{1}{2}{\Vert \mathbf{FM}-\mathbf{u}\Vert }_2^2+\frac{\rho }{2}{\Vert \mathbf{M}-{\mathbf{B}}^j+{\mathbf{C}}^j\Vert }_2^2,$$

$${\mathbf{B}}^{j+1}={argmin}_{\mathbf{B}}\lambda {\Vert \mathbf{B}\Vert }_1+\frac{\rho }{2}{\Vert {\mathbf{M}}^{j+1}-\mathbf{B}+{\mathbf{C}}^j\Vert }_2^2,$$

$${\mathbf{C}}^{j+1}={\mathbf{C}}^j+{\mathbf{M}}^{j+1}-{\mathbf{B}}^{j+1}.$$

$${argmin}_{\hat{\mathbf{M}}}\frac{1}{2}{\Vert \hat{\mathbf{F}}\hat{\mathbf{M}}-\hat{\mathbf{u}}\Vert }_2^2+\frac{\rho }{2}{\Vert \hat{\mathbf{M}}-\hat{\mathbf{B}}+\hat{\mathbf{C}}\Vert }_2^2,$$

$$\left({\hat{\mathbf{F}}}^H\hat{\mathbf{F}}+\rho \mathbf{I}\right)\hat{\mathbf{M}}={\hat{\mathbf{F}}}^H\hat{\mathbf{u}}+\rho \left(\hat{\mathbf{B}}-\hat{\mathbf{C}}\right),$$

$${\mathbf{B}}^{j+1}=S_{\lambda /\rho }\left({\mathbf{M}}^{j+1}+{\mathbf{C}}^j\right),$$

$$S_l\left(\mathbf{x}\right)=sign\left(\mathbf{x}\right)max\left(0,\left|\mathbf{x}\right|-l\right).$$

$$\begin{array}{rl}& {argmin}_{\mathbf{u},\{{\mathbf{M}}_i\},\{{\mathbf{f}}_i\}}\frac{1}{2}(\mathbf{y}-\mathbf{Au}{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{y}-\mathbf{Au})+\\ & \beta (\frac{1}{2}{\Vert \sum _{i=1}^N{\mathbf{f}}_i\ast {\mathbf{M}}_i-\mathbf{u}\Vert }_2^2+\lambda \sum _{i=1}^N{\Vert {\mathbf{M}}_i\Vert }_1+\\ & \frac{\tau }{2}\sum _{i=1}^N{\Vert \sqrt{({\mathbf{g}}_0\ast {\mathbf{M}}_i{)}^2+({\mathbf{g}}_1\ast {\mathbf{M}}_i{)}^2}\Vert }_2^2).\end{array}$$

$$\begin{array}{rl}{argmin}_{\{{\mathbf{M}}_i\}}& \frac{1}{2}{\Vert \sum _{i=1}^N{\mathbf{f}}_i\ast {\mathbf{M}}_i-\mathbf{u}\Vert }_2^2+\lambda \sum _{i=1}^N{\Vert {\mathbf{M}}_i\Vert }_1+\\ & \frac{\tau }{2}\sum _{i=1}^N{\Vert \sqrt{({\mathbf{g}}_0\ast {\mathbf{M}}_i{)}^2+({\mathbf{g}}_1\ast {\mathbf{M}}_i{)}^2}\Vert }_2^2,\end{array}$$

$$\frac{\tau }{2}\sum _{i=1}^N\left({\Vert {\mathbf{G}}_0{\mathbf{M}}_i\Vert }_2^2+{\Vert {\mathbf{G}}_1{\mathbf{M}}_i\Vert }_2^2\right).$$

$$\begin{array}{rl}{argmin}_{\mathbf{M}}& \frac{1}{2}{\Vert \mathbf{FM}-\mathbf{u}\Vert }_2^2+\lambda {\Vert \mathbf{M}\Vert }_1+\\ & \frac{\tau }{2}{\Vert {\mathbf{\Phi }}_0\mathbf{M}\Vert }_2^2+\frac{\tau }{2}{\Vert {\mathbf{\Phi }}_1\mathbf{M}\Vert }_2^2,\end{array}$$

$${\mathbf{\Phi }}_l=\left(\begin{array}{ccc}{\mathbf{G}}_l& 0& \cdots \\ 0& {\mathbf{G}}_l& \cdots \\ ⋮& ⋮& \ddots \end{array}\right).$$

$$\begin{array}{rl}{argmin}_{\mathbf{M},\mathbf{B}}& \frac{1}{2}{\Vert \mathbf{FM}-\mathbf{u}\Vert }_2^2+\frac{\tau }{2}{\Vert {\mathbf{\Phi }}_0\mathbf{M}\Vert }_2^2+\\ & \frac{\tau }{2}{\Vert {\mathbf{\Phi }}_1\mathbf{M}\Vert }_2^2+\lambda {\Vert \mathbf{B}\Vert }_1,s.t.\mathbf{M}=\mathbf{B}.\end{array}$$

$$\begin{array}{rl}{\mathbf{M}}^{j+1}=& {argmin}_{\mathbf{M}}\frac{1}{2}{\Vert \mathbf{FM}-\mathbf{u}\Vert }_2^2+\frac{\tau }{2}{\Vert {\mathbf{\Phi }}_0\mathbf{M}\Vert }_2^2\\ & +\frac{\tau }{2}{\Vert {\mathbf{\Phi }}_1\mathbf{M}\Vert }_2^2+\frac{\rho }{2}{\Vert \mathbf{M}-{\mathbf{B}}^j+{\mathbf{C}}^j\Vert }_2^2,\end{array}$$

$${\mathbf{B}}^{j+1}={argmin}_{\mathbf{B}}\lambda {\Vert \mathbf{B}\Vert }_1+\frac{\rho }{2}{\Vert {\mathbf{M}}^{j+1}-\mathbf{B}+{\mathbf{C}}^j\Vert }_2^2,$$

$${\mathbf{C}}^{j+1}={\mathbf{C}}^j+{\mathbf{M}}^{j+1}-{\mathbf{B}}^{j+1}.$$

$$\begin{array}{rl}& \left({\hat{\mathbf{F}}}^H\hat{\mathbf{F}}+\tau {\hat{\mathbf{\Phi }}}_0^H{\hat{\mathbf{\Phi }}}_0+\tau {\hat{\mathbf{\Phi }}}_1^H{\hat{\mathbf{\Phi }}}_1+\rho \mathbf{I}\right)\hat{\mathbf{M}}\\ & ={\hat{\mathbf{F}}}^H\hat{\mathbf{u}}+\rho \left(\hat{\mathbf{B}}-\hat{\mathbf{C}}\right),\end{array}$$

Full text

I INTRODUCTION

Over the past few decades, computed tomography (CT) has been very successfully used for clinical diagnosis and intervention. However, X-rays can be absorbed partially by the human body and may cause genetic or cancerous diseases [1]. As a result, the famous ALARA (as low as reasonably achievable) principle has been widely utilized to avoid an excessive radiation dose in clinical practice. Reducing the X-ray flux towards each detector element and decreasing the number of sampling views are two effective ways to reduce the radiation dose, resulting in noisy measurements or insufficient projection data, respectively. However, it is chanllenging to generate high-quality images with degraded projection data using traditional methods, such as filtered back projection (FBP). Thus, reconstructing high-quality images from contaminated and/or undersampled projection data has attracted increased attention in the field of CT imaging. In this work, we focus on the sparse-view problem.

To recover from the situation of undersampling, iterative reconstruction methods have been developed, which can efficiently improve the image quality when the projection views are incomplete. Classical iterative reconstruction algorithms, such as the algebraic reconstruction technique (ART) [2], simultaneous the algebraic reconstruction technique (SART) [3] and expectation maximization (EM) [4], can mitigate this problem to a certain degree. However, it is difficult to achieve a satisfactory result when projection views are highly sparse without extra priors. Introducing reasonable priors can significantly improve the imaging quality of traditional iterative methods. In recent years, compressed sensing (CS) theory has been developed and has been proven to be a powerful tool, which indicates that a signal can be exactly reconstructed with a very high probability if the signal can be represented sparsely with a certain sparsifying transform [5, 6]. Inspired by CS theory, Sidky et al. [7, 8] coupled total variation (TV) and projection onto convex sets (POCS) to solve incomplete projection data reconstruction problems and achieved good performance. However, TV assumes that the signal is piecewise smooth, resulting in undesired patchy effects on the reconstructed images. Following this path, many variants of TV have been proposed, such as adaptive-weighted TV [9], fractional-order TV [10, 11], and nonlocal TV [12, 13]. Furthermore, Niu et al. [14] combined penalized weighted least-squares (PWLS) [15] with total generalized variation (TGV) [16] to reduce the patchy effects. Chen and Lang indicated that a high-quality image can be used as the prior to accurately reconstruct CT dynamic image sequences from highly undersampled projection data [17].

Recently, inspired by the rapid development of deep learning, several pioneering studies combining the ideas of deep learning and CS were presented [18, 19]. Typically, Chen et al. proposed to unfold the iterative reconstruction into a residual convolutional neural network, and the parameters and the regularization terms can be learned from the external dataset [20]. Alder and Öktem employed the primal dual hybrid gradient (PDHG) algorithm for non-smooth convex optimization and unrolled the optimization process to a neural network [21]. Zhu et al. proposed a unified image reconstruction framework–automated transform by manifold approximation (AUTOMAP)–which directly learned a potential mapping function from the sensor domain to the image domain [22]. In [23], a parameterized plug-and-play alternating direction method of multipliers (3pADMM) was developed for the PWLS model, and then the prior terms and related parameters were optimized with a neural network with a large quantity of training data. Except for the reconstruction problem, deep learning was also used in low-dose CT (LDCT) denoising. As a pioneering work, Chen et al. used a three-layer convolutional neural network (CNN) to handle the noise in the LDCT images [24]. Kang et al. introduced contourlet-based multi-scale U-Net to improve the image quality of LDCT [25]. Motivated by the idea of autoencoders, Chen et al. proposed a residual encoder-decoder CNN (RED-CNN) for LDCT [26]. These methods obtained promising results, but the requirement of a large number of training samples may limit the application scenarios.

As another direction of CS, in recent years, dictionary learning (DL)-based methods have been proven useful in various image restoration problems. In 2006, Elad and Aharon adopted a dictionary-based method to address the image denoising problem [27]. Mairal et al. explored how to apply this method to color image restoration [28]. In the field of medical imaging, the DL method was first introduced into magnetic resonance imaging (MRI) by Chen et al. [29]. Ravishankar and Bresler utilized an adaptive dictionary-based method to reconstruct MR images from the highly undersampled $k$-space data [30]. Later, Xu et al. first introduced DL into LDCT reconstruction [31]. Two types of dictionaries were learned as bases: (1) a global dictionary was learned from an external training set that was fixed during the reconstruction procedure, and (2) an adaptive dictionary was learned from the intermediate result that continued to update during the iterations. Inspired by the work of super-resolution, which uses two dictionaries to link low-resolution images and the corresponding high-resolution training images [32], Lu et al. proposed using two dictionaries to predict full-sampled CT images [33], and Zhao et al. extended this method to spectral CT [34]. Chen et al. further improved the DL-based artifact reduction method with tissue feature atoms and artifact atoms to build discriminative dictionaries [35]. Liu et al. extended this idea to 3D sinogram restoration for LDCT [36]. Zheng et al. proposed to learn a union of sparsifying transforms instead of dictionary atoms for LDCT image reconstruction [37].

However, most dictionary learning methods are patch-based, and the learned features often contain shifted versions of the same features. The learned dictionary may be over-redundant, and the results will suffer from patch aggregation artifacts on the patch boundaries. To address these problems, convolutional sparse coding (CSC) methods have been proposed in which the shift invariance is directly modeled in the objective function. CSC can be traced back to [38] and was first proposed by Zeiler et al. [39]. The effectiveness of CSC has been demonstrated in several computer vision tasks. Gu et al. introduced CSC into super-resolution and achieved promising results [40]. Wohlberg approached impulse noise via CSC aided by gradient constraints [41]. Liu et al. applied CSC for image fusion [42]. More recently, Zhang and Patel combined CSC and low-rank decomposition to address rain streak removal and image decomposition problems [43, 44]. Thanh *et al.*proposed to adopt a shared 3D multi-scale convolutional dictionary to recover the high-frequency information of the MRI images [45].

Inspired by the pioneering studies of CSC, in this work, two models based on CSC are proposed to address the sparse-view CT reconstruction problem. In the first model, which is called PWLS-CSC, a CSC regularization is adapted with the PWLS model. The filters are trained with an external full-sampled CT dataset. Although PWLS-CSC can work well for sparse-view CT, ringing artifacts may appear in some reconstructed images due to inaccurate filters. On the other hand, gradient constraint is an effective way to suppress these types of artifacts. Therefore, in the second model, an improved version of PWLS-CSC, which imposes gradient regularization on feature maps, is proposed. For simplicity, we called this method PWLS-CSCGR.

The remainder of this paper is organized as follows. In Section II, we introduce the background on PWLS and CSC. In Section III, the proposed reconstruction schemes are elaborated. In Section IV, the experimental designs and representative results are given. Finally, we will discuss relevant issues and conclude this paper in Section V.

II BACKGROUND

II-A Penalized Weighted Least-Squares for CT Image Reconstruction

The calibrated and log-transformed projection data approximately follow a Gaussian distribution, which can be described by the following analytic formula [46]:

[TABLE]

where $\bar{\bm{y}}_{p}$ and $\sigma_{p}^{2}$ are the mean and variance of the measured projection data at $p$-th bin, respectively, $r_{p}$ is a parameter adaptive to different bins and $\varepsilon$ is a scaling parameter. Based on these properties of the projection data, the PWLS-based CT reconstruction problem can be modeled as follows:

[TABLE]

where $\bm{y}$ denotes the projection data, $\bm{u}$ is the vectorized attenuation coefficients to be reconstructed and $(\cdot)^{T}$ denotes the transpose operation. $\bm{A}$ denotes the system matrix with a size of $P\times Q$ ($P$ is the total number of projection data and $Q$ is the total number of image pixels). $\bm{\Sigma}$ is a diagonal matrix with the $p$-th element of $\sigma_{p}^{2}$ as calculated by (1). $R(\bm{u})$ represents the regularization term, and $\beta$ is a balancing parameter to control the tradeoff between the data-fidelity and regularization terms.

II-B Convolutional Sparse Coding

In contrast to traditional patch-based DL methods, CSC assumes that the image can be represented as a summation of convolutions between a set of filters and its corresponding feature maps:

[TABLE]

where $\bm{x}$ denotes the image, $\ast$ denotes the convolution operator, $\{\bm{f}_{i}\}_{i=1,2,\dots,N}$ is a set of filters, $\{\bm{M}_{i}\}_{i=1,2,\dots,N}$ are the feature maps corresponding to filters $\bm{f}_{i}$, and $\lambda$ is the regularization parameter. In the CSC model, the feature map $\bm{M}_{i}$ has the same size as image $\bm{x}$. Since the reconstructed image is obtained by the summation of convolution operations $\bm{f}_{i}\ast\bm{M}_{i}$, the patch aggregation artifacts in traditional patch-based sparse coding methods can be avoided. Meanwhile, the features learned by patch-based sparse coding methods often contain shifted versions of the same features, and the filters learned by CSC can obtain rotation-invariant features, which makes this method more efficient.

III METHOD

III-A PWLS-CSC Model

Current TV- or DL-based methods suffer from blocky effects or patch aggregation artifacts. To circumvent these problems, we propose introducing CSC as the regularization term in (2). Combining CSC with (2), we aim to solve the following optimization problem:

[TABLE]

In this model, we use an external full-sampled CT dataset to train the filters $\{\bm{f}_{i}\}$ by the method proposed in [47], so (4) becomes the following optimization problem:

[TABLE]

An alternating minimization scheme can be used to solve (5). First, an intermediate reconstructed image $\bm{u}$ can be obtained with a set of fixed feature maps $\{\widetilde{\bm{M}}_{i}\}$. Thus, the objective function (5) becomes:

[TABLE]

With the separable paraboloid surrogate method [48], the solution of (6) can be obtained by:

[TABLE]

where $t=0,1,2,\dots,T$ represents the iteration index and $a_{pq}$ is the element of system matrix $\bm{A}$. The second step is to represent the intermediate result $\bm{u}$ obtained in the last step with the fixed filters $\{\bm{f}_{i}\}$, which means calculating $\{\bm{M}_{i}\}$. Since CSC does not provide a good representation of the low-frequency component of the image [41], only the high-frequency component of $\bm{u}$ is processed by CSC, and the results can be obtained by simply summing the recovered high-frequency component and the low-frequency component. The high-frequency component of the image is precomputed by a high-pass filtering operation. Note that the symbol $\bm{u}$ in the remainder of this section denotes the high frequency component. Then, we have the following optimization problem:

[TABLE]

Let $\bm{F}_{i}\bm{M}_{i}=\bm{f}_{i}\ast\bm{M}_{i}$, $\bm{F}=\left(\bm{F}_{1}\ \bm{F}_{2}\ \cdots\ \bm{F}_{N}\right)$ and $\bm{M}=\left(\bm{M}_{1}\ \bm{M}_{2}\ \cdots\ \bm{M}_{N}\right)^{T}$, (8) can be rewritten as:

[TABLE]

The alternating Direction Method of Multipliers (ADMM) [49] is adopted to solve (9) by introducing a dual variable $\bm{B}$ that is constrained to be equal to $\bm{M}$, leading to the following problem:

[TABLE]

$\bm{M}$ and $\bm{B}$ are optimized alternately as follows:

[TABLE]

[TABLE]

[TABLE]

(11) can be solved by transforming it into Fourier domain [50, 51], which gives us:

[TABLE]

where $\hat{\bm{F}}$, $\hat{\bm{M}}$, $\hat{\bm{u}}$, $\hat{\bm{B}}$ and $\hat{\bm{C}}$ denote the corresponding expressions in the Fourier domain. (14) can be solved with the following closed-form solution:

[TABLE]

where $\bm{I}$ is the identity matrix. (15) can be computed efficiently with the Sherman Morrison formula [52]. The closed-form solution of (12) can be obtained by:

[TABLE]

where $S(\cdot)$ denotes the soft-thresholding function [53], which can be calculated as follows:

[TABLE]

The overall PWLS-CSC algorithm for sparse-view CT reconstruction is summarized in Algorithm 1.

III-B PWLS-CSCGR

As claimed in [31], for DL-based methods, the structure loss or new artifacts may appear due to the inaccurate dictionary atoms. Similarly, CSC-based methods will suffer this problem if inaccurate filters exist. The constraint with gradient regularization (isotropic total variation), which is usually used to suppress outliers, is an effective way to overcome this problem [41, 54]. Furthermore, as suggested by [54], imposing the gradient constraint on the feature maps is superior to applying the constraint on the image domain. As a result, to further improve the proposed PWLS-CSC, in this subsection, we propose to utilize the gradient regularization (PWLS-CSCGR) on feature maps as an additional constraint in PWLS-CSC, which yields the following optimization problem:

[TABLE]

where $\bm{g}_{0}$ and $\bm{g}_{1}$ are the filters that compute the gradient along the x- and y-axes respectively, and $\tau$ is a hyperparameter. Similar to (5), given the predetermined filters $\{\bm{f}_{i}\}$, (18) can be decomposed into two optimization problems. One problem is exactly the same as (6), and the other can be given as:

[TABLE]

The linear operators $\bm{G}_{0}$ and $\bm{G}_{1}$ are defined such that $\bm{G}_{l}\bm{M}_{i}=\bm{g}_{l}\ast\bm{M}_{i}$ and the gradient term in (19) can be rewritten as:

[TABLE]

By introducing block matrix notation, (19) can be rewritten as follows:

[TABLE]

where

[TABLE]

To apply ADMM to (21), we introduce a dual variable $\bm{B}$ and have the following optimization problem:

[TABLE]

Then we have an iterative scheme for (23) as follows:

[TABLE]

[TABLE]

[TABLE]

The only difference between solving (23) and (10) is solving $\bm{M}$. The solution for (24) in the Fourier domain is given by:

[TABLE]

that can also be solved by the Sherman Morrison formula [47]. The main steps of the proposed PWLS-CSCGR algorithm are summarized in Algorithm 2.

IV EXPERIMENTAL RESULTS

In this section, extensive experiments were conducted to demonstrate the performance of the proposed methods. Two types of data, including simulated and real clinical data, were used in our experiments. All experiments were performed in MATLAB 2017a on a PC (equipped with an AMD Ryzen 5 1600 CPU at 3.2 GHz and 16 GB RAM). To accelerate the algorithm, we performed the CSC and CSCGR phases with a graphics processing unit (GPU) GTX 1080 Ti. All filters and the regularization parameters used in this section were the same if there was no special instruction, including 32 filters with sizes of 10 $\times$ 10, $\beta$ = 0.005, $\lambda$ = 0.005, $\rho$ = 100 $\times\ \lambda$ + 1 and $\tau$ = 0.06. Notably, in our numerical scheme, the computation of convolution was transformed into the frequency domain, which can be treated as the dot product operation, so there was no need to keep the filter size to an odd number. All the parameters were constants except $\rho$, which was self-adaptive according to [55]. The peak signal-to-noise ratio (PSNR), root-mean-square error (RMSE) and structural similarity (SSIM) [56] were utilized to quantitatively evaluate the performance of the methods.

IV-A Simulated Results

To demonstrate the performance of the proposed methods, a dataset of clinical CT images, authorized by Mayo Clinic and downloaded from ”the NIH-AAPM-Mayo Clinic Low Dose CT Grand Challenge,” was used [57]. This dataset contains 5,936 1 mm thick full dose CT slices collected from 10 patients. The high-frequency components of the ten images randomly selected from this dataset were used to train the predetermined filters $\{\bm{f}_{i}\}$. Fig. 1 shows the images in the training set. The projection data in fan-beam geometry were simulated by the Siddon’s ray-driven algorithm [58]. A total of 64 projection views were evenly distributed over $360^{o}$, and a flat detector with a length of 41.3 cm had 512 bins. The image arrays were 20 $\times$ 20 cm2. The distance from the source to the center of rotation was 40 cm, and the distance from the center of the detector to the center of rotation was 40 cm. The outer iteration number was set to 2000. The inner iteration numbers were set to 20 and 100 for PWLS and CSCGR, respectively. The iteration numbers for PWLS-CSC were the same as those for PWLS-CSCGR.

Several state-of-the-art methods were compared with our proposed algorithms, including FBP, PWLS-TGV [14], PWLS-DL [31], PWLS-ULTRA[37] and LEARN [20].

IV-A1 Abdominal Case

The original abdominal image and the results reconstructed from different methods are shown in Fig. 2. The result of FBP contains severe streak artifacts and is clinically almost useless, while all other algorithms efficiently remove the artifacts. However, in Fig. 2(c), blocky effects are noticed to a certain degree. PWLS-DL effectively avoids the blocky effects in Fig. 2(d), but the result appears smooth in some organs. This effect is because the DL method is patch-based, which usually averages all the overlapped patches to produce the final result. This procedure may efficiently suppress the noise, but some small details may also be smoothened. PWLS-ULTRA suffers the same problem as PWLS-DL. In Fig. 2(f), the recent deep learning-based method LEARN also fails to recover small details. In Fig. 2(g) and (h), the proposed methods maintain more details while the artifacts are removed, especially some contrast-enhanced blood vessels in the liver. The differences between PWLS-CSC and PWLS-CSCGR are almost invisible to visual inspection, but the absolute difference images, which are illustrated in Fig. 3, show that the PWLS-CSCGR exhibits better performance.

To better visualize the details, Fig. 4 shows the results of a magnified region of interest (ROI), which is indicated by the red box in Fig. 2(a). The red arrows indicate several contrast-enhanced blood vessels in the liver, which can only be well identified by our methods. In particular, in Fig. 4(d)-(f), the liver region is smoothened, and the contrast-enhanced blood vessels are almost lost. In Fig. 4(g) and (h), the proposed methods preserve most details and have the most coherent visual effect on the reference image, even for the mottle-like structures in the liver.

IV-A2 Thoracic Case

Fig. 5 presents the thoracic images reconstructed by different methods. It can be easily observed that the FBP result is covered by the artifacts, and it is difficult to clearly discriminate the artifacts and blood vessels in the lung. In Fig. 5(c)-(f), most artifacts are suppressed, but the results are smoothened to different degrees in some organs, and some important details are lost. In Fig. 5(g) and (h), more details in some regions are preserved while the artifacts are effectively removed.

In particular, two small regions, which are indicated by the blue and red boxes in Fig. 5(a), are enlarged in Fig. 6 and Fig. 7. In Fig. 6(b), the heart could not be well recognized. PWLS-TGV distorted the boundaries and the structures in the heart in Fig. 6(c). As indicated by the red arrow, in Figs. 6(d)-(f), the tiny vessels is blurred, and only PWLS-CSC and PWLS-CSCGR satisfactorily recover this structural information. In Fig. 6(g), the result of PWLS-CSC introduces ringing artifacts, which are probably caused by inaccurate filters, while PWLS-CSCGR effectively suppresses these artifacts with the aid of gradient regularization. In Fig. 7, the curve-like structures indicated by the red arrow are clearly visible in the results of PWLS-CSC and PWLS-CSCGR, and all the other methods cannot adequately maintain these structures.

IV-A3 Quantitative Evaluation

To demonstrate the robustness of the proposed methods at various sampling rates, experiments with 48, 64 and 80 projection views were performed. Table I tabulates the quantitative results of the abdominal image, which is shown in Fig. 2(a) and reconstructed from 48, 60 and 80 projection views with different methods. When the number of projection views is set to 48, PWLS-DL and PWLS-ULTRA demonstrate better performance than PWLS-CSC, but PWLS-CSC outperforms them when the number of projection views is set to 64 or 80 and PWLS-CSCGR achieves the best performance in all situations. Table II displays the quantitative evaluations for the thoracic image shown in Fig. 5(a) and reconstructed from 48, 60 and 80 projection views with different methods. A similar trend is shown in Table II, in which PWLS-CSCGR obtains the highest scores in all three metrics and PWLS-CSC also outperforms PWLS-TGV and PWLS-DL for these conditions. The plausible reason is that the thoracic image has more small details, which can be better represented by CSC. In summary, the proposed PWLS-CSC and PWLS-CSCGR methods are more robust to different sampling rates than other methods.

IV-B Real Data

To evaluate the proposed methods in a real clinical configuration, in this subsection, projection data acquired on a UEG Dental CT scanner using a normal-dose protocol of 110 kVp and 13 mA/10 ms per projection were used. Five hundred sampling angles were evenly distributed over 360 degrees. In each projection view, 950 bins were distributed defining a field-of-view (FOV) of 18 $\times$ 18 $\times$ 25 cm3. The reconstructed image is a matrix of 256 $\times$ 256 pixels. The outer iteration number was set to 40, and the inner iteration numbers for PWLS and CSCGR were set to 2 and 100, respectively. The iteration numbers for PWLS-CSC are the same as those for PWLS-CSCGR.

Fig. 8 shows the reference image reconstructed by ART with full-sampled (500 views) projection data and the results reconstructed by various methods with downsampled projection data with 125 views. It can be observed that due to the limitation of hardware, noticeable noise and artifacts still exist in the full-sampled result in Fig. 8(a). In Fig 8(b), FBP suffers from severe artifacts that fulfill every area of the mouth. In Fig. 8(c), PWLS-TGV completely removes most of the artifacts, but the blocky effects at the edges of tissue, which are indicated by the blue and yellow arrows, are obvious. PWLS-DL efficiently suppresses the artifacts shown in the reference image. However, some artifacts caused by patch aggregation can be observed in Fig. 8(d), which are indicated by the red arrow, and the details of the teeth are blurred. In Fig. 8(e) and (f), PWLS-ULTRA and LEARN cannot remove all the artifacts, and the details of teeth are blurred. In Fig. 8(g) and (h), PWLS-CSC and PWLS-CSCGR efficiently eliminate the noise and artifacts and avoid the side effect caused by patch aggregation in PWLS-DL, and PWLS-CSCGR clearly preserves the details of teeth.

IV-C Effect of Filter Parameters

In this subsection, we evaluate the impact of filter parameters on PWLS-CSCGR, including the number of filters, the filter sizes, the number of training samples and the different types of training sets.

IV-C1 Number of Filters

To evaluate the impact of the number of filters, experiments were performed with $N$ = 4, 8, 16, 24, 32, 64 and 150 on the selected images in subsection IV-A. The quantitative results are summarized in Table III. It can be seen that the metrics for both images rise at first and then begin to stabilize when $N$ continues to increase. This phenomenon may rely on the fact that when $N$ is small, the performance can be improved significantly by adding more filters due to the increasing precision of representation, but when $N$ is large enough, the benefit will also decay rapidly. Based on this observation, to balance the performance and computational time, the number of filters was set to 32 in this paper.

IV-C2 Filter Sizes

To sense the impact of filter sizes, different sizes including 6, 8, 10, 12 and 14 were tested. The results are given in Fig. 9. The results show that the performance of our method is not quite sensitive to the filter size and that there are no distinct differences for different sizes. As a result, we simply set filter size to 10, which is approximately the maximum in Fig. 9.

IV-C3 Number of Training Samples

In this subsection, different numbers of training samples are used to validate the impact of the size of the training set. All images used were randomly selected from the Mayo dataset. The quantitative results are given in Table IV. Notably, as the number of training samples increases, the improvement in the performance is small, which can be treated as the evidence that for the current number of filters, the performance is robust to the size of the training set.

IV-C4 Filter Trained with Different Training Sets

In this subsection, we consider different types of images to form the training set. Three datasets, including undersampled CT images, natural images and the intermediate images in the reconstruction procedure, were tested. Fig. 10 shows the undersampled CT images reconstructed by FBP. Fig. 11 shows the natural image dataset. The undersampled CT images and natural images were used to train predetermined filters. The adaptive filters are trained with the intermediate images in the reconstruction procedure and updated in each iteration, which is similar to the adaptive dictionary learning in [31]. Thirty-two filters with a size of 10 $\times$ 10 are trained with different training sets, and the filters are shown in Fig. 12 (the adaptive filters are trained with the intermediate result of the abdominal image). It can be observed that the learned adaptive filters become more sharper and cleaner in Fig. 12(c)-(e), and the two pretrained filters set in Fig. 12(a) and (b) also capture the main structures although the filters look noisy. Moreover, in Fig. 13, the reconstruction results with different filters trained with various image sets appear to be visually similar, and it can be seen that the quantitative results shown in Table V are also close. Notably, although the results reconstructed by the adaptive filters achieve high scores, it is time consuming to update the filters in each iteration. In summary, our algorithm is robust to the choice of filters (either predetermined or adaptively learned).

IV-D Convergence analysis and the Effect of Regularization Parameters

IV-D1 Convergence analysis

Fig. 14 shows the PSNR and RMSE versus the iteration steps for the PWLS-CSCGR. It can be seen that with an increase in the iteration number, both PSNR and RMSE curves converge to a stable position. These observations demonstrate that the proposed numerical scheme can efficiently optimize the energy functions to a satisfactory solution.

IV-D2 Effect of $\beta$

To investigate the sensitivity of $\beta$, experiments were performed for various $\beta$ values. The results are presented in Fig. 15. The PSNR and SSIM decrease with increasing $\beta$. However, according to our experiments, the smaller $\beta$ is, the slower the algorithm converges. When $\beta$ is less than or equal to 0.0025, the computational cost is twice or more when $\beta$ is equal to 0.005. Therefore, we set $\beta$ to be 0.005 in our algorithms.

IV-D3 Effect of $\lambda$

To investigate the sensitivity of our method to $\lambda$, experiments were performed for various $\lambda$ values. The results are shown in Fig. 16. It can be observed that PSNR and SSIM for both images demonstrate similar trends. The values of PSNR and SSIM rise rapidly with an increase in $\lambda$ and reach a peak when $\lambda$ is equal to 0.005. PSNR and SSIM begin to slowly decrease when $\lambda$ is larger than 0.005. Thus, $\lambda$ was set to 0.005 in this paper.

IV-D4 Effect of $\tau$

To explore the impact of $\tau$, experiments were performed for various $\tau$ values. The results are given in Fig. 17. It is obvious that PSNR and SSIM change slowly with $\tau$, which means that PWLS-CSCGR is not sensitive to $\tau$. According to the experimental results, $\tau$ was set to 0.06 in our experiments.

V DISCUSSION AND CONCLUSION

With the development of CSC in recent years, CSC has been proven useful in many imaging problems, including super-resolution, image fusion, image decomposition and so on. Instead of dividing an image into overlapped patches, CSC directly works on the whole image, which maintains more details and avoids artifacts caused by patch aggregation. In this paper, we propose two methods based on CSC. The basic version introduces CSC into the PWLS reconstruction framework. To further improve the performance and preserve more structural information, gradient regularization on feature maps is imposed into the basic version. Qualitative and quantitative results demonstrate the merits of our methods.

In the experiments of Section IV-A and IV-B, the filters and parameters were the same, showing the generalization of the proposed methods and that there is no need to adjust the filters or the parameters patient by patient. We also examined the impacts of filters on our method. The experimental results show that PWLS-CSCGR can work well even with only four filters. PWLS-CSCGR is also robust to the training set or even without the training set and can be treated as an unsupervised learning method.

Importantly, another issue is the computational time. The main cost of our methods depends on two parts: training the filters and the reconstruction. Training 32 filters with 10 images costs 85 s of GPU. Because this operation is offline and there is no need for a large training set, this part will not be the main problem. On the other hand, the reconstruction is time-consuming. Although our methods have a similar heavy computational burden to PWLS-DL, several techniques, including parallel computing and advanced optimization methods, can be applied for acceleration.

One of the most important deep learning models is CNN, which is also based on the convolution operator. For CSC, a signal can be represented by a summation of convolutions between a set of filters and the corresponding feature maps, and the key point is to calculate the feature maps with certain (predetermined or adaptive) filters. CNN trains the cascaded filters to convolve with the inputs. Furthermore, current CNN-based methods still lack theoretical proof. Most deep learning methods are data-driven, and the results cannot be guaranteed without sufficient training data. However, CSC, as an unsupervised learning method, has a strict mathematical proof. This method is robust to the number of training samples (as shown in Sec. IV-C.3 and IV-C.4) and even without training data. On the other hand, the same groups analyzed the relationship between the CSC and CNN methods in [59, 60] and found that assuming that our signals originate from the multi-layer CSC model, the layered-thresholding pursuit algorithm for decomposing a given measurement vector $Y$ completely equals the forward propagation in CNNs. This interesting finding provides a new way to explore the interpretability of deep learning.

In conclusion, inspired by successful applications of CSC in the field of signal processing, we explored the potential of this method incorporating a PWLS image reconstruction framework, resulting in two novel algorithms referred to as PWLS-CSC and PWLS-CSCGR. We evaluated the proposed algorithms with simulated and real data. In the experimental results, our methods have been shown to be competitive with several state-of-art methods. The robustness of our methods was also investigated by extensive analysis with experimental configurations. In our future work, we will extend our methods to other CT imaging topics, such as metal artifact reduction and LDCT. Furthermore, the combination with deep learning-based methods is also an interesting direction.