PWLS-ULTRA: An Efficient Clustering and Learning-Based Approach for Low-Dose 3D CT Image Reconstruction

TL;DR

This paper introduces PWLS-ULTRA, a novel low-dose CT reconstruction method that leverages learned transform unions and efficient optimization to produce higher quality images faster than existing techniques.

Contribution

The paper presents a new PWLS-based reconstruction method using a union of learned transforms, improving image quality and computational efficiency in low-dose CT imaging.

Findings

Significantly improved image quality over traditional PWLS-EP.

Better reconstruction quality than single transform methods.

Faster computation compared to learned dictionaries.

Abstract

The development of computed tomography (CT) image reconstruction methods that significantly reduce patient radiation exposure while maintaining high image quality is an important area of research in low-dose CT (LDCT) imaging. We propose a new penalized weighted least squares (PWLS) reconstruction method that exploits regularization based on an efficient Union of Learned TRAnsforms (PWLS-ULTRA). The union of square transforms is pre-learned from numerous image patches extracted from a dataset of CT images or volumes. The proposed PWLS-based cost function is optimized by alternating between a CT image reconstruction step, and a sparse coding and clustering step. The CT image reconstruction step is accelerated by a relaxed linearized augmented Lagrangian method with ordered-subsets that reduces the number of forward and back projections. Simulations with 2-D and 3-D axial CT scans of the extended cardiac-torso phantom and 3D helical chest and abdomen scans show that for both normal-dose and low-dose levels, the proposed method significantly improves the quality of reconstructed images compared to PWLS reconstruction with a nonadaptive edge-preserving regularizer (PWLS-EP). PWLS with regularization based on a union of learned transforms leads to better image reconstructions than using a single learned square transform. We also incorporate patch-based weights in PWLS-ULTRA that enhance image quality and help improve image resolution uniformity. The proposed approach achieves comparable or better image quality compared to learned overcomplete synthesis dictionaries, but importantly, is much faster (computationally more efficient).

Tables

Table 1. TABLE I: RMSE (HU) and SSIM of 2D (fan-beam) image reconstructions with FBP, PWLS-EP, PWLS-DL, PWLS-ST, PWLS-ULTRA ( K = 15 𝐾 15 K=15 ), and PWLS-ULTRA ( K = 15 𝐾 15 K=15 ) with patch-based weights ( τ j subscript 𝜏 𝑗 \tau_{j} ), for two incident photon intensities

Intensity	FBP	EP	DL	ST	ULTRA	ULTRA-$\left\{{\tau }_j\right\}$
$1\times {10}^4$	73.7	39.4	33.6	36.5	34.4	33.1
	0.547	0.892	0.966	0.966	0.967	0.969
$5\times {10}^3$	89.0	49.7	39.1	43.9	39.8	38.9
	0.472	0.884	0.958	0.955	0.953	0.956

Table 2. TABLE II: RMSE (HU) in three ROIs of 2D (fan-beam) image reconstructions with FBP, PWLS-EP, PWLS-DL, PWLS-ST, PWLS-ULTRA ( K = 15 𝐾 15 K=15 ), and PWLS-ULTRA ( K = 15 𝐾 15 K=15 ) with patch-based weights ( τ j subscript 𝜏 𝑗 \tau_{j} ), for two incident photon intensities

Intensity	Methods	ROI-$1$	ROI-$2$	ROI-$3$
$1\times {10}^4$	FBP	21.8	15.6	39.6
	EP	6.6	10.9	14.7
	DL	3.7	9.9	16.6
	ST	3.9	10.8	14.1
	ULTRA	4.2	9.6	13.8
	ULTRA-$\left\{{\tau }_j\right\}$	4.2	9.3	12.6
$5\times {10}^3$	FBP	51.7	36.4	39.0
	EP	7.1	14.9	28.5
	DL	7.0	14.5	20.7
	ST	6.3	14.3	21.4
	ULTRA	5.8	13.7	17.5
	ULTRA-$\left\{{\tau }_j\right\}$	5.9	13.7	18.1

Table 3. TABLE III: RMSE (HU) and SSIM of 3D (cone-beam) reconstructions with FDK, PWLS-EP, PWLS-ST, PWLS-ULTRA ( K = 15 𝐾 15 K=15 ), and PWLS-ULTRA ( K = 15 𝐾 15 K=15 ) with patch-based weights ( τ j subscript 𝜏 𝑗 \tau_{j} ), for two incident photon intensities

Intensity	FDK	EP	ST	ULTRA	ULTRA-$\left\{{\tau }_j\right\}$
$1\times {10}^4$	67.8	34.6	32.1	30.7	29.2
	0.536	0.940	0.976	0.978	0.981
$5\times {10}^3$	89.0	41.1	37.3	35.7	34.2
	0.463	0.921	0.967	0.970	0.974

Equations

$$\underset{\mathbf{\Omega },\mathbf{Z}}{\min }\Vert \mathbf{\Omega X}-\mathbf{Z}{\Vert }_F^2+\lambda Q(\mathbf{\Omega })+\sum _{i=1}^{N^{\prime }}{\eta }^2\Vert {\mathbf{Z}}_i{\Vert }_0$$

$$\underset{\mathbf{x}⪰0}{\min }\frac{1}{2}\Vert \mathbf{y}-\mathbf{Ax}{\Vert }_{\mathbf{W}}^2+\beta \mathsf{R}(\mathbf{x})$$

$$\mathsf{R}(\mathbf{x})\triangleq\min_{\{\mathbf{z}_{j}\}}\sum_{j=1}^{\tilde{N}}\tau_{j}\bigg{\{}\|\mathbf{\Omega}\mathbf{P}_{j}\mathbf{x}-\mathbf{z}_{j}\|^{2}_{2}+\gamma^{2}\|\mathbf{z}_{j}\|_{0}\bigg{\}}$$

$${\tau }_j\triangleq \Vert {\mathbf{P}}_j\mathbf{\kappa }{\Vert }_1/l$$

$$\underset{\{{\mathbf{\Omega }}_k,{\mathbf{Z}}_i,C_k\}}{\min }$$

$$+\sum _{k=1}^K{\lambda }_{k}Q({\mathbf{\Omega }}_k)s.t.\{C_k\}\in \mathcal{G}.$$

$$\mathsf{R}(\mathbf{x})\triangleq \underset{\{{\mathbf{z}}_j,C_k\}}{\min }$$

$$s.t.\{C_k\}\in \mathcal{G}.$$

$$\underset{\{{\mathbf{\Omega }}_k\}}{\min }\sum _{k=1}^K\sum _{i\in C_k}\Vert {\mathbf{\Omega }}_k{\mathbf{X}}_i-{\mathbf{Z}}_i{\Vert }_2^2+\sum _{k=1}^K{\lambda }_{k}Q({\mathbf{\Omega }}_k).$$

$$\underset{{\mathbf{\Omega }}_k}{\min }\sum _{i\in C_k}\Vert {\mathbf{\Omega }}_k{\mathbf{X}}_i-{\mathbf{Z}}_i{\Vert }_2^2+{\lambda }_{k}Q({\mathbf{\Omega }}_k).$$

$$\hat{\mathbf{\Omega}}_{k}=0.5\mathbf{R}\big{(}{\mathbf{\Sigma}}+({\mathbf{\Sigma}}^{2}+2\lambda_{k}\mathbf{I})^{\frac{1}{2}}\big{)}\mathbf{Q}^{T}\mathbf{L}^{-1}.$$

$$\min_{\{\mathbf{Z}_{i},C_{k}\}}\sum_{k=1}^{K}\sum_{i\in C_{k}}\bigg{\{}\|\mathbf{\Omega}_{k}\mathbf{X}_{i}-\mathbf{Z}_{i}\|_{2}^{2}+\eta^{2}\|\mathbf{Z}_{i}\|_{0}+\lambda_{0}\|\mathbf{X}_{i}\|_{2}^{2}Q(\mathbf{\Omega}_{k})\bigg{\}}$$

$$\underset{\mathbf{x}⪰0}{\min }\frac{1}{2}\Vert \mathbf{y}-\mathbf{Ax}{\Vert }_{\mathbf{W}}^2+{\mathsf{R}}_2(\mathbf{x})$$

$${\mathbf{D}}_{\mathbf{A}}\triangleq \mathsf{diag}\{{\mathbf{A}}^{T}WA1\}⪰{\mathbf{A}}^T\mathbf{WA}.$$

$$\rho_{r}(\alpha)=\begin{cases}1,&r=0\\ \frac{\pi}{\alpha(r+1)}\sqrt{1-\big{(}\frac{\pi}{2\alpha(r+1)}\big{)}^{2}},&\text{otherwise,}\end{cases}$$

$${\mathbf{D}}_{\mathsf{R}}$$

$$⪰2\beta \sum _{k=1}^K\sum _{j\in C_k}{\tau }_j{\mathbf{P}}_j^T{\mathbf{\Omega }}_k^T{\mathbf{\Omega }}_k{\mathbf{P}}_j={\nabla }^2{\mathsf{R}}_2(\mathbf{x}).$$

$$\min_{\{\mathbf{z}_{j}\},\{C_{k}\}\in\mathcal{G}}\sum_{k=1}^{K}\bigg{\{}\sum_{j\in C_{k}}\tau_{j}\big{\{}\|\mathbf{\Omega}_{k}\mathbf{P}_{j}\mathbf{x}-\mathbf{z}_{j}\|^{2}_{2}+\gamma^{2}\|\mathbf{z}_{j}\|_{0}\big{\}}\bigg{\}}.$$

$${\hat{k}}_j=\underset{1\le k\le K}{\mathrm{argmin}}\Vert {\mathbf{\Omega }}_k{\mathbf{P}}_j\mathbf{x}-H_{\gamma }({\mathbf{\Omega }}_k{\mathbf{P}}_j\mathbf{x}){\Vert }_2^2+{\gamma }^2\Vert H_{\gamma }({\mathbf{\Omega }}_k{\mathbf{P}}_j\mathbf{x}){\Vert }_0$$

$$\{\begin{array}{rl}{\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{R}}{)}^{-1}({\mathbf{s}}^{(r+1)}+\nabla {\mathsf{R}}_2({\mathbf{x}}^{(r)})){]}_{+}\\ {\mathbf{\zeta }}^{(r+1)}& \triangleq M{\mathbf{A}}_m^T{\mathbf{W}}_m({\mathbf{A}}_m{\mathbf{x}}^{(r+1)}-{\mathbf{y}}_m)\\ {\mathbf{g}}^{(r+1)}& =\frac{\rho }{\rho +1}(\alpha {\mathbf{\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)}-{\mathbf{\zeta }}^{(r+1)})+(1-\alpha ){\mathbf{h}}^{(r)}\end{array}$$

$$\mathbf{D}_{\mathsf{R}}\triangleq 2\beta\bigg{\{}\max_{k}\lambda_{\max}(\mathbf{\Omega}_{k}^{T}\mathbf{\Omega}_{k})\bigg{\}}\sum_{k=1}^{K}\sum_{j\in C_{k}}\tau_{j}\mathbf{P}_{j}^{T}\mathbf{P}_{j}.$$

Full text

I Introduction

There is a growing interest in techniques for computed tomography (CT) image reconstruction that significantly reduce patient radiation exposure while maintaining high image quality. Dictionary learning based techniques have been proposed for low-dose CT (LDCT) imaging, but often involve expensive computation. This paper proposes a new penalized weighted least aquares (PWLS) reconstruction approach that exploits regularization based on an efficient Union of Learned TRAnsforms (PWLS-ULTRA). In the following, we briefly review recent methods for LDCT image reconstruction and summarize the contributions of this work.

I-A Background

Various methods have been proposed for image reconstruction in LDCT imaging. When radiation dose is reduced, analytical filtered back-projection (FBP) image reconstruction methods (e.g., the Feldkamp-Davis-Kress or FDK method [1]) typically provide unacceptable image quality. For example, streak artifacts increase severely as radiation dose is reduced [2]. Model-based image reconstruction (MBIR) methods, aka statistical image reconstruction (SIR) methods, can provide high-quality reconstructions from low-dose scans [3, 4]. These methods iteratively find the image based on the system (physical) model, the measurement statistical model, and (assumed) prior information about the unknown object. A typical MBIR method for CT uses a penalized weighted-least squares (PWLS) cost function with a statistically weighted quadratic data-fidelity term and a penalty term (regularizer) modeling prior knowledge of the underlying unknown object [5, 6, 7].

Many current LDCT reconstruction methods use simple prior information. Adopting better image priors in MBIR could substantially improve image reconstruction quality for LDCT scans. The prior image constrained compressed sensing (PICCS) method was first proposed to enable accurate reconstruction of CT images from highly undersampled projection data sets [8, 9, 10]. Since a normal-dose CT image scanned previously may be available in some clinical applications, dose reduction using prior image constrained compressed sensing (DR-PICCS) was proposed to reduce image noise [11]. Ma et al. [12] proposed the previous normal-dose scan induced nonlocal means (ndiNLM) method to utilize the normal-dose image to enable low dose CT image reconstruction. The ndiNLM method expects that the normal-dose and the current low-dose scans are spatially aligned, and determines optimal local weights from the normal-dose image to improve the NLM weighted average [12, 13]. The PICCS and ndiNLM class of methods incorporate prior information from corresponding normal-dose CT images, assumed available. We propose a method that differs from these approaches in that it does not require prior normal-dose images of the same patient or object, and can rather learn general CT image features or filters from diverse image sets and datasets.

Extracting prior information from big datasets of CT images has great potential to enable MBIR methods to produce significantly improved reconstructions from LDCT measurements. Images are often sparse in certain transform domains (such as wavelets, discrete cosine transform, and discrete gradient) or dictionaries. The synthesis dictionary model approximates a signal by a linear combination of a few columns or atoms of a pre-specified dictionary [14]. The choice of the synthesis dictionary is critical for the success of sparse representation modeling and other applications [15]. The data-driven adaptation of dictionaries, or dictionary learning [16, 17, 18, 19, 20] yields dictionaries with better sparsifying capability for specific classes of signals than analytic dictionaries based on mathematical models. Such learned dictionaries have been widely exploited in various applications in recent years, including super-resolution imaging, image or video denoising, classification, and medical image reconstruction [21, 22, 23, 24, 25, 26, 27]. Some recent works also studied parametrized models such as adaptive tight frames [28], multivariate Gaussian mixture distributions [29], and shape dictionaries [30].

Recently, Xu et al. [31] applied dictionary learning to 2D LDCT image reconstruction by proposing a PWLS approach with an overcomplete synthesis dictionary-based regularizer. Their method uses either a global dictionary trained from 2D image patches extracted from a normal-dose FBP image, or an adaptive dictionary jointly estimated with the low-dose image. The trained global dictionary worked better than the adaptively estimated dictionary for highly limited (e.g., with very few views, or ultra-low dose) data. Several works proposed 3D CT reconstruction by learning either a 3D dictionary from 3D image patches, or learning three 2D dictionaries (dubbed 2.5D) from image patches extracted from slices along the x-y, y-z, and x-z directions, respectively [32, 33].

Dictionary learning methods typically alternate between estimating the sparse coefficients of training signals or image patches (sparse coding step) and updating the dictionary (dictionary update step). The sparse coding step in both synthesis dictionary learning [21, 18] and analysis dictionary learning [34] is NP-Hard (Non-deterministic Polynomial-time hard) in general, and algorithms such as K-SVD [21, 18] involve relatively expensive computations for sparse coding. A recent generalized analysis dictionary learning approach called sparsifying transform learning [35, 36] more efficiently learns a transform model for signals. The transform model assumes that a signal $\mathbf{x}\in\mathbb{R}^{n}$ is approximately sparsifiable using a transform $\mathbf{\Omega}\in\mathbb{R}^{m\times n}$, i.e., $\mathbf{\Omega}\mathbf{x}=\mathbf{z}+\mathbf{e}$ where $\mathbf{z}\in\mathbb{R}^{m}$ is sparse in some sense, and $\mathbf{e}\in\mathbb{R}^{m}$ denotes the modeling error in the transform domain. Transform learning methods typically alternate between sparse approximation of training signals in the transform domain (sparse coding step) and updating the transform operator (transform update step). In contrast to dictionary learning methods, the sparse coding step in transform learning involves simple thresholding [35, 36]. Transform learning methods have been recently demonstrated to work well in applications [37, 38, 39, 40]. Pfister and Bresler [41, 42, 43] showed the promise of PWLS reconstruction with adaptive square transform-based regularization, wherein they jointly estimated the square transform (ST) and the image. Pre-training a (global) transform from a large dataset would save computations during CT image reconstruction, and may also be well-suited for highly limited data (evidenced earlier for dictionary learning in [31]).

Wen *et al. *recently extended the single ST learning method to learning a union of square transforms model, also referred to as an overcomplete transform with block cosparsity (OCTOBOS) [44]. This transform learning approach jointly adapts a collection (or union) of $K$ square transforms and clusters the signals or image patches into $K$ groups. Each (learned) group of signals is well-matched to a corresponding transform in the collection. Such a learned union of transforms outperforms the ST model in applications such as image denoising [44].

I-B Contributions

Incorporating the efficient square transform (ST) model, we propose a new PWLS approach for LDCT reconstruction that exploits regularization based on a pre-learned square transform (PWLS-ST). We also extend this approach to a more general PWLS scheme involving a Union of Learned TRAnsforms (PWLS-ULTRA). The transform models are pre-learned from numerous patches extracted from a dataset of CT images or volumes. We also incorporate patch-based weights in the proposed regularizer to help improve image resolution or noise uniformity. We propose an efficient iterative algorithm for the PWLS costs that alternates between a sparse coding and clustering step (which reduces to a sparse coding step for PWLS-ST) that uses closed-form solutions, and an iterative image update step. There are several iterative algorithms that could be used for the image update step such as the preconditioned conjugate gradient (PCG) method [45], the separable quadratic surrogate method with ordered-subsets based acceleration (OS-SQS) [46], iterative coordinate descent (ICD) [47], splitting-based algorithms [48], and the optimal gradient method (OGM) [49]. We chose the relaxed linearized augmented Lagrangian method with ordered-subsets (relaxed OS-LALM) [50] for the image update step.

The proposed PWLS-ULTRA approach clusters the voxels into different groups. These groups often capture features such as bones, specific soft tissues, edges, etc. Experiments with 2D and 3D axial CT scans of the extended cardiac-torso (XCAT) phantom and 3D helical chest and abdomen scans show that for both normal-dose and low-dose levels, the proposed methods significantly improve the quality of reconstructed images compared to conventional reconstruction methods such as filtered back-projection or PWLS reconstruction with a nonadaptive edge-preserving regularizer (PWLS-EP). The union of learned transforms provides better image reconstruction quality than using a single learned square transform. The proposed PWLS-ULTRA achieves comparable or better image quality compared to learned overcomplete synthesis dictionaries, but importantly, is much faster (computationally more efficient).

We presented a brief study of PWLS-ST for low-dose fan-beam (2D) CT image reconstruction in [51]. This paper investigates the more general PWLS-ULTRA framework, and presents experimental results illustrating the properties of the PWLS-ST and PWLS-ULTRA algorithms and demonstrating their performance for low-dose fan-beam, cone-beam (3D) and helical (3D) CT.

I-C Organization

Section II describes the formulations for pre-learning a square transform or a union of transforms, and the formulations for PWLS reconstruction with regularization based on learned sparsifying transforms. Section III derives efficient optimization algorithms for the proposed problems. Section IV presents experimental results illustrating properties of the proposed algorithms and demonstrating their promising performance for LDCT reconstruction compared to numerous recent methods. Section V presents our conclusions and mentions areas of future work.

II Problem Formulations for Transform Learning and Image Reconstruction

II-A PWLS-ST Formulation for LDCT Reconstruction

Given $N^{\prime}$ vectorized image patches (2D or 3D) extracted from a dataset of CT images or volumes, we learn a square transform $\mathbf{\Omega}\in\mathbb{R}^{l\times l}$ by solving the following (training) optimization problem:

[TABLE]

where $l$ is the number of pixels in each patch, $\lambda=\lambda_{0}\|\mathbf{X}\|_{F}^{2}$ ($\lambda_{0}>0$ is a constant) and $\eta>0$ are scalar parameters, and $\{\mathbf{Z}_{i}\}_{i=1}^{N^{\prime}}$ denote the sparse codes of the training signals (vectorized patches) $\{\mathbf{X}_{i}\}_{i=1}^{N^{\prime}}$. Matrices $\mathbf{X}\in\mathbb{R}^{l\times N^{\prime}}$ and $\mathbf{Z}\in\mathbb{R}^{l\times N^{\prime}}$ have the training signals and sparse codes respectively, as their columns. The $\ell_{0}$ “norm” counts the number of non-zeros in a vector. The term $\|\mathbf{\Omega}\mathbf{X}-\mathbf{Z}\|_{F}^{2}$ is called the sparsification error and measures the deviation of the signals in the transform domain from their sparse approximations. Regularizer $Q(\mathbf{\Omega})\triangleq\|\mathbf{\Omega}\|_{F}^{2}-\log|\det\mathbf{\Omega}|$ prevents trivial solutions and controls the condition number of $\mathbf{\Omega}$ [36].

After a transform $\mathbf{\Omega}$ is learned, we reconstruct an (vectorized) 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 optimization problem [51]:

[TABLE]

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

[TABLE]

where $\tilde{N}$ is the number of image patches, the operator $\mathbf{P}_{j}\in\mathbb{R}^{l\times N_{p}}$ extracts the $j$th patch of $l$ voxels of $\mathbf{x}$ as $\mathbf{P}_{j}\mathbf{x}$, and vector $\mathbf{z}_{j}\in\mathbb{R}^{l}$ denotes the transform-sparse representation of $\mathbf{P}_{j}\mathbf{x}$. The regularizer includes a sparsification error term and a $\ell_{0}$ “norm”-based sparsity penalty with weight $\gamma^{2}$ ($\gamma>0$).

We also include patch-based weights $\{\tau_{j}\}$ in (1) to encourage uniform spatial resolution or uniform noise in the reconstructed image [52] as follows:

[TABLE]

with $\bm{\kappa}$ (of same size as $\mathbf{x}$) whose elements $\kappa_{j}$ are defined in terms of the entries of $\mathbf{A}$ (denoted $a_{ij}$) and $\mathbf{W}$ as $\kappa_{j}\triangleq\sqrt{{\sum_{i=1}^{N_{d}}a_{ij}w_{i}}/{\sum_{i=1}^{N_{d}}a_{ij}}}$ [53, eq(39)]. While (2) uses the $\ell_{1}$ norm, corresponding to the mean value of $\mathbf{P}_{j}\bm{\kappa}$, to define $\tau_{j}$, we have observed that other alternative norms also work well in practice for LDCT reconstruction.

Fig. 1 shows example transforms (rows of $\mathbf{\Omega}$ are reshaped as $8\times 8\times 8$ patches and the first $8\times 8$ slices of $256$ such 3D patches are shown) learned from $8\times 8\times 8$ patches of an XCAT phantom [54] volume. The transform learned with $\eta=100$ in (P0) has more oriented features whereas the transform learned with $\eta=50$ shows more gradient (or finite-difference) type features (pointed by the green arrows). This behavior suggests that a single ST may not be rich enough to capture the diverse features, edges, and other properties of CT volumes. Therefore, next we consider the extension of the ST approach to a rich union of learned transforms scheme.

II-B Learning a Union of Sparsifying Transforms

To learn a union of sparsifying transforms $\left\{\mathbf{\Omega}_{k}\right\}_{k=1}^{K}$ from $N^{\prime}$ (vectorized) patches, we solve

[TABLE]

This formulation groups the training signals $\{\mathbf{X}_{i}\}$ into $K$ classes according to the transform they best match, and $C_{k}$ denotes the set of indices of signals matched to the $k$th class. Set $\mathcal{G}$ denotes all possible partitionings of $\left\{1,2,..,N^{\prime}\right\}$ into $K$ disjoint subsets. We use $K$ regularizers $Q(\mathbf{\Omega}_{k})\triangleq\|\mathbf{\Omega}_{k}\|_{F}^{2}-\log|\det\mathbf{\Omega}_{k}|$, $1\leq k\leq K$, to control the properties of the transforms. We set these regularizer weights as $\lambda_{k}=\lambda_{0}\|\mathbf{X}_{C_{k}}\|_{F}^{2}$ [44], where $\lambda_{0}>0$ is a constant and $\mathbf{X}_{C_{k}}$ is a matrix whose columns are the training signals in the $k$th cluster. This choice of $\{\lambda_{k}\}$ together with $\eta=\eta_{0}\|\mathbf{X}\|_{F}$ for $\eta_{0}>0$ allows the terms in (P2) to scale appropriately with the data. Problem (P2) learns a collection of transforms and a clustering for the image patches, together with the patches’ sparse coefficients $\{\mathbf{Z}_{i}\}$. The next section uses these transforms for image reconstruction.

II-C LDCT Reconstruction with ULTRA Regularization

We propose a PWLS-ULTRA framework, where we solve (P1) but with the regularizer $\mathsf{R}(\mathbf{x})$ defined based on a union of sparsifying transforms as

[TABLE]

This regularizer measures the sparsification error of each patch using its best-matched transform. Using (3), (P1) estimates the image $\mathbf{x}$, the sparse coefficients of image patches $\{\mathbf{z}_{j}\}$, and the cluster assignments $\{C_{k}\}$ from LDCT sinogram data $\mathbf{y}$.

III Algorithms and Properties

The square transform learning and the PWLS-ST formulations are special cases (corresponding to $K=1$) of the ULTRA-based formulations. Therefore, this section describes algorithms for solving (P1) with regularizer (3) and (P2).

III-A Algorithm for Training a Union of Transforms

We adopt an alternating minimization algorithm for (P2) that alternates between a transform update step (solving for $\{\mathbf{\Omega}_{k}\}$) and a sparse coding and clustering step (solving for $\{\mathbf{Z}_{i},C_{k}\}$). These steps are described next.

III-A1 Transform Update Step

With $\{\mathbf{Z}_{i},C_{k}\}$ fixed, we solve the following optimization problem for $\{\mathbf{\Omega}_{k}\}$ [44]:

[TABLE]

Since the objective is in summation form, the above problem separates into $K$ independent single transform learning problems that we solve in parallel. The $k$th such optimization problem is as follows:

[TABLE]

We update the transform $\mathbf{\Omega}_{k}$ following prior work [44, 36]. Let $\mathbf{Q}\mathbf{\Sigma}\mathbf{R}^{T}$ denote the full singular value decomposition of $\mathbf{L}^{-1}\mathbf{X}_{C_{k}}\mathbf{Z}_{C_{k}}^{T}$, with $\mathbf{L}\mathbf{L}^{T}\triangleq\mathbf{X}_{C_{k}}\mathbf{X}_{C_{k}}^{T}+\lambda_{k}\mathbf{I}$ (i.e., $\mathbf{L}$ is a matrix square root). Then, the minimizer of (5) is

[TABLE]

III-A2 Sparse Coding and Clustering Step

With $\{\mathbf{\Omega}_{k}\}$ fixed, we solve the following sub-problem for $\{\mathbf{Z}_{i},C_{k}\}$:

[TABLE]

For given cluster memberships, the optimal sparse codes are ${\mathbf{Z}}_{i}=H_{\eta}(\mathbf{\Omega}_{k}\mathbf{X}_{i}),\forall i\in C_{k}$, where the hard-thresholding operator $H_{\eta}(\cdot)$ zeros out vector entries with magnitude less than $\eta$. Using this result, it follows that the optimal cluster membership for each $\mathbf{X}_{i}$ in (7) is \hat{k}_{i}=\underset{1\leq k\leq K}{\operatorname{arg\,min}}\big{\{}\|\mathbf{\Omega}_{k}\mathbf{X}_{i}-H_{\eta}(\mathbf{\Omega}_{k}\mathbf{X}_{i})\|^{2}_{2}+\eta^{2}\|H_{\eta}(\mathbf{\Omega}_{k}\mathbf{X}_{i})\|_{0}+\lambda_{0}\|\mathbf{X}_{i}\|_{2}^{2}Q(\mathbf{\Omega}_{k})\big{\}}, and the corresponding optimal sparse code is $\hat{\mathbf{Z}}_{i}=H_{\eta}(\mathbf{\Omega}_{\hat{k}_{i}}\mathbf{X}_{i})$.

III-B PWLS-ULTRA Image Reconstruction Algorithm

We propose an alternating algorithm for the PWLS-ULTRA formulation (i.e., (P1) with regularizer (3)) that alternates between updating $\mathbf{x}$ (image update step), and $\{\mathbf{z}_{j},C_{k}\}$ (sparse coding and clustering step).

III-B1 Image Update Step

With $\{\mathbf{z}_{j},C_{k}\}$ fixed, (P1) for PWLS-ULTRA reduces to the following weighted least squares problem:

[TABLE]

where $\mathsf{R}_{2}(\mathbf{x})\triangleq\beta\sum_{k=1}^{K}\sum_{j\in C_{k}}\tau_{j}\|\mathbf{\Omega}_{k}\mathbf{P}_{j}\mathbf{x}-\mathbf{z}_{j}\|^{2}_{2}$.

We solve (8) using the recent relaxed OS-LALM [50], whose iterations are shown in Algorithm 1. Here, for each iteration $n$, we further iterate over $1\leq m\leq M$ corresponding to $M$ ordered subsets. The matrices $\mathbf{A}_{m}$, $\mathbf{W}_{m}$, and the vector $\mathbf{y}_{m}$ in Algorithm 1 are sub-matrices of $\mathbf{A}$, $\mathbf{W}$, and sub-vector of $\mathbf{y}$, respectively, for the $m$th subset. Matrix $\mathbf{D}_{\mathbf{A}}\succeq\mathbf{A}^{T}\mathbf{W}\mathbf{A}$ is a diagonal majorizing matrix of $\mathbf{A}^{T}\mathbf{W}\mathbf{A}$; specifically we use [46]

[TABLE]

The gradient $\nabla\mathsf{R}_{2}(\mathbf{x})=2\beta\sum_{k=1}^{K}\sum_{j\in C_{k}}\tau_{j}\mathbf{P}_{j}^{T}\mathbf{\Omega}_{k}^{T}(\mathbf{\Omega}_{k}\mathbf{P}_{j}\mathbf{x}-\mathbf{z}_{j})$, the (over-)relaxation parameter $\alpha\in[1,2)$, and the parameter $\rho>0$ decreases gradually with iteration [50],

[TABLE]

where $r$ indexes the total number of $n$ and $m$ iterations. Lastly, $\mathbf{D}_{\mathsf{R}}$ in Algorithm 1 is a diagonal majorizing matrix of the Hessian of the regularizer $\mathsf{R}_{2}(\mathbf{x})$, specifically:

[TABLE]

Since this $\mathbf{D}_{\mathsf{R}}$ is independent of $\mathbf{x}$, $\{\mathbf{z}_{j}\}$, and $\{C_{k}\}$, we precompute it using patch-based operations [25] (cf. the supplement111Supplementary material is available in the supplementary files/multimedia tab. for details) prior to iterating.

III-B2 Sparse Coding and Clustering Step

With $\mathbf{x}$ fixed, we solve the following sub-problem to determine the optimal sparse codes and cluster assignments for each patch:

[TABLE]

For each patch $\mathbf{P}_{j}\mathbf{x}$, with (optimized) ${\mathbf{z}}_{j}=H_{\gamma}(\mathbf{\Omega}_{k}\mathbf{P}_{j}\mathbf{x})$, the optimal cluster assignment is computed as follows:

[TABLE]

Minimizing over $k$ above finds the best-matched transform. Then, the optimal sparse codes are $\hat{\mathbf{z}}_{j}=H_{\gamma}(\mathbf{\Omega}_{\hat{k}_{j}}\mathbf{P}_{j}\mathbf{x})$.

III-B3 Overall Algorithm

The proposed method for the PWLS-ULTRA problem is shown in Algorithm 1. The algorithm for the PWLS-ST formulation is obtained by setting $K=1$ and skipping the clustering procedure in the sparse coding and clustering step. Algorithm 1 uses an initial image estimate and the union of pre-learned transforms $\{\mathbf{\Omega}_{k}\}$. It then alternates between the image update, and sparse coding and clustering steps until a convergence criterion (such as $\|\tilde{\mathbf{x}}^{(t+1)}-\tilde{\mathbf{x}}^{(t)}\|_{2}<\epsilon$ for some small $\epsilon>0$) is satisfied, or alternatively until some maximum inner/outer iteration counts are reached.

III-B4 Computational Cost

Each outer iteration of the proposed Algorithm 1 involves the image update and the sparse coding and clustering steps. The cost of the sparse coding and clustering step scales as $O(l^{2}N)$ and is dominated by matrix-vector products. Importantly, unlike prior dictionary learning-based works [31], where the computations for the sparse coding step (involving orthogonal matching pursuit (OMP) [55]) can scale worse as $O(l^{3}N)$ (assuming synthesis sparsity levels of patches $\propto l$), the exact sparse coding and clustering in PWLS-ULTRA is cheaper, especially for large patch sizes. Similar to prior works [31], the computations in the image update step are dominated by the forward and back projection operations. Section IV compares the proposed method to synthesis dictionary learning-based approaches, and shows that our transform approach runs much faster.

IV Experimental Results

This section presents experimental results illustrating properties of the proposed algorithms and demonstrating their promising performance for LDCT reconstruction compared to numerous recent methods. We include additional experimental results in the supplement. A link to software to reproduce our results is provided at http://web.eecs.umich.edu/~fessler/irt/reproduce/.

IV-A Framework and Data

We evaluate the proposed PWLS-ULTRA and PWLS-ST (i.e., with $K=1$) methods for 2D fan-beam and 3D axial cone-beam CT reconstruction of the XCAT phantom [54]. We also apply the proposed methods to helical CT clinical data of the chest and abdomen.

Section IV-B discusses the role and intuition of each parameter in the proposed methods. Section IV-C illustrates the properties of the transform learning and image reconstruction methods. Sections IV-D and IV-E show results for 2D fan-beam and 3D axial cone-beam CT, respectively, for the XCAT phantom data. We used the “Poisson + Gaussian” model, i.e., $\tilde{k}\,\text{Poisson}\{I_{0}\exp(-[\mathbf{A}\mathbf{x}]_{i})\}+\text{Normal}\{0,\sigma^{2}\}$ to simulate CT measurements of the XCAT phantom, where $I_{0}$ is the incident X-ray intensity incorporating X-ray source illumination and the detector gain, the parameter $\tilde{k}=1$ models the conversion gain from X-ray photons to electrons, and $\sigma^{2}=5^{2}$ is the variance of electronic noise [56]. We compare the image reconstruction quality obtained with PWLS-ST and PWLS-ULTRA with those of:

• •

FBP: conventional FBP method with a Hanning window.

• •

PWLS-EP: PWLS reconstruction with the edge-preserving regularizer $\mathsf{R}(\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 [53], and ${\varphi}{(t)}\triangleq\delta^{2}(|t/\delta|-\log(1+|t/\delta|))$. We optimized this PWLS cost function using the relaxed OS-LALM [50].

• •

PWLS-DL: PWLS reconstruction with a learned overcomplete synthesis dictionary based regularization, whose image update step is optimized by relaxed OS-LALM instead of the SQS-OS used in [31].

Section IV-F reports the reconstructions from helical CT clinical data of the chest and abdomen (low-dose). Finally, Section IV-G compares the performance of PWLS-ULTRA to an oracle scheme that uses cluster memberships estimated directly from the reference or ground truth images.

To compare various methods quantitatively for the case of the XCAT phantom, we calculated the Root Mean Square Error (RMSE) and Structural Similarity Index Measurement (SSIM) [57] of the reconstructions in a region of interest (ROI). RMSE in Hounsfield units 222Modified Hounsfield units, where air is [math] HU and water is $1000$ HU. (HU) is defined as RMSE $=\sqrt{\sum_{i=1}^{N_{p,\text{ROI}}}(\hat{x}_{i}-x^{*}_{i})^{2}/{N_{p,\text{ROI}}}}$, where $\mathbf{x}^{*}$ is the ground truth image and $N_{p,\text{ROI}}$ is the number of pixels (voxels) in the ROI. Unless otherwise noted, we tuned the parameters of various methods for each experiment to achieve good RMSE and SSIM. For the clinical chest and low-dose abdomen data, the reconstructions were evaluated visually using voxel profiles. We display all reconstructions in this section using a display window $[800,1200]$ HU, unless otherwise noted.

In the 2D fan-beam CT experiments, we pre-learned square transforms and union of square transforms from $8\times 8$ overlapping image patches extracted from five $512\times 512$ XCAT phantom slices, with a patch stride $1\times 1$. We ran $1000$ iterations of the alternating minimization transform learning algorithm in Section III-A (or in [36] when $K=1$) to ensure convergence, and used $\lambda_{0}=31$. The transforms were initialized with the 2D DCT, and k-means clustering (of patches) was used to initialize the clusters for learning a union of transforms. We simulated a 2D fan-beam CT scan using an $840\times 840$ XCAT phantom slice (air cropped) that differs from the training slices, and $\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 $1\times 10^{4}$ and $5\times 10^{3}$ incident photons per ray and no scatter, respectively. We reconstructed a $420\times 420$ image 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.

In the 3D cone-beam CT reconstruction experiments, we pre-learned STs and union of square transforms from $8\times 8\times 8$ patches ($N^{\prime}\approx 1\times 10^{6}$) extracted from a $420\times 420\times 54$ XCAT phantom (air cropped) with a patch stride $2\times 2\times 2$. We set $\lambda_{0}$ large enough, e.g., $\lambda_{0}=31$, to ensure well-conditioned learned transforms. We ran the alternating minimization transform learning algorithms for $1000$ iterations. The transforms were initialized with the 3D DCT, and a random initialization was used for the clusters (because k-means produced some empty clusters for large $K$) for learning a union of square transforms. We simulated an axial cone-beam CT scan using an $840\times 840\times 96$ XCAT phantom with $\Delta_{x}=\Delta_{y}=0.4883$ mm and $\Delta_{z}=0.625$ mm. We generated sinograms of size $888\times 64\times 984$ using GE LightSpeed cone-beam geometry corresponding to a monoenergetic source with $1\times 10^{4}$ and $5\times 10^{3}$ incident photons per ray and no scatter, respectively. We reconstructed a $420\times 420\times 96$ volume with a coarser grid, where $\Delta_{x}=\Delta_{y}=0.9766$ mm and $\Delta_{z}=0.625$ mm. For PWLS-ST and PWLS-ULTRA reconstructions, the patch size was $8\times 8\times 8$ with a patch stride $2\times 2\times 2$ ($\tilde{N}\approx 2\times 10^{6}$ patches). The ROI for the 3D case consisted of the central $64$ of $96$ axial slices and a circular (around center) region in each slice (cylinder in 3D). The diameter of the circle was $420$ pixels, which is the width of each slice.

For the clinical chest data, we reconstructed a $420\times 420\times 222$ image volume (air cropped) with patch size $8\times 8\times 8$ and patch stride $3\times 3\times 3$ ($\tilde{N}\approx 1.5\times 10^{6}$ patches), where $\Delta_{x}=\Delta_{y}=1.1667$ mm and $\Delta_{z}=0.625$ mm, from a helical CT scan. The size of the sinogram was $888\times 64\times 3611$ and pitch was $1.0$ (about $3.7$ rotations with rotation time $0.4$ seconds). The tube current and tube voltage of the X-ray source were $750$ mA and $120$ kVp, respectively. To further evaluate the proposed method, we reconstructed $512\times 512\times 200$ abdomen region volumes with patch size $8\times 8\times 8$, patch stride $3\times 3\times 3$, $\Delta_{x}=\Delta_{y}=1$ mm and $\Delta_{z}=0.625$ mm, from low-dose helical CT patient scans. The size of the sinogram was $888\times 64\times 2952$ and pitch was $1.375$ ($3$ rotations with rotation time $0.8$ seconds). The tube voltage was $120$ kVp, and the tube currents were $150$ mA and $35$ mA (scanned twice for the same patient).

IV-B Parameter Selection

The $\{\tau_{j}\}$ parameters are designed using the $\bm{\kappa}$ information as per (2), so no additional tuning is needed. Since the transforms are pre-learned once from a given dataset and used to reconstruct new data, the parameters $\lambda$ and $\eta$ are tuned during training. As mentioned in prior work [36], the parameter $\lambda$ controls the condition number and larger values of $\lambda$ encourage well-conditioned transforms that work well for image reconstruction. The $\eta$ parameter can be set to achieve low sparsity (e.g., $5-10\%$) and a good trade-off with sparsification error (the transform-domain residual in the training objective) for training data. In our experiments, we learned transforms for a couple different $\eta$ values (training sparsities) and compared their effectiveness in some test reconstructions before picking the best learned model.

During reconstruction, mainly the parameters $\beta$ and $\gamma$ (Section IV-C discusses about $K$) need to be tuned. These parameters are tuned to achieve a good trade-off between image resolution and noise. For example, large values of $\gamma$ would achieve very low sparsities and reduce the noise but potentially oversmooth the image. For a given learned transform, we tuned $\beta$ and $\gamma$ together to achieve good RMSE and SSIM of the reconstruction. Since the PWLS-ST and PWLS-ULTRA formulations are quite similar, except for the richer model and implicit clustering in the latter case, one could tune $\beta$ and $\gamma$ for ST first, and use these optimized values for ULTRA. In our experiments, we tuned parameters separately for ST and ULTRA, and found the tuned values to be typically similar.

Likewise, standard methods like the PWLS-EP method have an overall regularization parameter $\beta$ and an edge-preserving parameter $\delta$, so the number of parameters that one must tune during reconstruction (after training is done) is similar for EP and ULTRA. Similarly as for PWLS-ULTRA, the parameters (maximum patch-wise sparsity level and error threshold for sparse coding) for the prior PWLS-DL were selected carefully (by sweeping over values in a grid) to achieve good RMSE and SSIM in each case, for fair comparison.

IV-C Behavior of the Learning and PWLS-ULTRA Algorithms

We evaluate the behavior of the PWLS-ULTRA method (with $\tau_{j}=1\,\forall j$) for 3D cone-beam CT data with $I_{0}=1\times 10^{4}$. Fig. 2 shows the central slices along three directions for the underlying (true) XCAT phantom volume. We reconstruct the volume from low-dose CT measurements. Fig. 2 shows the RMSE and SSIM of PWLS-ULTRA for various choices of $K$, the number of clusters (patch size $8\times 8\times 8$ and patch stride $2\times 2\times 2$). Rich models (large $K$) produce better reconstructions compared to using a single ST ($K=1$). For the piece-wise constant phantom, $K=5$ clusters works well enough, with only a small additional RMSE or SSIM improvement observed for larger $K$. Larger values of $K$ led to sharper image edges.

Fig. 3 presents an example of the pixel-level clustering in the central axial slice achieved with the PWLS-ULTRA method for $K=5$. Since PWLS-ULTRA clusters patches, we cluster individual pixels using a majority vote among the 3D patches that overlap the pixel. Class $1$ contains most of the soft tissues; class $2$ comprises most of the bones and blood vessels; classes $3$ and $4$ have some high-contrast edges oriented along specific directions; and class $5$ mainly includes low-contrast edges. Since the clustering step (during both training and reconstruction) is unsupervised, i.e., different anatomical structures were not labeled manually, there are also a few edges with high pixel intensities included in class $2$. The trained (3D) transforms (with $\eta=50$) for each cluster are also displayed in a similar manner as in Fig. 1. The transforms show features (e.g., with specific orientations) that clearly reflect the properties of the patches/tissues in each class.

IV-D 2D LDCT Reconstruction Results and Comparisons

IV-D1 Reconstruction Quality

We evaluate the performance of various algorithms for image reconstruction from low-dose fan-beam CT data. Initialized with FBP reconstructions, we ran the PWLS-EP algorithm for $50$ iterations using relaxed OS-LALM with $24$ subsets, and set $\delta=10$ (HU) and the regularization parameter $\beta=2^{16.0}$ and $\beta=2^{16.5}$ for $I_{0}=1\times 10^{4}$ and $I_{0}=5\times 10^{3}$, respectively. For PWLS-DL, PWLS-ST, and PWLS-ULTRA, we initialized with the PWLS-EP reconstruction, and ran $200$ outer iterations with $2$ iterations of the image update step with $4$ ordered subsets, i.e., $N=2$, $M=4$. For PWLS-DL, we pre-learned a $64\times 256$ overcomplete dictionary from $8\times 8$ patches extracted from five XCAT phantom slices (same slices as used for transform learning) with a patch stride $1\times 1$, using a maximum patch-wise sparsity level of $20$ and an error threshold or tolerance for sparse coding of $10^{-1}$. During reconstruction with PWLS-DL, we used a maximum sparsity level of $25$, an error tolerance of $55$, and a regularization parameter of $7.0\times 10^{4}$ and $6.0\times 10^{4}$ for $I_{0}=1\times 10^{4}$ and $I_{0}=5\times 10^{3}$, respectively. For PWLS-ST and PWLS-ULTRA ($K=15$), we chose $\left(\beta,\gamma,\eta\right)$ for the two incident photon intensities as follows: $\left(2.0\times 10^{5},20,75\right)$ and $\left(1.3\times 10^{5},20,75\right)$ for PWLS-ST ($\tau_{j}=1$); $\left(2.0\times 10^{5},20,125\right)$ and $\left(1.0\times 10^{5},25,125\right)$ for PWLS-ULTRA ($\tau_{j}=1$), and $\left(1.3\times 10^{4},22,125\right)$ and $\left(1.0\times 10^{4},25,125\right)$ for PWLS-ULTRA with the weights $\tau_{j}$.

Table I lists the RMSE and SSIM values for reconstructions with FBP, PWLS-EP, PWLS-DL, PWLS-ST ($\tau_{j}=1$), PWLS-ULTRA ($K=15$, $\tau_{j}=1$), and PWLS-ULTRA ($K=15$) with the weights $\tau_{j}$. The adaptive PWLS methods outperform the conventional FBP and the non-adaptive PWLS-EP. Both PWLS-DL that uses an overcomplete dictionary and PWLS-ULTRA using a union of learned transforms lead to better reconstruction quality than PWLS-ST. Importantly, PWLS-ULTRA achieves comparable or better image quality than PWLS-DL. Table II lists the RMSE values in various ROIs (corresponding to specific tissues) for reconstructions with the six methods. The three zoom-ins from left to right in Fig. 4 correspond to ROI-1 to ROI-3 in Table  II, respectively. ULTRA achieve lower RMSE in most of these ROIs compared to DL. Fig. 4 compares the reconstructions for PWLS-DL and PWLS-ULTRA without the weights $\tau_{j}$ at $I_{0}=1\times 10^{4}$. The ULTRA reconstruction shows fewer artifacts and better clarity of bone and soft tissue edges in the selected ROIs.

IV-D2 Runtimes

To compare the runtimes of various data-driven methods, we ran PWLS-DL, PWLS-ST, and PWLS-ULTRA ($K=15$) (all initialized with the FBP reconstruction) for $200$ outer iterations with $2$ iterations of the image update step and $4$ ordered subsets. For PWLS-ULTRA, we performed the clustering step once every outer iteration. While the total runtime for the 200 iterations (using a machine with two $2.80$ GHz 10-core Intel Xeon E5-2680 processors) was $95$ minutes for PWLS-DL, it was only $20$ minutes for PWLS-ST and $27$ minutes for PWLS-ULTRA. We observed that PWLS-DL and the proposed methods had similar convergence rates, but the latter were much faster per iteration, thus leading to much lower net runtimes. The runtime of PWLS-DL was quite equally dominated by the sparse coding (with OMP [55]) and image update steps, whereas for the transform-based methods, the sparse coding and clustering involving simple closed-form solutions and thresholding operations required negligible runtime. The advantage in runtime was achieved despite using an unoptimized Matlab implementation of PWLS-ST and PWLS-ULTRA, and using an efficient MEX/C implementation for sparse coding with OMP [55] in PWLS-DL. PWLS-DL is far slower for 3D reconstructions with large 3D patches. Hence, we focus our comparisons between the transform learning and dictionary learning-based schemes for 2D LDCT reconstruction.

IV-E Low-dose Cone-beam CT Results and Comparisons

We evaluate the performance of various algorithms for reconstructing CT volumes from simulated low-dose cone-beam data. Initialized with FDK reconstructions, we ran the PWLS-EP algorithm with edge-preserving parameter $\delta=10$ (HU) and regularization parameter $\beta=2^{14.5}$ for $50$ iterations with $24$ subsets for both $I_{0}=1\times 10^{4}$ and $I_{0}=5\times 10^{3}$. We evaluate PWLS-ST and PWLS-ULTRA without the patch-based weights. We also evaluate PWLS-ULTRA with such weights. Initialized with the PWLS-EP reconstruction, we ran $2$ iterations of the image update step for the proposed methods with $4$ subsets. We performed the clustering step once every $20$ outer iterations, which worked well and saved computation. We chose $\left(\beta,\gamma,\eta\right)$ for $I_{0}=1\times 10^{4}$ and $I_{0}=5\times 10^{3}$ as follows: $\left(2.0\times 10^{5},18,50\right)$ and $\left(1.5\times 10^{5},20,50\right)$ for PWLS-ST ($\tau_{j}=1$); $\left(2.5\times 10^{5},18,75\right)$ and $\left(1.5\times 10^{5},20,75\right)$ for PWLS-ULTRA ($\tau_{j}=1$); and $\left(1.5\times 10^{4},18,75\right)$ and $\left(1.2\times 10^{4},20,75\right)$ for PWLS-ULTRA with the weights $\tau_{j}$.

Table III lists the RMSE and SSIM values of the reconstructions with FDK, PWLS-EP, PWLS-ST ($\tau_{j}=1$), PWLS-ULTRA ($K=15$, $\tau_{j}=1$), and PWLS-ULTRA ($K=15$) with patch-based weights $\tau_{j}$. Both PWLS-ST and PWLS-ULTRA significantly improve the RMSE and SSIM compared to FDK and the non-adaptive PWLS-EP. Importantly, PWLS-ULTRA with a richer union of learned transforms leads to better reconstructions than PWLS-ST with a single learned ST. Incorporating the patch-based weights in PWLS-ULTRA leads to further improvement in reconstruction quality compared to PWLS-ULTRA with uniform weights $\tau_{j}=1$ for all patches. In particular, the patch-based weights lead to improved resolution for soft tissues in 3D LDCT reconstructions.

Fig. 5 shows the reconstructions and the corresponding error (magnitudes) images (shown for the central axial, sagittal, and coronal planes) for FDK, PWLS-EP, and PWLS-ULTRA ($K=15$) with the patch-based weights. Compared to FDK and PWLS-EP, PWLS-ULTRA significantly improves image quality by reducing noise and preserving structural details (see zoom-ins). Fig. 6 shows the RMSE for each axial slice in the PWLS-EP and PWLS-ULTRA (with the weights $\tau_{j}$) reconstructions. PWLS-ULTRA clearly provides large improvements in RMSE for many slices, with greater improvements near the central slice.

IV-F Results for Clinical Data: Chest and Abdomen Scans

We reconstructed the chest volume from helical CT data. For PWLS-EP, we used the same parameter settings as used for this data in prior work [50]. Initializing with the PWLS-EP reconstruction, we ran the PWLS-ULTRA ($K=5$) method with the weights $\tau_{j}$ for $78$ outer iterations with $3$ iterations of the image update step and $4$ subsets. We performed clustering once every $10$ outer iterations. We chose $\beta=2\times 10^{5}$ and $\gamma=25$ for PWLS-ULTRA to obtain good visual quality of the reconstruction. We used the transforms learned from the XCAT phantom volume with $\eta=100$ to obtain reconstructions with PWLS-ULTRA for the clinical chest CT data. The supplement shows that transforms learned from the XCAT phantom provide similar visual reconstructions as transforms learned from the PWLS-EP reconstruction of the chest data. This suggests that the transform learning algorithm may extract quite general and effective image features without requiring a very closely matched training dataset, which is a key distinction from the PICCS and ndiNLM-type methods [8, 9, 10, 11, 12, 13].

Fig. 7 shows the reconstructions (shown for the central axial plane in the 3D volume) for FDK (provided by GE Healthcare), PWLS-EP (corresponds to Fig. 8(a)), and PWLS-ULTRA with $K=5$ (corresponds to Fig. 9(a)). The PWLS-ULTRA reconstruction has lower artifacts and noise. Moreover, the image features and edges are better reconstructed by PWLS-ULTRA than by PWLS-EP or FDK.

Fig. 8 shows the reconstructions (shown for the central axial, sagittal, and coronal planes in the 3D volume) for PWLS-EP with different regularization strengths $\beta$, denoted as a multiplicative factor of the parameter value in Fig. 7. Fig. 9 shows the reconstructions for PWLS-ULTRA (with patch-based weights) with different parameter combinations. For the sagittal and coronal planes, we show the central $135$ out of $222$ axial slices. Larger regularization strengths $\beta$ would achieve more noise reduction but simultaneously lower spatial resolution in PWLS-EP and PWLS-ULTRA, e.g., compare Fig. 8 and Figs. 9(a) and (d). Larger values of $\gamma$ would achieve lower sparsities and more noise reduction but potentially oversmooth the image, e.g., compare Figs. 9(c) and (d). Small values of $\gamma$ may introduce additional spurious noise in the PWLS-ULTRA reconstruction (compare Figs. 9(a) and (b)). Fig. 11 shows profiles of chest reconstructions (plotted from the central axial slice) for the PWLS-EP and PWLS-ULTRA methods. The profile locations are shown in green lines in Fig. 7. Both PWLS-EP with regularization strength $2$X and PWLS-ULTRA (with patch-based weights) in Fig. 9(a) have lower noise than the PWLS-EP with regularization strength $1$X. Though the spatial resolution of PWLS-EP with regularization strength $2$X is close to PWLS-ULTRA in the selected soft-tissue regions, PWLS-ULTRA reconstructs bone and spine areas with higher resolution, and preserves small features better (compare the zoomed-in areas in Fig. 8 and Fig. 9).

We reconstructed the abdomen volume from low-dose helical CT data. With an initialization of zeros, we ran the PWLS-EP algorithm with $\beta=2^{18.0}$ and $\beta=2^{19.0}$ for $20$ iterations with $12$ subsets for the $150$ mA and $35$ mA scans, respectively. For PWLS-ULTRA, we chose $\beta=1\times 10^{5},\gamma=25$ for the $150$ mA scan, $\beta=1.5\times 10^{5},\gamma=30$ for the $35$ mA scan, and ran it for $50$ outer iterations. The other parameter settings and the transform were the same as those used for the chest scan.

Fig. 10 shows the reconstructions (shown for the central axial, sagittal, and coronal planes in the 3D volume) for PWLS-EP and PWLS-ULTRA with patch-based weights ($K=5$) from low-dose abdomen scans. For the sagittal and coronal planes, we show the central $160$ out of $200$ axial slices. The supplement provides PWLS-EP reconstructions with different regularization strengths. The PWLS-ULTRA reconstructions in Fig. 10 have reduced noise as well as higher resolution, better structural details and shaper image edges than the PWLS-EP results. These results are further example of the potential performance of the proposed PWLS-ULTRA method in clinical settings.

IV-G Comparison to Oracle Clustering Scheme

We consider the 3D cone-beam CT data in Section IV-E with $I_{0}=1\times 10^{4}$, and compare the PWLS-ULTRA ($K=15$) method without patch-based weights to an oracle PWLS-ULTRA scheme without patch-based weights, where the cluster memberships are pre-determined (and fixed during reconstruction) by performing the sparse coding and clustering step (with the learned transforms) on the patches of the reference or ground truth volume. The oracle scheme thus uses the best possible estimate of the cluster memberships. Otherwise, we used the same parameters for the two cases. Fig. 12 compares the reconstructions for the two cases. The proposed PWLS-ULTRA underperforms the oracle scheme by only $1.7$ HU. The more precise clustering leads to sharper edges for the latter method. This also suggests that there is room for potentially improving the proposed clustering-based PWLS-ULTRA scheme, which could be pursued in future works.

V Conclusions

We presented the PWLS-ST and PWLS-ULTRA methods for low-dose CT imaging, combining conventional penalized weighted least squares reconstruction with regularization based on pre-learned sparsifying transforms. Experimental results with 2D and 3D axial CT scans of the XCAT phantom and 3D helical chest and abdomen scans show that for both normal-dose and low-dose levels, the proposed methods provide high quality image reconstructions compared to conventional techniques such as FBP or PWLS reconstruction with a nonadaptive edge-preserving regularizer. The ULTRA scheme with a richer union of transforms model provides better reconstruction of various features such as bones, specific soft tissues, and edges, compared to the proposed PWLS-ST. Finally, the proposed approach achieves comparable or better image quality compared to learned overcomplete synthesis dictionaries, but importantly, is much faster (computationally more efficient). We leave the investigation of convergence guarantees and automating the parameter selection for the proposed PWLS algorithms to future work. The field of transform learning is rapidly growing, and we hope to investigate new transform learning-based LDCT reconstruction methods, such as involving rotationally invariant transforms [39], or online transform learning [58, 59], etc., in future work.

VI Acknowledgments

The authors thank GE Healthcare for supplying the helical chest and abdomen data used in this work. The authors also thank Dr. Hung Nien for his feedback.

VII Computing $\mathbf{D}_{\mathsf{R}}$ in Algorithm $1$

Recall the following definition of $\mathbf{D}_{\mathsf{R}}$ in (11):

[TABLE]

When $\{\tau_{j}\}$ are all identical (or by replacing them with $\max_{j}\tau_{j}$ above for a looser majorizer), $\sum_{k=1}^{K}\sum_{j\in C_{k}}\mathbf{P}_{j}^{T}\mathbf{P}_{j}=\sum_{j=1}^{\tilde{N}}\mathbf{P}_{j}^{T}\mathbf{P}_{j}\in\mathbb{R}^{N_{p}\times N_{p}}$ is a diagonal matrix with the diagonal entries corresponding to image voxel locations and their values being the total number of image patches overlapping each voxel. Moreover, if the patches are periodically positioned with a stride of 1 voxel along each dimension and wrap around at image boundaries, then $\sum_{j=1}^{\tilde{N}}\mathbf{P}_{j}^{T}\mathbf{P}_{j}=l\mathbf{I}$. In this case, $\mathbf{D}_{\mathsf{R}}=2\beta\{\max_{k}\lambda_{\max}(\mathbf{\Omega}_{k}^{T}\mathbf{\Omega}_{k})\}\{\max_{j}\tau_{j}\}l\mathbf{I}$. More generally, when $\tau_{j}$ values differ, we compute $\sum_{k=1}^{K}\sum_{j\in C_{k}}\tau_{j}\mathbf{P}_{j}^{T}\mathbf{P}_{j}$ voxel-wise by summing the $\tau_{j}$ values of the patches overlapping the voxel.

VIII Additional Experimental Results

Section IV.E and Table III of [60] compared the performance of various methods for low-dose cone-beam (3D) CT reconstruction, for the XCAT phantom volume. Fig. 13 shows the reconstructions and the corresponding error (magnitudes) images (shown for the central axial, sagittal, and coronal planes) at $I_{0}=5\times 10^{3}$ for FDK, PWLS-EP, PWLS-ST (with $\tau_{j}=1$ $\forall$ $j$), and PWLS-ULTRA ($K=15$) with patch-based weights $\tau_{j}$. PWLS-ULTRA provides better reconstructions and reconstruction errors compared to the conventional FDK and the non-adaptive PWLS-EP. PWLS-ULTRA also outperforms the proposed PWLS-ST scheme, and provides sharper reconstructions of image edges (see zoom-ins).

Recall that in Section IV.F, we used the transforms learned from the patches of the XCAT phantom volume to perform reconstruction of the chest volume from helical CT data. Alternatively, one could learn the transforms from the patches of the PWLS-EP reconstruction of the helical CT data. Fig. 14 shows the union of transforms ($K=5$) learned from $8\times 8\times 8$ patches of the XCAT phantom and the PWLS-EP chest reconstruction, with $\eta=100$. These two union of transforms display some similar types of features, and provide similar visual reconstructions in PWLS-ULTRA (with patch-based weights $\tau_{j}$) in Fig. 14. Thus, the transform learning algorithm extracts quite general and effective sparsifying features for images, without requiring a very closely matched training dataset.

Fig. 15 provides abdomen reconstructions (shown for the central axial, sagittal, and coronal planes) from low-dose (120kVp, 150mA and 35mA) helical CT data for PWLS-EP with different regularization strengths. We have labeled the reconstruction with good trade-off between image resolution and noise in bold for both doses. These images were used to initialize the PWLS-ULTRA reconstructions in Section IV.F.