Two-layer Residual Sparsifying Transform Learning for Image Reconstruction

TL;DR

This paper introduces a two-layer residual sparsifying transform model for improved image reconstruction, demonstrating enhanced performance in low-dose CT imaging through efficient optimization algorithms.

Contribution

It proposes a novel two-layer transform learning framework that better captures residual sparsity, advancing image reconstruction methods with efficient algorithms.

Findings

Two-layer model outperforms previous schemes in CT reconstruction.

Efficient block coordinate descent algorithms enable practical implementation.

Preliminary results show improved image quality from limited measurements.

Abstract

Signal models based on sparsity, low-rank and other properties have been exploited for image reconstruction from limited and corrupted data in medical imaging and other computational imaging applications. In particular, sparsifying transform models have shown promise in various applications, and offer numerous advantages such as efficiencies in sparse coding and learning. This work investigates pre-learning a two-layer extension of the transform model for image reconstruction, wherein the transform domain or filtering residuals of the image are further sparsified in the second layer. The proposed block coordinate descent optimization algorithms involve highly efficient updates. Preliminary numerical experiments demonstrate the usefulness of a two-layer model over the previous related schemes for CT image reconstruction from low-dose measurements.

Tables

Table 1. Table 1 : RMSE in HU (first row) and PSNR in dB (second row) of fan-beam reconstructions with FBP, PWLS-EP, PWLS-ST, and PWLS-MRST2 for two slices and three incident photon intensities.

		Slice $1$			Slice $2$
Intensity	10000	5000	3000	10000	5000	3000
FBP	73.7	89.0	101.0	72.5	86.1	112.2
	27.3	25.7	23.5	27.5	26.0	23.7
EP	39.4	49.7	56.9	37.1	45.5	53.5
	32.8	30.8	29.6	33.3	31.5	30.1
ST	36.5	43.9	49.4	33.7	41.5	49.0
	33.4	31.9	30.8	34.1	32.3	30.9
MRST2	$35.7$	$42.7$	$48.9$	$33.0$	$40.8$	$48.6$
	$33.6$	$32.0$	$30.9$	$34.3$	$32.5$	$31.0$

Equations

$$\displaystyle\min_{\{\mathbf{\Omega}_{l},\mathbf{Z}_{l}\}}\sum_{l=1}^{2}\bigg{\{}\|\mathbf{\Omega}_{l}\mathbf{R}_{l}-\mathbf{Z}_{l}\|^{2}_{F}+\eta_{l}^{2}\|\mathbf{Z}_{l}\|_{0}\bigg{\}}$$

$$s.t.{\mathbf{R}}_2={\mathbf{\Omega }}_1{\mathbf{R}}_1-{\mathbf{Z}}_1,{\mathbf{\Omega }}_l^T{\mathbf{\Omega }}_l=\mathbf{I}\forall l,$$

$$\underset{\mathbf{x}\ge 0}{\min }\frac{1}{2}\Vert \mathbf{y}-\mathbf{Ax}{\Vert }_{\mathbf{W}}^2+\beta \mathsf{S}(\mathbf{x}),\text{vspace}-0.05in$$

$$\displaystyle\mathsf{S}(\mathbf{x})\triangleq\min_{\{\mathbf{Z}_{l}\}}\sum_{l=1}^{2}\bigg{\{}\|\mathbf{\Omega}_{l}\mathbf{R}_{l}-\mathbf{Z}_{l}\|^{2}_{F}+\gamma_{l}^{2}\|\mathbf{Z}_{l}\|_{0}\bigg{\}}$$

$$s.t.{\mathbf{R}}_2={\mathbf{\Omega }}_1{\mathbf{R}}_1-{\mathbf{Z}}_1,{\mathbf{r}}_1^j={\mathbf{P}}^j\mathbf{x}\forall j.$$

$$\underset{{\mathbf{Z}}_1}{\min }\Vert {\mathbf{\Omega }}_1{\mathbf{R}}_1-{\mathbf{Z}}_1{\Vert }_F^2+\Vert {\mathbf{\Omega }}_2{\mathbf{R}}_2-{\mathbf{Z}}_2{\Vert }_F^2+{\eta }_1^2\Vert {\mathbf{Z}}_1{\Vert }_0.\text{vspace}-0.05in$$

$$\underset{{\mathbf{\Omega }}_1}{\min }\Vert {\mathbf{\Omega }}_1{\mathbf{R}}_1-{\mathbf{Z}}_1{\Vert }_F^2+\Vert {\mathbf{\Omega }}_2({\mathbf{\Omega }}_1{\mathbf{R}}_1-{\mathbf{Z}}_1)-{\mathbf{Z}}_2{\Vert }_F^{2}s.t.{\mathbf{\Omega }}_1^T{\mathbf{\Omega }}_1=\mathbf{I}.\text{vspace}-0.05in$$

$$\underset{{\mathbf{Z}}_2}{\min }\Vert {\mathbf{\Omega }}_2({\mathbf{\Omega }}_1{\mathbf{R}}_1-{\mathbf{Z}}_1)-{\mathbf{Z}}_2{\Vert }_F^2+{\eta }_2^2\Vert {\mathbf{Z}}_2{\Vert }_0.$$

$$\underset{{\mathbf{\Omega }}_2}{\min }\Vert {\mathbf{\Omega }}_2({\mathbf{\Omega }}_1{\mathbf{R}}_1-{\mathbf{Z}}_1)-{\mathbf{Z}}_2{\Vert }_F^{2}s.t.{\mathbf{\Omega }}_2^T{\mathbf{\Omega }}_2=\mathbf{I}.\text{vspace}-0.05in$$

$$\underset{\mathbf{x}\ge 0}{\min }\frac{1}{2}\Vert \mathbf{y}-\mathbf{Ax}{\Vert }_{\mathbf{W}}^2+{\mathsf{S}}_2(\mathbf{x}),\text{vspace}-0.05in$$

$$\underset{{\mathbf{Z}}_1}{\min }\Vert {\mathbf{\Omega }}_1\mathbf{X}-{\mathbf{Z}}_1{\Vert }_F^2+{\gamma }_1^2\Vert {\mathbf{Z}}_1{\Vert }_0+\Vert {\mathbf{\Omega }}_2({\mathbf{\Omega }}_1\mathbf{X}-{\mathbf{Z}}_1)-{\mathbf{Z}}_2{\Vert }_F^2.\text{vspace}-0.05in$$

$$\underset{{\mathbf{Z}}_2}{\min }\Vert {\mathbf{\Omega }}_2({\mathbf{\Omega }}_1\mathbf{X}-{\mathbf{Z}}_1)-{\mathbf{Z}}_2{\Vert }_F^2+{\gamma }_2^2\Vert {\mathbf{Z}}_2{\Vert }_0.\text{vspace}-0.05in$$

Full text

Index Terms—  Low-dose CT, statistical image reconstruction, sparse representation, transform learning, unsupervised learning.

1 Introduction

Methods for image reconstruction from limited or corrupted data often exploit various inherent properties or models of the images. A variety of models such as sparsity, tensor, manifold, and convolutional models, etc. [1, 2, 3, 4], have been exploited for explaining or reconstructing images in computational imaging applications. In this work, we focus our investigations on generalization of a subset of models called sparsifying transform models [5, 6], their learning, and application to low-dose computed tomography (CT) image reconstruction.

A major challenge in CT imaging is to reduce the radiation exposure to patients while maintaining the high quality of reconstructed images. This is typically done by reducing the X-ray dose to low or ultra-low levels or by reducing the number of projection views (sparse-view CT). In such cases, conventional filtered back-projection (FBP) [7] reconstructions suffer from artifacts that degrade image quality.

Model-based image reconstruction methods produce accurate reconstructions from reduced dose CT measurements [8]. In particular, penalized weighted-least squares (PWLS) approaches, which have shown promise for CT reconstruction optimize a weighted-least squares data fidelity or measurement modeling term (for the logarithm of the measurements) along with added regularization exploiting prior knowledge of the underlying object [9].

Learning signal models or priors from datasets of images or image patches is an attractive way to obtain adaptive CT image features to improve reconstruction. Recent works have proposed learning various models, including dictionary and sparsifying transform models [10, 11, 12], as well as supervised learning for reconstruction [13]. The learning of sparsifying transform models offers numerous advantages [5] over synthesis dictionary models. In particular, sparse coding in the dictionary model can be expensive, whereas in the sparsifying transform (ST) model, sparse coefficient maps are computed exactly and cheaply by thresholding-type operations (i.e., transform sparse coding even with the $\ell_{0}$ “norm” is not NP-hard). Thus, transform learning-based approaches, including those for image reconstruction, can offer significant computational benefits [12] and often come with convergence guarantees [6, 14, 15]. Recent work has also shown that they can generalize better to unseen data than supervised deep learning schemes [16, 15].

In this work, we investigate the model-based learning of a two-layer extension of the transform model [17] from datasets for image reconstruction. The transform domain or filtering residuals for the data are further sparsified in the second layer. The method in [17] exploited downsampling/pooling operations for image denoising that cannot be readily incorporated in the general inverse problem optimization explored here. Here, we propose pre-learning the two-layer transform (and estimating corresponding sparse coefficient maps) in a model-based fashion to minimize the aggregated transform domain residuals in the second layer, which is used as a regularizer for reconstruction. An efficient block coordinate descent algorithm is derived for learning and for reconstruction with the pre-learned regularizer. Unlike the recent multi-layer convolutional sparse coding (ML-CSC) approach [18, 19, 20], which uses the general synthesis dictionary model and involves expensive sparse coding, exact and cheap sparse coefficients can be computed in our models. Moreover, ML-CSC sparsified the sparse coefficients over layers rather than reducing the modeling residuals. With the transform model, optimizing the residuals significantly improved performance. Finally, ML-CSC has not been investigated for imaging inverse problems. Here, we present numerical experiments demonstrating potential for our approach for low-dose CT reconstruction compared to recent learned single layer transform and nonadaptive methods.

2 Learning and Reconstruction Formulations

This section discusses the proposed two-layer framework and formulations for learning and image reconstruction.

2.1 Two-Layer Residual Transform Learning

For a signal $x\in\mathbb{R}^{N_{p}}$ and operator $W\in\mathbb{R}^{p\times N_{p}}$, the sparsifying transform model suggests that $Wx\approx z$, where $z$ has many zeros. Given the signal $x$ and operator $W$, the transform sparse coding problem finds the best sparse approximation $z$ by minimizing the approximation error or residual in $Wx\approx z$, and the solution is obtained in closed-form by thresholding $Wx$ [5]. When the transform is applied to all the overlapping patches of the image, the model is equivalent to a sparsifying filterbank for images [17, 21].

Here, we study a two-layer extension of the transform model, in which the transform domain residuals or sparse approximation errors in the first layer are further sparsified in the second layer. We propose a patch-based formulation for learning, which could also be equivalently cast in a convolutional form [17]. Given $N^{\prime}$ vectorized (2D or 3D) image patches extracted from a dataset of CT images or volumes, we learn transforms $\{\mathbf{\Omega}_{1},\mathbf{\Omega}_{2}\in\mathbb{R}^{p\times p}\}$ by solving the following training optimization problem:

[TABLE]

where $\mathbf{R}_{1}\in\mathbb{R}^{p\times N^{\prime}}$ denotes the matrix whose columns are the initial vectorized training image patches, $\mathbf{R}_{2}\in\mathbb{R}^{p\times N^{\prime}}$ denotes the residual maps in the second layer, and $\{\mathbf{Z}_{1},\mathbf{Z}_{2}\in\mathbb{R}^{p\times N^{\prime}}\}$ denote sparse coefficient maps in two layers. The non-negative parameters $\{\eta_{l}\}$ control the sparsity of the coefficient maps, with the $\ell_{0}$ “norm” counting the number of non-zero entries in a matrix or vector. The transforms $\{\mathbf{\Omega}_{1},\mathbf{\Omega}_{2}\}$ are assumed to be unitary [6], which simplifies the optimization, and $\mathbf{I}$ denotes the identity matrix.

2.2 CT Image Reconstruction Formulation

We propose using a pre-learned two-layer transform model as a prior for image reconstruction. We reconstruct the image or volume $\mathbf{x}\in\mathbb{R}^{N_{p}}$ from noisy sinogram data $\mathbf{y}\in\mathbb{R}^{N_{d}}$ by solving the following PWLS optimization problem:

[TABLE]

where $\mathbf{A}\in\mathbb{R}^{N_{d}\times N_{p}}$ is the system matrix of the CT scan, $\mathbf{W}=\mathsf{diag}\{w_{i}\}\in\mathbb{R}^{N_{d}\times N_{d}}$ is the diagonal weighting matrix with elements being the estimated inverse variance of $y_{i}$ [9], parameter $\beta>0$ controls the trade-off between noise and resolution, and the regularizer $\mathsf{S}(\mathbf{x})$ based on (P0) is

[TABLE]

Here, $\{\gamma_{l}\}$ are non-negative scalar parameters, the operator $\mathbf{P}^{j}\in\mathbb{R}^{p\times N_{p}}$ extracts the $j$th patch of $p$ voxels of $\mathbf{x}$ as $\mathbf{P}^{j}\mathbf{x}$, and $\mathbf{r}^{j}_{1}$ denotes the $j$th column of $\mathbf{R}_{1}$. The columns of $\mathbf{Z}_{l}$ are $\{\mathbf{z}_{l}^{j}\in\mathbb{R}^{p}\}_{j=1}^{N_{r}}$ and denote the transform-sparse coefficients in the $l$th layer, where $N_{r}$ is the number of extracted patches.

3 Algorithms

3.1 Algorithm for Learning

We solve (P0) using an exact block coordinate descent algorithm that alternates between sparse coding steps (solving for $\mathbf{Z}_{1}$ or $\mathbf{Z}_{2}$) and transform update steps (solving for $\mathbf{\Omega}_{1}$ or $\mathbf{\Omega}_{2}$). The transforms $\mathbf{\Omega}_{1}$ and $\mathbf{\Omega}_{2}$ and the coefficients $\mathbf{Z}_{2}$ need to be first initialized. In our experiments, we used the 2D DCT and identity matrices to initialize $\mathbf{\Omega}_{1}$ and $\mathbf{\Omega}_{2}$ respectively, and the initial $\mathbf{Z}_{2}$ was an all-zero matrix.

3.1.1 Sparse Coding Step for $\mathbf{Z}_{1}$

Here, we solve the following sub-problem for $\mathbf{Z}_{1}$ with all other variables fixed:

[TABLE]

Substituting $\mathbf{R}_{2}=\mathbf{\Omega}_{1}\mathbf{R}_{1}-\mathbf{Z}_{1}$ and using the unitary property of $\mathbf{\Omega}_{2}$, we rewrite (2) as $\min_{\mathbf{Z}_{1}}2\|\mathbf{Z}_{1}-(\mathbf{\Omega}_{1}\mathbf{R}_{1}-0.5\mathbf{\Omega}_{2}^{T}\mathbf{Z}_{2})\|^{2}_{F}+\eta_{1}^{2}\|\mathbf{Z}_{1}\|_{0}$. Then the optimal solution is obtained as $\hat{\mathbf{Z}}_{1}=H_{\eta_{1}/\sqrt{2}}(\mathbf{\Omega}_{1}\mathbf{R}_{1}-0.5\mathbf{\Omega}_{2}^{T}\mathbf{Z}_{2})$, where the hard-thresholding operator $H_{\eta}(\cdot)$ zeros out vector entries with magnitude less than $\eta$.

3.1.2 Transform Update Step for $\mathbf{\Omega}_{1}$

With $\mathbf{\Omega}_{2}$, $\mathbf{Z}_{2}$, and $\mathbf{Z}_{1}$ fixed, we update $\mathbf{\Omega}_{1}$ by solving the following sub-problem:

[TABLE]

This is equivalent to minimizing the cost $2||\mathbf{\Omega}_{1}\mathbf{R}_{1}-(\mathbf{Z}_{1}+0.5\mathbf{\Omega}_{2}^{T}\mathbf{Z}_{2})||_{F}^{2}$. Denoting the full singular value decomposition (SVD) of $\mathbf{R}_{1}\mathbf{Z}_{1}^{T}+0.5\mathbf{R}_{1}\mathbf{Z}_{2}^{T}\mathbf{\Omega}_{2}$ as $\mathbf{U}_{1}\mathbf{\Sigma}_{1}\mathbf{V}_{1}^{T}$ (cf. [6]), the optimal solution to (3) is $\hat{\mathbf{\Omega}}_{1}=\mathbf{V}_{1}\mathbf{U}_{1}^{T}$.

3.1.3 Sparse Coding Step for $\mathbf{Z}_{2}$

With $\mathbf{\Omega}_{2}$, $\mathbf{Z}_{1}$, and $\mathbf{\Omega}_{1}$ fixed, we update $\mathbf{Z}_{2}$ by solving the following sub-problem:

[TABLE]

The optimal sparse coefficients for the second layer are readily computed in closed-form by hard-thresholding as $\hat{\mathbf{Z}}_{2}=H_{\eta_{2}}(\mathbf{\Omega}_{2}(\mathbf{\Omega}_{1}\mathbf{R}_{1}-\mathbf{Z}_{1}))$.

3.1.4 Transform Update Step for $\mathbf{\Omega}_{2}$

Here, we update $\mathbf{\Omega}_{2}$ keeping the other variables fixed by solving:

[TABLE]

Denoting the full SVD of $(\mathbf{\Omega}_{1}\mathbf{R}_{1}-\mathbf{Z}_{1})\mathbf{Z}_{2}^{T}$ as $\mathbf{U}_{2}\mathbf{\Sigma}_{2}\mathbf{V}_{2}^{T}$, the optimal solution to (5) is $\hat{\mathbf{\Omega}}_{2}=\mathbf{V}_{2}\mathbf{U}_{2}^{T}$.

3.2 Image Reconstruction Algorithm

We propose an alternating-type algorithm for (P1) that alternates between updating $\mathbf{x}$ (image update step), and $\mathbf{Z}_{1}$ and $\mathbf{Z}_{2}$ (sparse coding steps).

3.2.1 Image Update Step

With the variables $\mathbf{Z}_{1}$ and $\mathbf{Z}_{2}$ fixed, we solve (P1) for $\mathbf{x}$, which reduces to the following weighted least squares problem:

[TABLE]

where \mathsf{S}_{2}(\mathbf{x})\triangleq\beta\sum_{j=1}^{N_{r}}\big{\{}\|\mathbf{\Omega}_{2}(\mathbf{\Omega}_{1}\mathbf{P}^{j}\mathbf{x}-\mathbf{z}_{1}^{j})-\mathbf{z}_{2}^{j}\|^{2}_{2}+\|\mathbf{\Omega}_{1}\mathbf{P}^{j}\mathbf{x}-\mathbf{z}_{1}^{j}\|^{2}_{2}\big{\}}. We solve (6) using the efficient relaxed OS-LALM algorithm [22]. The algorithmic details are similar to those in [12]. We precompute a diagonal majorizing matrix of the Hessian of the regularizer $\mathsf{S}_{2}(\mathbf{x})$ as $\mathbf{D}_{\mathsf{S}_{2}}\triangleq\nabla^{2}\mathsf{S}_{2}(\mathbf{x})=4\beta\sum_{j=1}^{N_{r}}(\mathbf{P}^{j})^{T}\mathbf{P}^{j}$.

3.2.2 Sparse Coding Steps

First, with $\mathbf{x}$ and $\mathbf{Z}_{2}$ fixed and $\mathbf{X}$ denoting the matrix with $\mathbf{P}^{j}\mathbf{x}$ as its columns, we update $\mathbf{Z}_{1}$ by solving

[TABLE]

Similar to the solution for (2), the optimal solution for (7) is $\hat{\mathbf{Z}}_{1}=H_{\gamma_{1}/\sqrt{2}}(\mathbf{\Omega}_{1}\mathbf{X}-0.5\mathbf{\Omega}_{2}^{T}\mathbf{Z}_{2})$.

Next, with $\mathbf{X}$ and $\mathbf{Z}_{1}$ fixed, coefficients $\mathbf{Z}_{2}$ are updated by solving the following sub-problem:

[TABLE]

The optimal solution is $\hat{\mathbf{Z}}_{2}=H_{\gamma_{2}}(\mathbf{\Omega}_{2}(\mathbf{\Omega}_{1}\mathbf{X}-\mathbf{Z}_{1}))$.

4 Experimental Results

We evaluated the proposed PWLS reconstruction method with a two-layer learned regularizer (referred to as PWLS-MRST2) and compared its image reconstruction quality with those of the FBP method with a Hanning window, and the PWLS-EP method that uses a non-adaptive edge-preserving regularizer $\mathsf{S}(\mathbf{x})=\sum_{j=1}^{N_{p}}\sum_{k\in N_{j}}\kappa_{j}\kappa_{k}\varphi(x_{j}-x_{k})$, where $N_{j}$ is the size of the neighborhood, $\kappa_{j}$ and $\kappa_{k}$ are the parameters encouraging uniform noise [23], and ${\varphi}{(t)}\triangleq\delta^{2}(|t/\delta|-\log(1+|t/\delta|))$ with $\delta=10$ Hounsfield units (HU)111Modified HU is used, where air is [math] HU and water is $1000$ HU.. We optimized the PWLS-EP problem using the relaxed OS-LALM method [22]. We also compared to the previous PWLS-ST method that uses a learned single-layer (square) transform [12, 11].

Various methods are compared quantitatively using the Root Mean Square Error (RMSE) and Peak Signal to Noise Ratio (PSNR) metrics in a region of interest (ROI). The RMSE of the reconstruction $\hat{\mathbf{x}}$ is defined as RMSE $=\sqrt{\sum_{i=1}^{N_{p}}(\hat{x}_{i}-x^{*}_{i})^{2}/{N_{p}}}$, where $\mathbf{x}^{*}$ is the ground truth image and $N_{p}$ is the number of pixels (voxels) in the ROI. We tuned the parameters of various methods for each experiment to achieve the lowest RMSE and highest PSNR.

We pre-learned two transforms (Fig. 1) for the proposed two-layer model from $8\times 8$ image patches extracted from five $512\times 512$ XCAT phantom [24] slices, with $\eta_{1}=80$, $\eta_{2}=60$, and a patch extraction stride $1\times 1$. We ran $1000$ iterations of the learning algorithm in Section 3.1 to ensure convergence. We simulated 2D fan-beam CT test scans using $840\times 840$ XCAT phantom slices (air cropped) that differ from the training slices, with pixel size $\Delta_{x}=\Delta_{y}=0.4883$ mm. Noisy sinograms of size $888\times 984$ were numerically simulated with GE LightSpeed fan-beam geometry corresponding to a monoenergetic source with $10000$, $5000$, and $3000$ incident photons per ray and no scatter, respectively. We reconstructed two $420\times 420$ images with a coarser grid, where $\Delta_{x}=\Delta_{y}=0.9766$ mm. The ROI here was a circular (around center) region containing all the phantom tissues.

Initialized with FBP reconstructions, we ran the PWLS-EP algorithm for $50$ iterations using relaxed OS-LALM with $24$ subsets. The PWLS-EP result was used to initialize the adaptive methods. The parameters for different methods for $I_{0}=10000$, $5000$, and $3000$ are as follows: $\beta=2^{16},2^{16.5}$, and $2^{16.5}$ respectively, for Slice $1$ and $\beta=2^{16}$ for Slice $2$ for PWLS-EP; $(\beta,\gamma_{1})=\left(2\times 10^{5},20\right)$ , $\left(1.3\times 10^{5},20\right)$, and $\left(1.3\times 10^{5},20\right)$ respectively for Slice $1$ and $\left(2.2\times 10^{5},20\right)$, $\left(2\times 10^{5},20\right)$, and $\left(1.5\times 10^{5},20\right)$ for Slice $2$ for PWLS-ST; and $(\beta,\gamma_{1},\gamma_{2})=\left(9\times 10^{4},30,10\right)$, $\left(4\times 10^{4},30,12\right)$ and $\left(3.5\times 10^{4},30,12\right)$, respectively for Slice $1$ and $\left(8\times 10^{4},30,12\right)$, $\left(5\times 10^{4},30,12\right)$ and $\left(5\times 10^{4},30,7\right)$ for Slice $2$ for PWLS-MRST2. For PWLS-ST and PWLS-MRST2, the image reconstruction algorithms were run for $1000$ and $1500$ outer iterations with $4$ and $2$ ordered subsets, respectively, and $2$ inner iterations of the image update step that ensured convergence.

Table 1 summarizes the RMSE and PSNR values for reconstructions with FBP, PWLS-EP, PWLS-ST, and the proposed PWLS-MRST2 for the three tested photon intensities. The adaptive PWLS methods significantly outperform the conventional FBP and the non-adaptive PWLS-EP. Moreover, PWLS-MRST2 with a learned two-layer model improves the reconstruction quality over the single-layer PWLS-ST scheme. It differs from PWLS-ST by only an additional simple sparse coding step and thus has a similar computational cost.

Fig. 2 shows representative reconstructions for FBP, PWLS-EP, PWLS-ST, and PWLS-MRST2. Compared to FBP and PWLS-EP, PWLS-MRST2 significantly improves image quality by reducing noise and preserving structural details. Furthermore, PWLS-MRST2 improves the quality of the central region and image edges compared to PWLS-ST.

5 Conclusion

We presented the learning of a two-layer extension of the sparsifying transform model for CT image reconstruction from low-dose measurements. The model is learned from datasets to sparsify the filtering or transform domain residuals in the second layer. The algorithms for both learning and reconstruction derived for the simple two-layer case are block coordinate descent-type algorithms and involve efficient updates. Our experimental results illustrated the superior performance of a learned two-layer scheme over the single layer adaptive transform scheme. The learned approaches significantly outperformed nonadaptive methods. Since unsupervised model-based learning of deep models for imaging is a new area, we plan to investigate the learning of more complex models and more layers for CT image reconstruction and other tasks in future work.

Acknowledgments

The authors thank Xikai Yang for his help with part of the experiments.