---
abstract: 'Signal models based on sparse representation have received considerable attention in recent years. Compared to synthesis dictionary learning, sparsifying transform learning involves highly efficient sparse coding and operator update steps. In this work, we propose a Multi-layer Residual Sparsifying Transform (MRST) learning model wherein the transform domain residuals are jointly sparsified over layers. In particular, the transforms for the deeper layers exploit the more intricate properties of the residual maps. We investigate the application of the learned MRST model for low-dose CT reconstruction using Penalized Weighted Least Squares (PWLS) optimization. Experimental results on Mayo Clinic data show that the MRST model outperforms conventional methods such as FBP and PWLS methods based on edge-preserving (EP) regularizer and single-layer transform (ST) model, especially for maintaining some subtle details.'
author:
- 'Xikai Yang, Xuehang Zheng, Yong Long$^\star$, Saiprasad Ravishankar, [^1] [^2] [^3]'
bibliography:
- 'refs.bib'
title: 'Learned Multi-layer Residual Sparsifying Transform Model for Low-dose CT Reconstruction'
---

Low-dose CT, Statistical image reconstruction, Sparse representation, Transform learning, Unsupervised learning.

Introduction
============

Signal models exploiting sparsity have been shown to be useful in a variety of applications such as compression, restoration, denoising, reconstruction, etc. Natural signals can be modeled as sparse in a synthesis dictionary (i.e., represented as linear combinations of few dictionary atoms or columns) or in a sparsifying transform domain.

Transforms such as wavelets [@pati:93:omp] and the discrete cosine transform (DCT) are well-known to sparsify images. Methods such as synthesis dictionary learning [@aharon:06:ksa] and analysis dictionary learning [@rubinstein:13:aka] allow adapting the model to the data. These dictionary learning problems are NP-hard in general and involve computationally expensive alternating-type algorithms that limit applicability to large-scale data. In contrast, the recently proposed sparsifying transform learning approaches [@ravishankar:13:lst] involve exact and highly efficient updates. In particular, the transform model approximates the signal as sparse in a transformed domain. Adopting a multi-layer sparsifying transform model enables sparsifying an input image successively over layers [@ravishankar:18:lml], creating a rich and more complete sparsity model, which forms the core of this work.

One of the most important application of such image models is for medical image reconstruction. In particular, an important problem in computed tomography (CT) is reducing the X-ray exposure to patients while maintaining good image reconstruction quality. A conventional method for CT reconstruction is the analytical filtered back-projection (FBP) [@feldkamp:84:pcb]. However, image quality degrades severely for FBP when the radiation dose is reduced. In contrast, model-based image reconstruction (MBIR) achieves comparable or better image quality [@elbakri:02:sir].

A typical MBIR method for low-dose CT (LDCT) is the penalized weighted least squares (PWLS) approach. The cost function for PWLS includes a weighted quadratic data-fidelity term and a penalty term or regularizer capturing prior information or model of the object [@sauer:93:alu; @thibault:07:atd]. Recent works have shown promising LDCT reconstruction quality by incorporating data-driven models into the regularizer, where the models are learned from datasets of images or image patches. In particular, PWLS reconstruction with adaptive sparsifying transform-based regularization has shown promise for tomographic reconstruction [@pfister:14:mbi; @zheng:18:pua; @chun:17:esv]. Recent work has also shown that they may generalize better to new data than supervised deep learning schemes [@ye:19:sld]. The adaptive transform-based image reconstruction algorithms can exploit a variety of image models [@pfister:14:mbi; @zheng:18:pua; @zhou:13:atf] learned in an unsupervised manner, and involve efficient closed-form solutions for sparse coding.

In this work, we propose a new formulation and algorithm for learning a multi-layer transform model [@ravishankar:18:lml], where the transform domain residuals (the difference between transformed data and their sparse approximations) are successively sparsified over several layers. We refer to the model as a multi-layer residual sparsifying transform (MRST) model. The transforms are learned over several layers from images to jointly minimize the transform domain residuals across layers, while enforcing sparsity conditions in each layer. Importantly, the filters in deeper layers can help better exploit finer features (e.g., edges and correlations) in the residual maps. We investigate the performance of the (unsupervised) learned MRST model for LDCT reconstruction using PWLS. We propose efficient alternating minimization algorithms for both learning and reconstruction. Experimental results on Mayo Clinic data illustrate that the learned MRST model outperforms the conventional FBP as well as PWLS methods based on the non-adaptive edge-preserving (EP) regularizer and learned single-layer transform model, especially for maintaining some subtle details. In the following sections, we will first discuss the proposed model and formulation, followed by our algorithms and experimental results on Mayo Clinic data.

Problem Formulation
===================

Here, we will introduce the proposed general multi-layer learning framework and the formulation for LDCT image reconstruction. The MRST learning cost and constraints are shown in Problem , which is an extension of simple single-layer transform learning. $$\label{eq:P0}
\begin{aligned}
& \quad \quad \quad \quad \min_{\{{\mathbf{\Omega}}_l,{\mathbf{Z}}_l\}}   \sum_{l=1}^L \bigg\{  \|{\mathbf{\Omega}}_l{\mathbf{R}}_l - {\mathbf{Z}}_l\|_F^2 + \eta_l^2\|{\mathbf{Z}}_l\|_0  \bigg\}
\\
& \quad \mathrm{s.t.} \quad {\mathbf{R}}_l = {\mathbf{\Omega}}_{l-1}{\mathbf{R}}_{l-1} - {\mathbf{Z}}_{l-1}, 2\leq l\leq L, {\mathbf{\Omega}}_l^T {\mathbf{\Omega}}_l  = {\mathbf{I}}, \forall l.
\end{aligned}
\tag{P0}$$ Here, $\{{\mathbf{\Omega}}_l\in \mathbb{R}^{p \times p}\}$ and $\{{\mathbf{Z}}_l\in \mathbb{R}^{p \times N}\}$ denote the learned transforms and sparse coefficient maps for the $1 \leq l \leq L$ layers. Parameter $\eta_l$ controls the maximum allowed sparsity level (computed using the $\ell_0$ “norm” penalty) of ${\mathbf{Z}}_l$. The residual maps $\{{\mathbf{R}}_l\in \mathbb{R}^{p \times N}\}$ are defined in recursive form over layers, with ${\mathbf{R}}_1$ denoting the input training data. Here, we assume $\mathbf{R}_1$ to be a matrix, whose columns are (vectorized) patches drawn from image data sets. The unitary constraints for $\{{\mathbf{\Omega}}_l\}$ enable closed-form solutions for the sparse coefficient and transform update steps in our algorithms. The learned MRST model can then be used to construct a data-driven regularizer in PWLS as shown in Problem . $$\label{eq:P1}   
\min_{{\mathbf{x}}\geq {\mathbf{0}}}  \frac{1}{2}\|{\mathbf{y}}- {\mathbf{A}}{\mathbf{x}}\|^2_{{\mathbf{W}}}  + \beta {\mathsf{S}}({\mathbf{x}}),   
\tag{P1}$$ $$\label{eq:Rx_DST}
\begin{aligned}
\quad \quad \quad &{\mathsf{S}}({\mathbf{x}}) \triangleq  \min_{\{{\mathbf{Z}}_l\}}  \sum_{l=1}^{L}  \bigg\{  \|{\mathbf{\Omega}}_l {\mathbf{R}}_l - {\mathbf{Z}}_l\|^2_F + \gamma_l^2\|{\mathbf{Z}}_l\|_0   \bigg\}  
\\ &\;   \mathrm{s.t.} \;  {\mathbf{R}}_l = {\mathbf{\Omega}}_{l-1}{\mathbf{R}}_{l-1} - {\mathbf{Z}}_{l-1}, 2 \leq l \leq L \, ,  {\mathbf{R}}_1^j = {\mathbf{P}}^j{\mathbf{x}}, \, \forall \, j.
\end{aligned}
\vspace{-0.05in}$$ In particular, we reconstruct the image ${\mathbf{x}}\in \mathbb{R}^{N_p}$ from noisy sinogram data ${\mathbf{y}}\in \mathbb{R}^{N_d}$ by solving , where ${\mathbf{A}}\in  \mathbb{R}^{  N_d  \times N_p}$ is the system matrix of the CT scan and ${\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$. Operator ${\mathbf{P}}^j\in \mathbb{R}^{p\times N_p}$ extracts the $j$th patch of ${\mathbf{x}}$ as ${\mathbf{P}}^j{\mathbf{x}}$. The $j$th colums of ${\mathbf{R}}_l$ and ${\mathbf{Z}}_l$ are denoted ${\mathbf{Z}}_l^j$ and ${\mathbf{R}}_l^j$. The non-negative parameters $\{\gamma_l\}$ control the sparsity of the coefficient maps in different layers, and $\beta>0$ captures the relative trade-off between the data-fidelity term and regularizer.

Algorithms for Learning and Reconstruction
==========================================

Transform Learning Algorithm
----------------------------

We propose an exact block coordinate descent (BCD) algorithm for the nonconvex Problem that cycles over updating each ${\mathbf{Z}}_l$ for $1 \leq l \leq L$ (sparse coding step) and each ${\mathbf{\Omega}}_l$ for $1 \leq l \leq L$ (transform update step). In each step, the remainder of the variables (that are not optimized) are kept fixed. The exact BCD algorithm ensures that the objective in is monotone decreasing over iterations and that it converges. In particular, under the unitarity condition on the transforms, every subproblem in the block coordinate descent minimization approach can be solved exactly. We initialize the algorithm with 2D DCT for ${\mathbf{\Omega}}_1$ and the identity matrix for $\{{\mathbf{\Omega}}_l\}  (l\geq2)$ respectively. The initial $\{{\mathbf{Z}}_l\}$ are all-zero matrices. Since the residuals are defined recursively in , to simplify the algorithmic description, we first define matrices ${\mathbf{B}}_p^q (p < q)$, which can be regarded as backpropagation matrices from the $q$th to $p$th layers. $$\label{eq:B_pq}
\begin{aligned}
&{\mathbf{B}}_p^q={\mathbf{\Omega}}_{p+1}^T{\mathbf{Z}}_{p+1}+{\mathbf{\Omega}}_{p+1}^T{\mathbf{\Omega}}_{p+2}^T{\mathbf{Z}}_{p+2}+...+{\mathbf{\Omega}}_{p+1}^T{\mathbf{\Omega}}_{p+2}^T...{\mathbf{\Omega}}_q^T{\mathbf{Z}}_q
\\
&=\sum_{k=p+1}^q \bigg( \prod_{s=p+1}^k {\mathbf{\Omega}}_s^T \bigg){\mathbf{Z}}_k
\end{aligned}$$

### Sparse Coding Step for ${\mathbf{Z}}_l$

Here, we solve (P0) for ${\mathbf{Z}}_l$ with all other variables fixed. The corresponding subproblem is as in . $$\label{eq:sub_pro_Z}
\vspace{-0.05in}
\min_{{\mathbf{Z}}_l}   \sum_{i=l}^L \bigg\{  \|{\mathbf{\Omega}}_i{\mathbf{R}}_i - {\mathbf{Z}}_i\|_F^2 \bigg\}  + \eta_l^2\|{\mathbf{Z}}_l\|_0$$ Although the cost function is nonconvex and the residual maps also depend on ${\mathbf{Z}}_l$ (as in ), we can exploit the unitarity of the transforms to simplify and rewrite Problem as $\min_{{\mathbf{Z}}_l}(L-l+1)\times\|{\mathbf{Z}}_l - ( {\mathbf{\Omega}}_l{\mathbf{R}}_l - \frac{1}{L-l+1}\sum_{i=l+1}^L {\mathbf{B}}_l^i )\|_F^2 + \eta_{l}^{2}\|{\mathbf{Z}}_l\|_0$. This problem has a similar form as the single-transform sparse coding problem [@ravishankar:13:lst], and the optimal solution $\hat{{\mathbf{Z}}}_l$ is obtained as in , where $H_{\eta}(\cdot)$ denotes the *hard-thresholding* operator that sets elements with magnitude less than the threshold $\eta$ to zero. $$\label{eq:Z_l}
\begin{aligned}
\hat{{\mathbf{Z}}}_l=
&\begin{cases}
H_{\eta_l/\sqrt{L-l+1}} \bigg( {\mathbf{\Omega}}_l{\mathbf{R}}_l - \frac{1}{L-l+1}\sum_{i=l+1}^L {\mathbf{B}}_l^i \bigg),&  1\leq l \leq L-1,\\
H_{\eta_L} ({\mathbf{\Omega}}_L{\mathbf{R}}_L),& l=L.
\end{cases}
\end{aligned}$$

### Transform Update Step for ${\mathbf{\Omega}}_l$

Here, we fix $\{{\mathbf{Z}}_l\}$ and all ${\mathbf{\Omega}}_l$ (except the target ${\mathbf{\Omega}}_l$ in ) and solve the following subproblem: $$\label{eq:sub_pro_mOmega}
\min_{{\mathbf{\Omega}}_l}   \sum_{i=l}^L \bigg\{ \|{\mathbf{\Omega}}_i{\mathbf{R}}_i - {\mathbf{Z}}_i\|_F^2 \bigg\}
\quad \mathrm{s.t.} \quad  {\mathbf{\Omega}}_l^T {\mathbf{\Omega}}_l  = {\mathbf{I}}.$$

Incorporating the recursive dependence of the residual maps on the transforms (as in ) along with the unitarity of the transforms simplifies the solution to . First, denoting the full singular value decomposition (SVD) of the matrix $\mathbf{G}_l$ below by ${\mathbf{U}}_l{\mathbf{\Sigma}}_l{\mathbf{V}}_l^T$, the optimal solution to is $\hat{{\mathbf{\Omega}}_l}={\mathbf{V}}_l{\mathbf{U}}_l^T$. $$\label{eq:mOmega_l}
\mathbf{G}_l =
\begin{cases}
{\mathbf{R}}_l\bigg( {\mathbf{Z}}_l + \frac{1}{L-l+1} \sum_{i=l+1}^L {\mathbf{B}}_l^i \bigg)^T, & 1\leq l \leq L-1,\\
{\mathbf{R}}_L{\mathbf{Z}}_L^T, & l=L.
\end{cases}
\vspace{-0.05in}$$

Image Reconstruction Algorithm
------------------------------

The proposed PWLS-MRST algorithm for low-dose CT image reconstruction exploits the learned MRST model. We reconstruct the image by solving the PWLS problem . We propose a block coordinate descent (BCD) algorithm for that cycles over updating the image ${\mathbf{x}}$ and each of the sparse coefficient maps ${\mathbf{Z}}_l$ for $1 \leq l \leq L$.

### Image Update Step for ${\mathbf{x}}$

First, with the sparse coefficient maps $\{{\mathbf{Z}}_l\}$ fixed, we optimize for ${\mathbf{x}}$ in by optimizing the following subproblem: $$\label{eq:sub_P1}
\min_{{\mathbf{x}}\geq {\mathbf{0}}} \frac{1}{2} \|{\mathbf{y}}- {\mathbf{A}}{\mathbf{x}}\|^2_{{\mathbf{W}}} + \beta {\mathsf{S}}_2({\mathbf{x}}),$$ where ${\mathsf{S}}_2({\mathbf{x}}) \triangleq   \sum_{l=1}^{L}  \bigg\{ \|{\mathbf{\Omega}}_l{\mathbf{R}}_l - {\mathbf{Z}}_l\|_F^2 \bigg\} $, with ${\mathbf{R}}_l={\mathbf{\Omega}}_{l-1}{\mathbf{R}}_{l-1} - {\mathbf{Z}}_{l-1}$, $2 \leq l \leq L$, and ${\mathbf{R}}_1^j = {\mathbf{P}}^j {\mathbf{x}}$. We use the *efficient* relaxed linearized augmented Lagrangian method with ordered-subsets (relaxed OS-LALM) algorithm [@nien:16:rla] to obtain the solution to . The algorithmic details are shown in **Algorithm 1**. In each iteration of the OS-LALM method, we update the image $M$ times corresponding to $M$ ordered subsets. The matrices ${\mathbf{A}}_m$ and ${\mathbf{W}}_m$, and the vector ${\mathbf{y}}_m$ are sub-matrices of ${\mathbf{A}}$ and ${\mathbf{W}}$, and sub-vector of ${\mathbf{y}}$ for the $m$th subset respectively. Matrix ${\mathbf{D}}_{{\mathbf{A}}}$ represents a diagonal majorizing matrix of ${\mathbf{A}}^T{\mathbf{W}}{\mathbf{A}}$. We precompute the Hessian matrix of ${\mathsf{S}}_2({\mathbf{x}})$ as ${\mathbf{D}}_{{\mathsf{S}}_2}$ in to accelerate the algorithm, and the gradient of ${\mathsf{S}}_2({\mathbf{x}})$ is as in . $$\label{eq:grad_R}       
\nabla {\mathsf{S}}_2({\mathbf{x}})= 2 \beta \sum_{j=1}^N ({\mathbf{P}}^T)^j \bigg\{ L{\mathbf{P}}^j{\mathbf{x}}- \sum_{k=1}^L ({\mathbf{B}}_0^k)^j \bigg\}$$ $$\label{eq:D_R}      
{\mathbf{D}}_{{\mathsf{S}}_2}  \triangleq  \nabla^2 {\mathsf{S}}_2({\mathbf{x}}) = 2L \beta \sum_{j=1}^{N}  ({\mathbf{P}}^T)^j{\mathbf{P}}^j \\
\vspace{-0.05in}$$

### Sparse Coding Step for ${\mathbf{Z}}_l$

Similar to the sparse coding step during transform learning, the solution of with respect to each sparse coefficient map ${\mathbf{Z}}_l$ is shown in , and is the solution of . $$\label{eq:recon_pro_z_l}
\begin{aligned}
\min_{{\mathbf{Z}}_l} \sum_{i=l}^L \bigg\{  \|{\mathbf{\Omega}}_i{\mathbf{R}}_i - {\mathbf{Z}}_i\|_F^2 \bigg\}  + \gamma_l^2\|{\mathbf{Z}}_l\|_0  
\\
\quad \mathrm{s.t.} \quad {\mathbf{R}}_i={\mathbf{\Omega}}_{i-1}{\mathbf{R}}_{i-1} - {\mathbf{Z}}_{i-1}, l \leq i \leq L
\end{aligned}$$ $$\label{eq:recon_formula_z_l}
\hat{{\mathbf{Z}}}_l=H_{\eta_l/\sqrt{L-l+1}} \bigg\{ {\mathbf{\Omega}}_l{\mathbf{R}}_l - \frac{1}{L-l+1}\sum_{i=l+1}^L {\mathbf{B}}_l^i \bigg\}$$

**Input:** initial image $\tilde{{\mathbf{x}}}^{(0)}$, all-zero initial $\{\tilde{{\mathbf{Z}}_l}^{(0)}\}$, pre-learned $\{{\mathbf{\Omega}}_l\}$, thresholds $\{\gamma_l\}$, $\alpha = 1.999$, ${\mathbf{D}}_{{\mathbf{A}}} $, ${\mathbf{D}}_{{\mathsf{S}}_2}$, number of outer iterations $T_{O}$, number of inner iterations $N$, and number of subsets $M$. **Output:** reconstructed image $\tilde{{\mathbf{x}}}^{(T)}$.

**1) Image Update**: $\{\tilde{{\mathbf{Z}}_l}^{(t)}\}$ fixed,

**Initialization:** $\rho=1$, ${\mathbf{x}}^{(0)} = \tilde{{\mathbf{x}}}^{(t)}$, ${\mathbf{g}}^{(0)} = {\bm{\zeta}}^{(0)}  = M{\mathbf{A}}_M^{T}{\mathbf{W}}_M({\mathbf{A}}_M{\mathbf{x}}^{(0)}-{\mathbf{y}}_M) $ and ${\mathbf{h}}^{(0)} = {\mathbf{D}}_{\mathbf{A}}{\mathbf{x}}^{(0)} - {\bm{\zeta}}^{(0)}$.

$r = nM +m$ $$\left\{         
        \begin{aligned}
        {\mathbf{s}}^{(r+1)} &= \rho({\mathbf{D}}_{\mathbf{A}}{\mathbf{x}}^{(r)} -{\mathbf{h}}^{(r)}) + (1-\rho){\mathbf{g}}^{(r)} \\
        {\mathbf{x}}^{(r+1)} &= [{\mathbf{x}}^{(r)} - (\rho{\mathbf{D}}_{\mathbf{A}}+{\mathbf{D}}_{{\mathsf{S}}_2})^{-1}({\mathbf{s}}^{(r+1)} +\nabla {\mathsf{S}}_2({\mathbf{x}}^{(r)}))]_+ \\
        {\bm{\zeta}}^{(r+1)} & \triangleq  M {\mathbf{A}}^T_m{\mathbf{W}}_m({\mathbf{A}}_m{\mathbf{x}}^{(r+1)}-{\mathbf{y}}_m)   \\
        {\mathbf{g}}^{(r+1)} &= \frac{\rho}{\rho+1}(\alpha {\bm{\zeta}}^{(r+1)} + (1-\alpha){\mathbf{g}}^{(r)}) +  \frac{1}{\rho+1}{\mathbf{g}}^{(r)}\\
        {\mathbf{h}}^{(r+1)}  &= \alpha({\mathbf{D}}_{{\mathbf{A}}} {\mathbf{x}}^{(r+1)} -{\bm{\zeta}}^{(r+1)}) + (1-\alpha){\mathbf{h}}^{(r)} 
        \end{aligned}
        \right.$$ decreasing $\rho$. $\tilde{{\mathbf{x}}}^{(t+1)} = {\mathbf{x}}^{(NM)}$. **2) Sparse Coding**: with $\tilde{{\mathbf{x}}}^{(t+1)}$ fixed, for each $1\leq l \leq L$, update $\tilde{{\mathbf{Z}}_l}^{(t+1)}$ by .

Experimental Results
====================

In this section, we evaluate the image reconstruction quality for the proposed PWLS-MRST algorithm and compare it with conventional **FBP** and the **PWLS-EP** (Edge Preserving) method [@cho:15:rdf].

Transform Learning
------------------

To better understand the potential of the MRST model, we vary the number of layers and pre-learn transforms for ST, MRST2, MRST3, MRST5, and MRST7, which possess $1$, $2$, $3$, $5$, and $7$ layers, respectively. We used 7 slices of the Mayo Clinic data to train the models. For each model, we run 1000 to 2000 iterations of the learning algorithm to ensure convergence. Fig. \[fig:lear\_tran\] shows some of the learned transforms, with each transform matrix row displayed as a square patch for simplicity. The single layer transform displays edge-like and directional structures that sparsify the image. However, with more layers, finer level features are learned to sparsify transform domain residuals in deeper layers. Nonetheless, transforms in deep layers could be more easily contaminated with noise in the training data, since the main image features are successively filtered out over layers.

Image Reconstruction
--------------------

\[spy using outlines=[rectangle,green,magnification=2,size=7mm, connect spies]{}\] ; on (-1.35,0.85) in node \[left\] at (-1.58,1.93); on (-0.80,-0.5) in node \[left\] at (-1.58,-1.9); on (-0.80,0.40) in node \[left\] at (2.25,1.92);

(-90,9)[ ]{}

\[spy using outlines=[rectangle,green,magnification=2,size=7mm, connect spies]{}\] ; on (-1.35,0.85) in node \[left\] at (-1.58,1.93); on (-0.80,-0.5) in node \[left\] at (-1.58,-1.9); on (-0.80,0.40) in node \[left\] at (2.25,1.92);

(-78,9)[ ]{}

\[spy using outlines=[rectangle,green,magnification=2,size=7mm, connect spies]{}\] ; on (-1.35,0.85) in node \[left\] at (-1.58,1.93); on (-0.80,-0.5) in node \[left\] at (-1.58,-1.9); on (-0.80,0.40) in node \[left\] at (2.25,1.92);

(-75,9)[ ]{}

\[spy using outlines=[rectangle,green,magnification=2,size=7mm, connect spies]{}\] ; on (-1.35,0.85) in node \[left\] at (-1.58,1.93); on (-0.80,-0.5) in node \[left\] at (-1.58,-1.9); on (-0.80,0.40) in node \[left\] at (2.25,1.92);

(-75,9)[ ]{}\

\[spy using outlines=[rectangle,green,magnification=2,size=7mm, connect spies]{}\] ; on (-1.35,0.85) in node \[left\] at (-1.58,1.93); on (-0.80,-0.5) in node \[left\] at (-1.58,-1.9); on (-0.80,0.40) in node \[left\] at (2.25,1.92);

(-85,9)[ ]{}

\[spy using outlines=[rectangle,green,magnification=2,size=7mm, connect spies]{}\] ; on (-1.35,0.85) in node \[left\] at (-1.58,1.93); on (-0.80,-0.5) in node \[left\] at (-1.58,-1.9); on (-0.80,0.40) in node \[left\] at (2.25,1.92);

(-85,9)[ ]{}

\[spy using outlines=[rectangle,green,magnification=2,size=7mm, connect spies]{}\] ; on (-1.35,0.85) in node \[left\] at (-1.58,1.93); on (-0.80,-0.5) in node \[left\] at (-1.58,-1.9); on (-0.80,0.40) in node \[left\] at (2.25,1.92);

(-85,9)[ ]{}

\[spy using outlines=[rectangle,green,magnification=2,size=7mm, connect spies]{}\] ; on (-1.35,0.85) in node \[left\] at (-1.58,1.93); on (-0.80,-0.5) in node \[left\] at (-1.58,-1.9); on (-0.80,0.40) in node \[left\] at (2.25,1.92);

(-85,9)[ ]{}

Here, we simulated low-dose CT measurements from full-dose CT images in Mayo Clinic dataset with GE 2D LightSpeed fan-beam geometry corresponding to a monoenergetic source with $10^4$ incident photons per ray and no scatter. We ran 100 iterations of the PWLS-EP algorithm with FBP reconstructions as initializations. We used the relaxed OS-LALM algorithm with $4$ ordered subsets and regularization parameter $\beta = 2^{15.5}$. For the MRST model, we used the OS-LALM algorithm for the image update step with $2$ inner iterations and $4$ subsets. We used $T_O = 1000$ and $T_O = 1500$ outer iterations for (ST, MRST2) and (MRST3, MRST5, MRST7), respectively. We set the regularization parameters (after tuning over ranges of values) as $(\beta$, $\gamma)$ $=$ $(7\times10^4$, $20)$ for ST, $(\beta$, $\gamma_1$, $\gamma_2)$ $=$ $(3\times10^4$, $30$, $12)$ for MRST2, $(\beta$, $\gamma_1$, $\gamma_2$, $\gamma_3)$ $=$ $(3\times10^4$, $30$, $12$, $10)$ for MRST3, $(\beta$, $\gamma_1$, $\gamma_2$, $\gamma_3$, $\gamma_4$, $\gamma_5)$ $=$ $(4\times10^4$, $30$, $20$, $10$, $10$, $10)$ for MRST5, and $(\beta$, $\gamma_1$, $\gamma_2$, $\gamma_3$, $\gamma_4$, $\gamma_5$, $\gamma_6$, $\gamma_7)$ $=$ $(7\times10^4$, $30$, $12$, $10$, $10$, $10$, $10$, $10)$ for MRST7, respectively.

To compare the performance between various models quantitatively, we compute the root mean square error (RMSE), peak signal to noise ratio (PSNR), structural similarity index measurement (SSIM) of the reconstructions in a region of interest (ROI). The RMSE in Hounsfield units (HU) computed between the ground truth image and reconstruction image is defined as RMSE $= \sqrt{\Sigma_{i \in \text{ROI}}(\hat{x}_i-x^*_i)^2/{N_{\text{ROI}}}}$, where $\hat{x}_i$ and $x^*_i$ denote the pixel intensities of the reconstructed and ground truth images, respectively, and $N_{\text{ROI}}$ is the number of pixels in the ROI. The ROI here was a circular (around center) region containing all the phantom tissues. Since the full-dose Mayo Clinic CT images still contain some considerable noise, we cannot solely rely on RMSE and PSNR as metrics of image quality, as they would include the noise. SSIM helps better evaluate the preservation of structural details in the reconstructed image.

We conduct experiments on three slices (L067-slice100, L192-slice100, L506-slice100) of the Mayo Clinic data. Fig. \[fig:recon\_mayo\] shows the reconstruction of L067-slice100 using FBP, PWLS-EP, PWLS-ST, PWLS-MRST2, PWLS-MRST3, PWLS-MRST5, and PWLS-MRST7, respectively at incident photon intensity $I_0=10^{4}$. TABLE \[tab:MRST\] lists the RMSE, PSNR, and SSIM values of reconstructions of the three test slices, with the best values bolded. The two-layer model (MRST2) provides the best RMSE and PSNR values among the methods. However, when we consider the SSIM criterion, MRST5 and MRST7 outperform ST and MRST2. So which MRST model is better? By observing the reconstructed images, we see that although MRST2 and ST have lower RMSE and higher PSNR values than MRST5 and MRST7, they sacrifice some sharpness of the central region and suffer from loss of details. The deeper models have a more positive effect in maintaining subtle features, which is clearly more essential to clinic medical diagnosis. Furthermore, after considerable parameter tuning, we have observed that the deeper models offer more stable image quality as $\beta$ is varied, i.e., they are more robust to oversmoothing.

Conclusion
==========

This paper proposes a general framework for Multi-layer Residual Sparsifying Transform (MRST) learning wherein the transform domain residual maps over several layers are jointly sparsified. Our work then applies MRST learning to low-dose CT (LDCT) image reconstruction by using a PWLS approach with a learned MRST-based regularizer. Experimental results illustrate the promising performance of the multi-layer scheme over single-layer learned sparsifying transforms. Learned MRST models also offer significant improvement over typical nonadaptive methods. In the future, we will investigate additional strategies such as pooling operations to reduce noise over layers.

Acknowledgment
==============

The authors thank Dr. Cynthia McCollough, the Mayo Clinic, the American Association of Physicists in Medicine, and the National Institute of Biomedical Imaging and Bioengineering for providing the Mayo Clinic data.

[^1]: X. Yang, X. Zheng and Y. Long are with the University of Michigan– Shanghai Jiao Tong University Joint Institute, Shanghai Jiao Tong University, Shanghai 200240, China (e-mail: yangxk@sjtu.edu.cn; zhxhang@sjtu.edu.cn; yong.long@sjtu.edu.cn)

[^2]: S. Ravishankar is with the Department of Computational Mathematics, Science and Engineering and Department of Biomedical Engineering, Michigan State University, East Lansing, MI 48824, USA (e-mail: ravisha3@msu.edu)

[^3]: $^\star$Yong Long is the corresponding author.
