# Discriminative Prior - Prior Image Constrained
Compressed Sensing Reconstruction for Low-Dose CT
Imaging

Yang Chen, Yong-Jin Liu, Lizhe Xie, Yining Hu, Huazhong Shu, Limin Luo, Libo Zhang, Zhiguo Gui, Gouenou Coatrieux

► To cite this version:

Yang Chen, Yong-Jin Liu, Lizhe Xie, Yining Hu, Huazhong Shu, et al. Discriminative Prior - Prior Image Constrained Compressed Sensing Reconstruction for Low-Dose CT Imaging. Scientific Reports, 2017, 7 (1), pp.13868. 10.1038/s41598-017-13520-y. hal-04160191

HAL Id: hal-04160191

https://hal.science/hal-04160191

Submitted on 12 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
SCIENTIFIC REPORTS

OPEN

# Discriminative Prior - Prior Image Constrained Compressed Sensing Reconstruction for Low-Dose CT Imaging

Received: 9 June 2017

Accepted: 25 September 2017

Published online: 24 October 2017

Yang Chen<sup>1,2,3</sup>, Jin Liu<sup>1,2,3</sup>, Lizhe Xie<sup>4</sup>, Yining Hu<sup>1,2,3</sup>, Huazhong Shu<sup>1,2,3</sup>, Limin Luo<sup>1,2,3</sup>, Libo Zhang<sup>5</sup>, Zhiguo Gui<sup>7</sup> & Gouenou Coatriex<sup>6</sup>

X-ray computed tomography (CT) has been widely used to provide patient-specific anatomical information in the forms of tissue attenuation. However, the cumulative radiation induced in CT scan has raised extensive concerns in recently years. How to maintain reconstruction image quality is a major challenge for low-dose CT (LDCT) imaging. Generally, LDCT imaging can be greatly improved by incorporating prior knowledge in some specific forms. A joint estimation framework termed discriminative prior-prior image constrained compressed sensing (DP-PICCS) reconstruction is proposed in this paper. This DP-PICCS algorithm utilizes discriminative prior knowledge via two feature dictionary constraints which built on atoms from the samples of tissue attenuation feature patches and noise-artifacts residual feature patches, respectively. Also, the prior image construction relies on a discriminative feature representation (DFR) processing by two feature dictionary. Its comparison to other competing methods through experiments on low-dose projections acquired from torso phantom simulation study and clinical abdomen study demonstrated that the DP-PICCS method achieved promising improvement in terms of the effectively-suppressed noise and the well-retained structures.

X-ray computed tomography (CT) is widely used in diagnostic imaging, image-guided surgeries and radiother-
apy tasks for clinic applications<fnref n="1" />⁻³. Recent years, more and more attention has been paid to the risk of ionizing
radiation received during CT scans in recent years. CT dose levels should be kept ALARA (as low as reasonably
achievable) if enough diagnostic information is provided. Reducing the exposure level (low-exposure protocol)
can be simply implemented to lower doses in clinic application⁴⁻⁶. Conventional analytical reconstruction algo-
rithms (e.g. Filter Back Projection (FBP) method, Feldkamp-Davis-Kress (FDK) method⁷) often lead to unsat-
isfactory diagnostic CT image quality with increased mottle noise and streak artifacts due to the deteriorated
projections.

Derived from statistical models of projection data, iterative methods have demonstrated better performance
in noise-artifacts suppression in reconstruction, and lead to improved image quality. With the rapid develop-
ment of compressed sensing (CS)⁸,⁹, sparsity-promotion (e.g., total variation (TV)¹⁰⁻¹³, tight frame (TF)¹⁴, wavelet
transform¹⁵, and dictionary learning¹⁶⁻²³) based constraint had been considered to alleviate the ill-posedness
in low-dose CT (LDCT) reconstruction. For example, the TV minimization based constraint (or regulariza-
tion) reconstruction, which was developed based on the sparse distribution of discrete gradient. Such constraint
has demonstrated its effectiveness in preserving edges and suppressing noise¹¹⁻²³. Recently, dictionary based
sparse representation was introduced to improve medical imaging¹⁶⁻²³. In¹⁶, Xu. *et al.* proposed a dictionary

<footnote n="1">¹Laboratory of Image Science and Technology, School of Computer Science and Engineering, Southeast University, Nanjing, 210096, China. ²Key Laboratory of Computer Network and Information Integration (Southeast University), Ministry of Education, Nanjing, 210096, China. ³International Joint Research Laboratory of Information Display and Visualization, Southeast University, Ministry of Education, Nanjing, 210096, China. ⁴Jiangsu Key Laboratory of Oral Diseases, Nanjing Medical University, No. 140, Hanzhong road, Nanjing, 210000, China. ⁵The General Hospital of Shenyang Military, Shenyang, Liaoning, 110016, China. ⁶The Institut Mines-Telecom, Telecom Bretagne; INSERM U1101 LaTIM, Brest, 29238, France. ⁷The National Key Laboratory for Electronic Measurement Technology, North University of China, Taiyuan, 030051, China. Correspondence and requests for materials should be addressed to H.S. (email: shu.list@seu.edu.cn) or L.Z. (email: zlb19782002@163.com)</footnote>
learning (DL) based method for LDCT reconstruction and demonstrated performance superior to the TV based method. Containing rich feature information, the composing atoms in dictionary work well in representing structured features, which lead to significantly enhanced reconstruction quality. In¹⁹, Liu. et al. proposed an iterative reconstruction by feature constrained to improve LDCT imaging, which relies on a pre-defined 3D feature dictionary constructed from standard-dose CT (SDCT) sample. The pre-defined 3D feature dictionary contains patient-specific anatomical structures, which have significant morphological discrepancy with the undesirable noise-artifacts features. This 3D feature constraint approach has good potential for suppressed noise-artifacts and the better retained anatomical structures.

Prior knowledge can also be utilized by directly including some available high quality prior image into the cost function, which leads to the algorithms termed prior image constrained compressed sensing (PICCS)²⁴. Some PICCS type algorithms have gotten successful applications in interventional imaging, treatment monitoring, and 4D cardiac reconstruction²⁴⁻²⁷. Nevertheless, the prior image term takes effect in PICCS algorithms via the minimization of subtracted residual function, which requires an exact position correspondence between the prior image and the current image. This requirement, however, is often hard to meet due to the unavailability of precisely matched high quality prior images, which often greatly limits the practical feasibility.

This paper presents a reconstruction approach termed discriminative prior - prior image constrained compressed sensing (DP-PICCS) for LDCT, which improves current PICCS reconstruction by utilizing discriminative prior knowledge. The discriminative prior knowledge works through feature dictionaries which containing rich tissue attenuation feature information and noise-artifacts residual feature information. The proposed DP-PICCS algorithm can be easily implemented via an alternative optimization scheme with a good parameter robustness. Experiments with simulated torso phantom and clinical abdomen data were conducted though the comparison with three state-of-art reconstruction algorithms. In summary, this paper is structured as follows: in section II, we present a detailed description of the proposed method and the implementation procedure for LDCT reconstruction. Experimental results and quantify the performance from simulated torso phantom and clinical abdomen projection are given and parameters setting of our proposed method are discussed in section III. Finally, conclusions the paper and plans for future work are sketched in section IV.

## Method

### Discriminative prior - prior image constrained compressed sensing model (DP-PICCS)
The standard CT reconstruction problem can be considered as an inverse problem. Directly solving the it may not be feasible due to the measurement incompletion or noise contamination in projection data. Approaches based on sparse representation can be used to overcome the ill-posedness by incorporating representation related knowledge in the form of a regularization term¹¹,²⁸. The formula for CT image reconstruction with a constraint term can be expressed as:

$$ 
\begin{array}{l}
u^* = \arg \min_{u} \lambda R_g(u) \\
\text{s.t.} \quad (Gu - p)^T W(Gu - p) = 0
\end{array}
\qquad (1) $$

where $G = \{g_{mn}\} \in R^{M \times N}$ is the system matrix representing the contribution of the n-th voxel to the m-th X-ray path; $u = \{u_n\} \in R^N$ denotes the discrete vector of target image reconstructed from the projection data $p = \{p_m\} \in R^M$. N is the total voxel number of the volume to reconstruct; W is projection noise statistics weight matrix, $R_g$ denotes some specific global constraint term. Specifically, further reconstruction improvement can be brought via the PICCS approach by directly incorporating a high-quality prior image $u_{prior}$ into cost function²⁴. This is accomplished by incorporating an image similar to that which we want to reconstruct into the reconstruction procedure. The PICCS algorithm can be formulated as a constrained optimization procedure as follows:

$$ 
\begin{array}{l}
u^* = \arg \min_{u} \left[ \lambda R_p(u - u_{prior}) + (1 - \lambda) R_g(u) \right] \\
\text{s.t.} \quad (Gu - p)^T W(Gu - p) = 0
\end{array}
\qquad (2) $$

Here, $u_{prior}$ is the vector of attenuation coefficients of prior image and $\lambda \in [0, 1]$ is the weight parameter of prior image constraint terms in the objective function. The first prior image constraint term $R_p$ in Eq. (2) constrains the reconstructed image toward the prior image $u_{prior}$, while the second one $R_g$ works as a global constraint term to overcome the ill-posedness in LDCT. Effect of the two constraints is modulated through the parameter $\lambda$ in Eq. (2).

Nevertheless, two limits exist for the PICCS model in Eq. (2). The first one is the position displacement between a prior image and the current target image. A high quality prior image with exactly matched position correspondence is often not available. The second one is that the unitary analysis transform model (i.e. global and prior image constraint are both TV constraint) is used in most PICCS methods. The global constraint $R_g$ and prior image constraint $R_p$ give penalization with respect to the information on anatomical structures and residuals, which inherently have significant morphological discrepancy. The L1 norm based TV constraint is not effective in tackling with the versatile features in practical reconstructions. Some blocky artifacts and smeared details tend to appear in the reconstructed LDCT images for the PICCS reconstruction with TV constraints.

In this study, to overcome these limits, we improve PICCS reconstructions by imposing two discriminative feature dictionary constraints built specifically for the desirable tissue attenuation features and the noise-artifacts residual features in the terms in the PICCS framework. The feature dictionary constraints to solve²⁹⁻³²:
<figure><img src="image_2.png" /><figcaption>Figure 1. (a) LDCT abdomen image; (b) Tissue attenuation component image after DFR; (c) Residual component after DFR.</figcaption></figure>

$$ \alpha^* = \arg \min_{\alpha} \sum_s (\|E_s u - D\alpha_s\|_2^2 + \nu \|\alpha_s\|_0) \quad (3) $$

where D is the pre-defined feature dictionary, and $E_s$ denotes an operator to extract the s-th 3D patch in the volume space, $\alpha_s$ is the sparse representation coefficients of $E_s u$ and $\nu$ is the Lagrange multiplier. Then the corresponding developed DP-PICCS framework is formulized in the following minimization problem:

$$ \begin{aligned} \{u^*, \alpha_r^*, \alpha_t^*\} = & \operatorname*{argmin}_{u, \alpha_r, \alpha_t} \left[ \lambda \sum_s \left( \|E_s(u - u_{\text{prior}}) - D^r \alpha_{r_s}\|_2^2 + \nu_r \|\alpha_{r_s}\|_0 \right) \right. \\ & \left. \qquad + (1 - \lambda) \sum_s \left( \|E_s u - D^t \alpha_{t_s}\|_2^2 + \nu_t \|\alpha_{t_s}\|_0 \right) \right] \\ \text{s.t.} \quad &(Gu - p)^T W(Gu - p) = 0 \end{aligned} \quad (4) $$

where $D^r$ and $D^t$ are the pre-defined feature dictionaries which composed of atoms featuring the anatomy structures for normal tissue attenuation and the noise-artifacts in the reconstruction image, $\alpha_{r_s}$ and $\alpha_{t_s}$ are the sparse coding, $\nu_r$ and $\nu_t$ are the Lagrange multiplier. Here, the **discriminative prior (DP)** is reflected in two aspects: the discriminative representation with the dictionaries $D^r$ and $D^t$, and the high-quality prior image $u_{\text{prior}}$ generated by the discriminative feature representation (DFR) method in$^{20}$.

**Implementation of the DP-PICCS method.** The prior image construction. In the implementation of the DP-PICCS method, the initial iteration is set to the FDK reconstruction, and then the DFR approach in$^{20}$ is applied as the post-processing method to generate the prior image. The discriminative composite dictionary $D = [D^r, D^t]$ is composed of a tissue attenuation feature dictionary $D^r$ and noise-artifacts residual dictionary $D^t$. It is noted that the atoms in $D$ need not to be inter-independent, which can provide over-complete representations of CT image features. All the overlapped patches in the 3D LDCT volume $u$ are represented by the linear combination of the atoms in dictionary $D$ with the linear coefficient vectors $\alpha = [\alpha_t, \alpha_r]^T$ calculated by the orthogonal matching pursuit (OMP) algorithm. After this, we obtain the processed volume (approximate tissue attenuation featured volume $\tilde{u}_t$) by DFR: $\tilde{u}_t = \sum_s (E_s^T (D^t \alpha_t))/\sum_s (E_s^T E_s)^{20}$.

Figure 1(a) depicts one slice of clinical LDCT image from one GE Discovery HD750 CT unit. Figure 1(b) and (c) are the corresponding DFR processed image and the noise-artifacts residual component. We can see the DFR approach works well in separating undesirable noise-artifacts features from the original LDCT images, and has the potential to provide a high quality prior image ($u_{\text{prior}} = \tilde{u}_t$).

**Iteration reconstruction.** With the initial prior image $u_{\text{prior}}^0$, tissue attenuation feature dictionary $D^r$ and noise-artifacts residual feature dictionary $D^t$, Eq. (4) is in fact a non-convex $L_0$-norm constrained optimization procedure with respect to $u$. We solve Eq. (4) using an alternative iteratively of image updating task and sparse coding task$^{33}$. Using the FDK reconstruction as the initial volume $u^0$, the implementation of the DP-PICCS method includes the following three steps:

(1) Image update:

$$ u^{k+1} = \arg \min_{u} \frac{\beta}{2} (Gu - p)^T W(Gu - p) + \left[ \lambda \sum_{s} \|E_s(u - u_{\text{prior}}^n) - D^r \alpha_{r_s}^k\|_2^2 + (1 - \lambda) \sum_{s} \|E_s u - D^t \alpha_{t_s}^k\|_2^2 \right] \quad (5) $$

where the superscript $k$ and $n$ are the update index, $\beta$ is the global parameter for the fidelity term.

(2) Tissue attenuation feature sparse coding:
<figure><img src="image_3.png" /><figcaption>Figure 2. A flowchart of the proposed DP-PICCS algorithm for LDCT reconstruction.</figcaption></figure>

$$ \alpha_{t_s}^{k+1} = \arg \min_{\alpha_{t_s}} \left\| E_s u^{k+1} - D^t \alpha_{t_s} \right\|_2^2 \quad \text{s.t.} \quad \left\| \alpha_{t_s} \right\|_0 \le T_t, \quad \forall s \qquad (6) $$

(3) Residual feature sparse coding:

$$ \alpha_{r_s}^{k+1} = \arg \min_{\alpha_{r_s}} \left\| E_s(u^{k+1} - u_{\text{prior}}^n) - D^r \alpha_{r_s} \right\|_2^2 \quad \text{s.t.} \quad \left\| \alpha_{r_s} \right\|_0 \le T_r, \quad \forall s \qquad (7) $$

where $T_r$ and $T_t$ are the sparsity level parameters which limiting the maximum atom numbers in the two dictionaries used for a 3D patch sparse coding.

The sup-problem of Eq. (5) is a typical quadratic form. By the separable paraboloid surrogate method, it can be optimized as²⁸:

$$ u^{k+1} = u^k - \frac{G^T W (Gu^k - p) + \frac{1}{\beta} \frac{\partial R(u)}{\partial u} \Big|_{u=u^k}}{G^T W G I + \frac{1}{\beta} \frac{\partial^2 R(u)}{\partial u^2} \Big|_{u=u^k}} \qquad (8) $$

where $R(u) = \lambda R_p(u - u_{\text{prior}}) + (1 - \lambda) R_g(u) = \lambda \sum_s \left\| E_s(u - u_{\text{prior}}^n) - D^r \alpha_{r_s}^k \right\|_2^2 + (1 - \lambda) \sum_s \left\| E_s u - D^r \alpha_{r_s}^k \right\|_2^2$, I is a unity vector. Then, the solution of volume u becomes:

$$ u^{k+1} = u^k - \frac{G^T W (Gu^k - p) + \frac{2(1-\lambda)}{\beta} \sum_s E_s^T (E_s u^k - D^t \alpha_{t_s}^k) + \frac{2\lambda}{\beta} \sum_s E_s^T (E_s u^k - u_{\text{prior}}^n) - D^r \alpha_{r_s}^k}{G^T W G I + \frac{2}{\beta} \sum_s E_s^T E_s} \qquad (9) $$

The sparse coefficients $\alpha_{t_s}^{k+1}$ and $\alpha_{r_s}^{k+1}$ in Eqs (6) and (7) can be solved by the greedy strategy based BP, MP or OMP algorithm²⁹⁻³²,³⁴. In this paper, the Batch-OMP algorithm is considered as the solver for improving high computational efficiency³⁵.

The overall DP-PICCS algorithm is implemented based on Algorithm 1, which contains two loops: the inside loop labeled by n for prior image update (the total inside iteration number is $N_{\text{max}}$), and the outside loop labeled by k. To get a stable solution, the whole reconstruction proceeds until the image update between two consecutive
<figure><img src="image_4.png" /><figcaption>Figure 3. Outline of the feature dictionaries construction.</figcaption></figure>

**Algorithm 1.** Discriminative prior and prior image constrained compressed sensing (DP-PICCS) for LDCT imaging

**Input:** projection data $p$; tissue attenuation dictionary $D^t$ and noise-artifacts residual dictionary $D^r$;

**Output:** prior image $u_{prior}$ and reconstruction image $u^*$;

**Reconstruction:**

**Initialization:** $\lambda, \beta, T_t, T_r, u^0 = \text{FDK}(p), u_{prior}^0$

**Loop:** for $k = 1, 2, \dots, K_{\text{max}}$ **do**

1) Image update: updating $u^k$ via Eq.(9);

2) Tissue attenuation feature sparse coding: updating $\alpha_{t_v}^k$ via Eq. (6);

3) Residual feature sparse coding: updating $\alpha_{r_v}^k$ via Eq. (7);

4) **IF:** update prior image

Apply DFR method to update the prior image $u_{prior}^n$;

**End IF**

**End for**

iterations falls below a pre-defined threshold $\theta$ ($\|u^{k+1} - u^k\|_2 / \|u^k\|_2 \le \theta$). In practice, the threshold parameter $\theta$ is hard to be normalized for different objective function and the measured projection. So, in this study we just use $K_{\text{max}}$ and $N_{\text{max}}$ as the simply stopping iteration numbers for the outer loop and inner loop, respectively. The full flowchart of the proposed DP-PICCS algorithm is shown in Fig. 2.

**Construction of the Feature Dictionaries.** The DFR approach considers a LDCT volume $u$ as the sum of a volume with normal tissue attenuation features and a volume with noise-artifacts residual features$^{20}$. So atoms with distinctly different features should be used to represent the different image constraint terms. With an aim to giving specific sparse representation, we build the two dictionaries using the SDCT projections acquired from the same CT scanner as the target LDCT images to reconstruct, in which the same scan protocol except the tube current is used. Performance of the proposed DP-PICCS method is highly determined by tissue attenuation feature dictionary $D^t$ and residual feature dictionary $D^r$, which are respectively the collections of the atoms learned from sample patches. Schemed in Fig. 3, this strategy of dictionary construction includes the following two steps:

**Step 1, Sample patches preparation:** From the same CT scanner geometry as the LDCT images to reconstruct, we acquire the SDCT image samples $V_{sd}$. Then, from this SDCT images, a set of LDCT images $V_{id}$ are simulated via Eq. (10) and FDK reconstruction with ramp kernel.

$$N_m = \text{Poisson}\{N_{0m} e^{-p_m} + r_m\}, m = 1, 2, \dots, M \quad (10)$$

where $N_m, N_{0m}$ and $r_m$ are the numbers of transmitted photons, incident x-ray photons and read-out noise along the $m$-th X-ray path, respectively. $P_m$ is the attenuation integral of the $m$-th X-ray path. Note that a lower number of incident photons $N_{0m}$ is used to ensure a rich information of residual features in Eq. (10). We then extract tissue attenuation features patches $F_t$ from the SDCT volume $V_{sd}$, and the extract residual patches $F_r$ from the difference volume $V_r$ between the two matched SDCT and LDCT. To avoid great many training samples for calculation, the extraction operation is performed with a 3D interval $\Delta_D$ specified by the intervals $\Delta_x, \Delta_y$, and $\Delta_z$ in the three directions (axes $x, y$ and $z$).

**Step 2, Dictionary construction:** Here, the dictionaries $D^t$ or $D^r$ are trained from the $F_t$ or $F_r$ by solving the following minimization problem via K-Singular Value Decomposition (K-SVD) method$^{29-32}$:
<table><caption>**Table 1.** Parameter Settings Of Different Iterative Reconstruction Methods.</caption><thead><tr><th>Dataset</th><th>PICCS<sub>DFR</sub></th><th>GDSIR</th><th>DP-PICCS</th></tr></thead><tbody><tr><td>Phantom Case D1</td><td>λ = 0.4, β = 4 × 10<sup>-3</sup></td><td>λ = 1.4, L = 8, ε = 1.6 × 10<sup>-3</sup></td><td>λ = 0.36, β = 0.24, T<sub>i</sub> = 8, T<sub>a</sub> = 6</td></tr><tr><td>Phantom Case D2</td><td>λ = 0.3, β = 2 × 10<sup>-3</sup></td><td>λ = 2.3, L = 8, ε = 2.6 × 10<sup>-3</sup></td><td>λ = 0.35, β = 0.21, T<sub>i</sub> = 8, T<sub>a</sub> = 8</td></tr><tr><td>Phantom Case D3</td><td>λ = 0.3, β = 1.4 × 10<sup>-3</sup></td><td>λ = 2.3, L = 8, ε = 3.6 × 10<sup>-3</sup></td><td>λ = 0.35, β = 0.19, T<sub>i</sub> = 8, T<sub>a</sub> = 8</td></tr><tr><td>Clinical Abdomen Data A and Data B</td><td>λ = 0.3, β = 2.5 × 10<sup>-3</sup></td><td>λ = 1.8, L = 8, ε = 3.6 × 10<sup>-3</sup></td><td>λ = 0.35, β = 0.2, T<sub>i</sub> = 8, T<sub>a</sub> = 8</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** The constructed 3D tissue attenuation feature dictionary and residual feature dictionary for experiments.</figcaption></figure>

$$ \min_{D', \Gamma_t} \|F_t - D^t \Gamma_t\|_2^2 \quad \text{s.t.} \quad \|\Gamma_t\|_{row-0} \le T_L \qquad (11) $$

$$ \min_{D', \Gamma_r} \|F_r - D^r \Gamma_r\|_2^2 \quad \text{s.t.} \quad \|\Gamma_r\|_{row-0} \le T_L \qquad (12) $$

where $\Gamma_t$ and $\Gamma_r$ are the sparse coefficient arrays for each block in $F_t$ and $F_r$, and are alternatively updated with the dictionary atoms in K-SVD algorithm. $T_L$ is the sparsity level limiting the term $\|\cdot\|_{row-0}$, a pseudo-norm counting the number of non-zero coefficients in each block representation. It should be noted that the atom number of the learned dictionaries is much smaller than the block number in $F_t$ and $F_r$.

## Experiment and Results

To evaluate the performance of the proposed DP-PICCS method, the FDK method (with ramp filter), the DFR post processing method (DFR-post)²⁰, the iterative TV minimization based reconstruction¹³, the PICCS method with TV constraint and DFR-post prior image (PICCSDFR) and a global feature dictionary-based statistical iterative reconstruction (GDSIR) approach¹⁶ were adopted for comparison.

The traditional PICCS algorithm under the TV constraint framework and DFR-post prior image (PICCSDFR) is formulated as²⁴:

$$ \begin{aligned} u^* = & \arg \min_{u} \left[ \lambda \|TV(u - u_{DFR-post})\|_1 + (1 - \lambda) \|TV(u)\|_1 \right] \\ \text{s.t.} \quad & (Gu - p)^T W(Gu - p) = 0 \end{aligned} \qquad (13) $$

where the sparsifying transform TV(•) is the TV constraint¹².

The GDSIR is a patch-based approach by extracting the prior information via a global dictionary trained from a high-quality CT image, and the associative image reconstruction is equivalent to solve the following minimization problem²⁵:
<figure><img src="image_6.png" /><figcaption>**Figure 5.** (a) The anthropomorphic torso phantom (Kyoto Kagaku co. LTD Japan); (b) Middle slice of this torso phantom volume axial view.</figcaption></figure>

$$ \begin{gathered} \{u^*, \alpha^*\} = \arg \min_{u, \alpha} \left( \sum_s (\|E_s u - D^g \alpha_s\|_2^2 + \nu_s \|\alpha_s\|_0) \right) \\ \text{st. } (Gu - p)^T W(Gu - p) = 0 \end{gathered} \quad (14) $$

In this study, the global dictionary $D^g$ was the same as tissue attenuation dictionary in DFR ($D^g = D^l$), and an alternating minimization scheme described in¹⁶ is applied to find the solution of Eq. (14).

All algorithms were implemented in Matlab 8.3 environment on a personal computer (Intel i7-4790k CPU and 32-GB RAM). The reconstruction parameters used for all the experiments are listed in Table 1. All reconstructions were stopped after 50 iterations and updated to a stable solution. Two dictionaries $D^l$ and $D^r$ were respectively built using four set clinical abdomen SDCT volumes containing total 100 slices acquired from the Siemens CT scanner following the steps in the Fig. 3. The sparsity level $T_L$ was set 8, the atom number $C$ and patch size $B$ in the dictionary construction stage were respectively set to 1000 and $8 \times 8 \times 5$. The built 3D feature dictionaries were illustrated in Fig. 4. We can observe that the atoms in the two dictionaries can well reflect the attenuation features and the noise-artifacts features in CT images. It is also noted that some atoms in the tissue attenuation dictionary $D^l$ present the textures with random intensity distribution, which were actually related to some background regions in the SDCT images when mean values are removed. Though not presenting obvious structural textures, these atoms with random intensity distribution in $D^l$ practically contribute to the restoration of realistic SDCT image textures.

**Anthropomorphic torso phantom simulation study.** *Data Acquisition.* In torso phantom simulation study, a set of high quality CT volume (SDCT images) was scanned at 600mAs protocol on a GE Discovery HD750 CT scanner (tube voltage: 120KVp). The high quality torso phantom volume was reconstructed by FDK method with standard ramp filter. The torso phantom and the middle slice of axial views are illustrated in Fig. 5. For simulation, the source-to-axial distance and source-to-detector distance of the trajectory was respectively set to 57.3 cm and 101 cm. The detector panel has $960 \times 256$ elements with element size $1.024 \times 1.10165 \text{mm}^2$. 360 projections covered one cycle axial scan were uniformly collected. Three different Poisson noise intensity were superimposed onto the raw projection data to synthesize low dose projection data (Photon intensities are $b_m = 5 \times 10^4$, $1 \times 10^4$ and $5 \times 10^3$ in Eq. (10) for Case D1, Case D2 and Case D3, respectively). For all the iterative algorithms, the whole projections were split into 20 subsets to accelerate computation.

**Visual Assessment.** Reconstructions with induced prior image depend heavily on the quality of prior image. From top to bottom, Fig. 6(a1–a3) provides the FDK reconstructed images for the scanner protocols of Case D1, Case D2 and Case D3. Fig. 6(b1–b3) provides the corresponding prior images produced by the DFR-post method from LDCT FDK images. All images are illustrated with the display window center 60HU and display window width 450HU. From the Fig. 6, we can observe the presence of strong noise and artifacts in the FDK reconstruction, which become severer as the doses decrease. From the second row in Fig. 6, the DFR-post method shows good performance in noise-artifacts suppression and structure preservation, and can provide as prior image in good quality and without deformation.

Figure 7 depicts the results from different reconstruction algorithms with three scanner protocols. From top to bottom, the scanner protocols are Case D1, Case D2 and Case D3, respectively. From left to right, the reconstruction algorithms are TV, PICCS<sub>DFR</sub>, GDSIR and DP-PICCS, respectively. We can observe that the smooth regions are plagued by some block and sharp artifacts, which may have resulted from the TV piecewise smooth constraint
<figure><img src="image_7.png" /><figcaption>**Figure 6.** FDK reconstruction images and corresponding DFR-post images. (a1–a3) are the FDK reconstructed images correspond to Case D1, Case D2 and Case D3, respectively; (b1–b3) are the DFR-post images correspond to Case D1, Case D2 and Case D3, respectively.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Selected axial views in the simulated phantom reconstruction. From left to right, images reconstructed using the TV, PICCS<sub>DFR</sub>, GDSIR and DP-PICCS methods from Case D1 (the first row), Case D2 (the second row) and Case D3 (the third row) projection, respectively.</figcaption></figure>

in TV and PICCS<sub>DFR</sub> methods results in Fig. 7(a3) and (b3). This is because the gradient variation based TV constraint fails to provide good discrimination ability between desirable tissue structure features and streak-artifacts. In the third columns in Fig. 7, the GDSIR method works well in suppressing noise and artifacts but at the cost of smoothing out some organ edges (see the zoomed region and yellow arrow in Fig. 7(c3)). Comparing the results in all the zoomed regions, we can see that the proposed DP-PICCS method achieves the best image quality in structure retention and noise-artifact suppression for all the three dose levels.

**PSNR and UQI Measure.** Quantitative evaluation is performed using two metrics the peak signal to noise ratio (PSNR) and universal quality index (UQI). Here the PSNR and UQI are calculated via Eqs (15) and (16)<sup>36</sup>:
<figure><img src="image_9.png" /><figcaption>Figure 8. PSNR and UQI measures on the different reconstruction images in Fig. 7. (a) PSNR; (b) UQI.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. Plots of PSNR and UQI versus DP-PICCS iteration number in Case D1, Case D2 and Case D3. (a) PSNR; (b) UQI.</figcaption></figure>

$$ \text{PSNR}(u, u_{GS}) = 10 \log_{10} \left( \frac{u_{\text{max}}^2}{\frac{1}{N} \sum_{n=1}^{N} (u_n - u_{GSn})^2} \right) \quad (15) $$

$$ \text{UQI}(u, u_{GS}) = \frac{4\sigma_{u,u_{GS}} \overline{u} \overline{u}_{GS}}{(\sigma_u^2 + \sigma_{u_{GS}}^2)[(\overline{u})^2 + (\overline{u}_{GS})^2]} \quad (16) $$

where u represents the reconstructed image from the LDCT, and $u_{GS}$ denotes a golden standard image, n indicates the voxel index and N is the number of voxels, $\overline{u}$ and $\overline{u}_{GS}$ are the mean pixel values of u and $u_{GS}$, respectively, $\sigma_u$ and $\sigma_{u_{GS}}$ are the standard deviation of image pixel values of u and $u_{GS}$, respectively, and $\sigma_{u,u_{GS}}$ is the correlation coefficient between u and $u_{GS}$. The UQI is unitless and has a dynamic range of [-1, 1] reflecting the similarity degree between the reconstructed and golden standard images.

In calculating the PSNR and UQI of the reconstructed images, FDK reconstructed SDCT images (tube voltage: 120KVP, tube current: 600mA, Fig. 5(b)) are used as the golden standard reference. The PSNR and UQI results are shown in Fig. 8. It can be seen in Fig. 8 that the FDK method obtains the worst scores for all the three cases dose levels, and the proposed method performs better than the competing methods in terms of the two metrics (with margins of 1–2 dB for PSNR and 0.01–0.02 for UQI). We can also see that such quantitative results are consistent with the visual performance in Fig. 7. After the LDCT FDK images processed by DFR-post method, the images quality gets significant improvement in both PSNR and UQI. As shown in Fig. 8, the PICCSDFR reconstructed CT images obtained higher quality scores than the TV reconstructed CT images, and such improvement is brought by the use of high quality prior images from the DFR post-processing. The PSNR indexes of GDSIR reconstructed images are higher than those of the TV and PICCSDFR reconstructed images, while the UQI indexes remain in the similar level.

*Convergence Analysis of the DP-PICCS Method.* As described in section II.B, we solve the DP-PICCS reconstruction using split alternating optimization iteration between feature sparse coding and reconstruction image
<figure><img src="image_11.png" /><figcaption>**Figure 10.** The axial views for Data C1. (a,b) are the images reconstructed by the FDK and TV methods for standard dose protocol; (c-h) are the images reconstructed by the FDK, DFR-post, TV, PICCS<sub>DFR</sub>, GDSIR, DP-PICCS methods for low dose protocol; (a1-h1) are the white zoomed regions of the rectangles delineated in (a-h).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** The axial views for Data C2. (a-b) are the images reconstructed by the FDK and TV methods for standard dose protocol; (c-h) are the images reconstructed by the FDK, DFR-post, TV, PICCS<sub>DFR</sub>, GDSIR, DP-PICCS methods for low dose protocol; (a1-h1) are the white zoomed regions of the rectangles delineated in (a-h).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** The sagittal views for Data C1. (a,b) are the images reconstructed by the FDK and TV methods for standard dose protocol; (c-h) are the images reconstructed by the FDK, DFR-post, TV, PICCS<sub>DFR</sub>, GDSIR, DP-PICCS methods for low dose protocol; (a1-h1) are the white zoomed regions of the rectangles delineated in (a-h).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** The coronary views for Data C2. (a-b) are the images reconstructed by the FDK and TV methods for standard dose protocol; (c-h) are the images reconstructed by the FDK, DFR-post, TV, PICCS<sub>DFR</sub>, GDSIR, DP-PICCS methods for low dose protocol; (a1-h1) are the white zoomed regions of the rectangles delineated in (a-h).</figcaption></figure>

update. In image updates step, sub-problem Eq. (5) is a quadric constraint minimization problem and we solved it with a separable paraboloid surrogate method. This step can be considering a convergence process.

In feature sparse coding step, the non-convex sub-problems of sparse coding in Eq. (6) and Eq. (7) are NP-hard problem. To solve them, the Batch-OMP algorithm employs greedy strategy only to obtain a local minimizer in most cases<sup>30,35</sup>. Such alternating process cannot be guaranteed to converge to a global minimum due to the non-convexity L0-norm of the objective function in Eq. (4). To analyze the iteration stability of the DP-PICCS method, the global PSNR and UQI measures on the entire to-be-reconstructed anthropomorphic torso phantom image were calculated. The plotted PSNR and UQI values in Fig. 9 show that the reconstruction quality increases as the iteration proceeds for the proposed algorithm. This observation is confirmed by the study on both torso phantom and clinical abdomen data. Considering the formulation of our problem is similar to that in DLMRI, we can make a similar statement of a stable iteration as in<sup>37</sup>.

**Clinical abdomen data study.** *Data Acquisition.* In the clinical abdomen imaging study, two sets of clinical abdomen projection data (Data C1 and Data C2) were provided by the Mayo Clinic (USA), which scanned from a Somatom Definition AS + CT scanner in a helical mode (tube voltage: 100KVp; tube current: 360mAs)<sup>38</sup>.
<table><caption>**Table 2.** The lSNR (Unit: dB) Values Of The Clinical Data Results.</caption><thead><tr><th rowspan="2">Scanner Protocol</th><th rowspan="2">Reconstruction Methods</th><th colspan="3">Data C1</th><th colspan="3">Data C2</th></tr><tr><th>ROI 1</th><th>ROI 2</th><th>ROI 3</th><th>ROI 1</th><th>ROI 2</th><th>ROI 3</th></tr></thead><tbody><tr><td rowspan="2">SDCT</td><td>FDK</td><td>50.15</td><td>49.23</td><td>54.89</td><td>25.43</td><td>33.87</td><td>38.79</td></tr><tr><td>TV</td><td><strong>90.69</strong></td><td><strong>114.49</strong></td><td><strong>109.95</strong></td><td><strong>111.52</strong></td><td><strong>118.27</strong></td><td><strong>126.54</strong></td></tr><tr><td rowspan="6">LDCT</td><td>FDK</td><td>24.22</td><td>21.94</td><td>28.47</td><td>13.77</td><td>15.60</td><td>16.78</td></tr><tr><td>DFR-post</td><td>72.88</td><td>77.63</td><td>85.91</td><td>42.23</td><td>49.36</td><td>53.19</td></tr><tr><td>TV</td><td>78.73</td><td>91.01</td><td>101.11</td><td>66.28</td><td>99.52</td><td>93.49</td></tr><tr><td>PICCS<sub>DFR</sub></td><td>75.53</td><td>105.05</td><td>97.41</td><td>94.82</td><td>109.31</td><td><strong>136.73</strong></td></tr><tr><td>GDSIR</td><td>79.45</td><td>99.93</td><td>98.29</td><td>87.91</td><td>94.34</td><td>103.54</td></tr><tr><td>DP-PICCS</td><td><strong>87.86</strong></td><td><strong>105.73</strong></td><td><strong>104.59</strong></td><td><strong>95.72</strong></td><td><strong>121.51</strong></td><td>114.23</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Calculated CNR values on the three ROIs in the two cases of (a) Data C1; and (b) Data C2.</figcaption></figure>

The protocol of this study (data collection and processing) was approved by the institutional ethical review board of the Mayo Clinic and was conducted in accordance with the principles of the 1964 Declaration of Helsinki. Two patients were involved in the experiments. All these patients have given their written informed consent to the participation. A non-conflict of interest for this work was declared. The proposed method was carried out in accordance with the approved guidelines. The data were analyzed fully anonymously. Specific Poisson noise was superimposed onto the raw projection data for each case in the library to synthesize a low dose level that corresponded to 1/4 the standard dose assuming 85mAs for tube current value³⁸. The detector has 736 × 64 elements with element size 1.2856 × 1.0947 mm². The source-to-axial distance and source-to-detector distance of the trajectory was respectively set to 59.5 cm and 108.56 cm. The helical trajectory of source was covered by 1152 view angles per cycle with a pitch equal to 0.6. For all the iterative algorithms, the whole projections were split into 32 subsets to accelerate computation.

**Visual Assessment.** Figure 10 and Fig. 11 display the selected axial views from the reconstructed results of Data C1 and Data C2. In Figs 10–11, images (a)–hepatic vein (see the zoomed region illustration(b) are respectively the illustrations for FDK and TV methods for the standard dose protocol, which the TV reconstructed images (Figs 10(b) and 11(b)) are used as reference; images (c)–(h) are respectively the illustrations for FDK, DFR-post, TV, PICCS<sub>DFR</sub>, GDSIR, DP-PICCS methods for the low dose protocol.

From the Figs 10 and 11, we can see that the DFR method works well with effective noise and artifacts suppression in images (b). It is found in Fig. 11(e,f) that the PICCS<sub>DFR</sub> approach leads to reconstruction in similar quality as the TV algorithm. Though demonstrating good performance in artifacts and noise suppression, the GDSIR method suffers from edge structure blurring of the hepatic vein (see the zoomed region illustration in Figs 10(g1) and 11(g1)). The illustrations in Figs 10 and 11 also show that, though still suffer less from noise-artifacts residual, the DP-PICCS reconstructed LDCT images present a better tiny structure identification than competing methods (see the red arrows in Figs 10 and 11).

Figures 12 and 13 display the selected sagittal and coronal views from the reconstructed CT volumes of Data C1 and Data C2. In Figs 12 and 13, it can be seen that the proposed DP-PICCS method performs better than the PICCS<sub>DFR</sub> and GDSIR methods, providing images with improved visual quality if using the TV reconstructed SDCT images (Figs 12(b) and 13(b)) as the references. Compared to other methods, the DP-PICCS method achieves a better preservation of anatomical features (see the zoomed region on liver tissue boundaries and arteries in Fig. 12(h1)) than the TV and PICCS<sub>DFR</sub> methods. From the Figs 12 and 13(g), we can see that the featured atoms in GDSIR method can be used to provide improved structure preservation in reconstruction. Also, it is found in Figs 12 and 13 that the proposed DP-PICCS method achieves the best performance in noise-artifacts
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Axial ROI illustrations of the DFR-post prior image when using dictionaries with different atom number and patch size. The rows from above to bottom correspond to the cases with patch sizes $B = \{6 \times 6 \times 3, 7 \times 7 \times 4, 8 \times 8 \times 5, 10 \times 10 \times 6\}$, respectively. The columns from left to right correspond to the cases with dictionaries in atom number (500, 500), (1000, 1000), (1500, 1500), (2000, 2000), respectively. The computation cost was tagged in unit second (s) in the right-bottom corner in each ROIs.</figcaption></figure>

suppression and anatomical features preservation among all competing reconstruction methods (see the structures pointed by red arrows in the Figs 12(h1) and 13(h1)).

**ISNR and CNR Based Quantification.** For clinical reconstruction data quantitative study, two metrics were utilized to give evaluation, namely, local Signal to Noise Ratio (ISNR) and Contrast-to-Noise Ratio (CNR). Here the ISNR and CNR are calculated by Eqs (17) and (18):

$$ \text{ISNR}(u) = \frac{1/N \sum_{n=1}^{N} u_n}{\sqrt{1/N \sum_{n=1}^{N} (u_n - 1/N \sum_{n=1}^{N} u_n)^2}} \quad (17) $$

$$ \text{CNR}(u, u_{BG}) = \frac{|\bar{u} - \bar{u}_{BG}|}{\sqrt{\sigma_u^2 - \sigma_{u_{BG}}^2}} \quad (18) $$

where $\bar{u}$ and $\bar{u}_{BG}$ are the mean intensities of the region of interest (ROI) and background region, respectively. $\sigma_u$ and $\sigma_{u_{BG}}$ are the associated standard deviations, respectively.

The ISNR values of the three different ROIs in Table 2, as indicated by the blue circles in Figs 10(a) and 11(a), were measured in the Data C1 and Data C2. The methods, i.e. DFR-post, TV, PICCSDFR, GDSIR, and DP-PICCS; yielded some gains on ISNR over the FDK reconstructions in the three different ROIs. It shows that, for the
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Evaluation of different parameter for torso phantom Case D2 data. (a) PSNR and UQI plots versus fidelity parameter $\beta$; (b) PSNR and UQI plots versus prior image parameter $\lambda$; (c) Axial ROI illustrations of reconstructed when using different $\beta$ and $\lambda$; (d) PSNR and UQI plots versus sparsity level $T_t$ and $T_r$. Here, the parameters were analyzed with the other parameters fixed to the values given in Table 1.</figcaption></figure>

low-dose case, the FDK method obtains the worst value, and the DP-PICCS method yields the highest value (in metric of ISNR for most ROIs) which are closer to the scores of TV SDCT reconstruction images than competing methods. We can also see that such quantitative results are consistent with the visual comparison in Figs 10-11.

Figure 14 lists the CNR values of the images reconstructed by the different cases for Data C1 and Data C2. The ROIs for the CNR measure are marked in Fig. 14(a) and (b) (red squares are ROIs and blue squares for backgrounds). It is also observed in Fig. 14 that the DFR-post processing significantly enhances the CNRs of LDCT images (most of them even higher than FDK reconstructed SDCT images) for Data C1 and Data C2, and the DP-PICCS algorithm leads to further higher CNRs. The plots in Fig. 14 also show that most of the CNRs of DP-PICCS indices are nearer to the TV reconstruction in SDCT case. This confirms the above visual comparison that the proposed method has good performance in terms of contrast preservation.

**Sensitivity of dictionary building parameters.** Performance of the proposed DP-PICCS is closely related to the parameters in dictionary construction, e.g. atom number C (dictionary size), and patch size B. Figure 15 displays the selected ROI (the LDCT images in Fig. 10(a) processed by DFR-post method) using different atom number C and patch size B. The computation cost (in unit second (s)) is tagged in the right-bottom corner in each ROI. The rows from above to bottom correspond to the results with $B = \{6 \times 6 \times 3, 7 \times 7 \times 4, 8 \times 8 \times 5, 10 \times 10 \times 6\}$, and the column from left to right respectively correspond to the results related to dictionaries with different atom numbers (500, 500), (1000, 1000), (1500, 1500), (2000, 2000) (in the format of $(C_t, C_r)$ for the dictionaries $D^t$ and $D^r$).

We can observe in Fig. 15 that an effective restoration of desirable attenuation information can be obtained if a (1000, 1000) sized dictionary with $8 \times 8 \times 5$ sized atoms are used in patch representation. It is also found in Fig. 15 that almost no visual difference can be discerned when the patch size and dictionary size increase over $8 \times 8 \times 5$ and (1000, 1000). Considering significantly larger computation cost is to be induced when enlarging the dictionary or atoms therein, in this study we use the dictionary size $C_t = C_r = 1000$ and $B = 8 \times 8 \times 5$ to maintain a good balance of representation accuracy and computation cost.

**Analysis of the parameters in the DP-PICCS reconstruction.** Several parameters need to be suitably set in the proposed DP-PICCS method, and these parameters include the fidelity parameter $\beta$, prior-image constraint parameter $\lambda$ and the sparsity parameters $T_t$ and $T_r$. In general, the fidelity parameter $\beta$ should be well set to balance the data fidelity and prior information constraint terms. The data fidelity term $\beta$ reflected the distance between the iterated images to the measured projection data, so should be lowered in the case of low dose noise projections$^{17,22}$. The parameter $\lambda$ for prior image term should be set according to the prior image quality. In feature patch sparse coding, the sparsity level parameters $T_t$ and $T_r$ are the maximum numbers of atom permitted
<table><caption>**Table 3.** Computation Cost Between The DFR-post TV PICCS<sub>DFR</sub>, GDSIR And DP-PICCS Methods (Unit: Second per Iteration).</caption><thead><tr><th>Method</th><th>Phantom data<br>(Size:512 × 512 × 30)</th><th>Clinical data<br>(Size:512 × 512 × 30)</th></tr></thead><tbody><tr><td>DFR-post</td><td>219 ± 2</td><td>225 ± 3</td></tr><tr><td>TV</td><td>3 ± 1</td><td>5 ± 1</td></tr><tr><td>PICCS<sub>DFR</sub></td><td>7 ± 1</td><td>11 ± 1</td></tr><tr><td>GDSIR</td><td>63 ± 2</td><td>77 ± 2</td></tr><tr><td>DP-PICCS</td><td>154 ± 2</td><td>167 ± 5</td></tr></tbody></table>

in representing each patch. Generally, larger $T_t$ and $T_r$ will give rise to a higher accuracy but increased computational cost. They should be set according to the measured projection noise level in reconstruction.

Figure 16(a) plots the PSNR and UQI values for different fidelity parameter $\beta$, using torso phantom Case D2 data with $\lambda = 0.35$ and $T_t = T_r = 8$. We can observe that the PSNR and UQI values are major affected by the fidelity parameter $\beta$ and attain the highest values when $\beta$ lies around 0.2. Figure 16(b) shows that prior image constraint parameter also have large impacts the reconstruction performance when $\lambda$ ranges from 0.31 to 0.43, and the highest PSNR and UQI value are reached when the $\lambda$ is 0.37. From the plots in Fig. 16(a,b) we can see that the performance of proposed method is quite sensitive to the fidelity parameter $\beta$ and prior image constraint parameter $\lambda$. The figures in Fig. 16(c) reflects that a large $\beta$ implies an increased fidelity weight of the measured noisy projections which results in increased noise and artifacts in the reconstruction, and a small $\beta$ leads to relatively increased weight of the regularization term which is related to increased smoothing effect. Figure 16(c) also shows that the reconstruction gets smoother when the value of $\lambda$ is increased. Figure 16(d) depicts the PSNR and UQI values for different sparsity level $T_t$ and $T_r$ with the other parameters set based on Table 1. The plots show that the proposed method is robust to sparsity level $T_t$ and $T_r$ from 8 to 12. The PSNR variation is less than 0.5 dB and UQI variation is less 0.01. The reason is due to fact that the same size dictionary is used and eight atoms to twelve atoms are enough to represent these volume patches. So, by taking into consideration both the reconstruction performance and the computation burden, we selected sparsity level $T_t = T_r = 8$ for the used dictionaries with 1000 atoms in this study.

We select parameters in the proposed approach following the assumption that the same parameters set can be used in reconstructing the data from the same type of CT scanner. Before large volume data reconstruction in our study, the feature dictionaries and parameters are trained to give satisfying results using a small amount of data in the same type.

**Computation cost.** Table 3. lists the computation overheads required per iteration (in seconds) for all reconstruction methods considered in our experiments. Both the GDSIR method and the proposed DP-PICCS method suffers from high computational operations to the sparse coding. In our experiments, it takes about 2148 s to train two 3D feature dictionaries (320 × 1000) from 1105142 sample patches for the proposed method. Specifically, we can notice that the GDSIR method and the proposed method require about 63 ± 2 seconds and 154 ± 3 seconds per iteration to reconstruct a 512 × 512 × 30 CT volume in torso phantom data experiment, respectively. In the clinical abdomen data experiment, one DP-PICCS iteration takes about 167 ± 5 seconds to reconstruct a 512 × 512 × 30 CT volume. We can see that the DP-PICCS method is about twice time consuming as much the GDSIR method due to the two 3D feature sparse coding for each volume patch. Besides, one of the biggest advantages of the TV method and PICCS<sub>DFR</sub> is its low computation complexity.

## Conclusion

The proposed DP-PICCS method is an extension of the PICCS method by using prior image knowledge from discriminative feature representation, which is reflected by a discriminative feature representation based constraint terms and high quality prior image. This discriminative feature representation is realized using a composite dictionary containing both desirable tissue attenuation features and undesirable noise-artifacts residual features.

In this study, three sparse constraint based methods were evaluated, i.e., TV, GDSIR and the proposed PICCS<sub>DFR</sub>. The TV and PICCS<sub>DFR</sub> methods can be regarded as voxel-wise updating methods while the GDSIR and DP-PICCS method are the patch-wise updating ones. Experiment results show that the PICCS<sub>DFR</sub> method cannot achieve satisfactory performance in suppressing noise-induced artifacts as the indistinguishability true tissue structures and image noise-artifacts are significant, and, the GDSIR method cannot also effectively deal with the tiny tissue structures in LDCT reconstruction. The proposed DP-PICCS method obtains PSNR values 1-2 dB higher than the PICCS<sub>DFR</sub> reconstruction, and the visual results also validate the improved ability of tissue structure discrimination. In the case that the projection data were seriously corrupted by excessive X-ray photon noise, the associative reconstructed images will suffer from serious noise-induced streak and mottle artifacts. In this case the DFR-post method cannot yield satisfactory results due to its ability in distinctive tissue and noise-artifacts structural.

The proposed DP-PICCS method can be feasibly implemented using the joint optimal reconstruction strategy in section II.B. However, due to the split alternating iteration scheme, the DP-PICCS method in sub-problem Eq. (6) and Eq. (7) are nonconvex for a global optimization. Similar to many existing dictionary learning based algorithms whose global convergence is an open issue, the proposed method also suffers from a lack of strict global convergence. Also, the DP-PICCS method needs to be accelerated to be more practically feasible. In the proposed
approach, parameters were selected following the assumption that the same parameter setting can be used in the reconstructing the data from the same type of CT scanner. In³⁹, a gradient projection method is applied in parameter selection where the regularization parameter is updated in an alternating mode. In⁴⁰, a reweighted objective function is defined to guide the selection of the regularization parameters. These approaches can be applied to guide the parameter selection of the proposed DP-PICCS reconstruction²²,³⁹,⁴⁰.

## References

1. Schubert, T. et al. CT-guided percutaneous biopsy of a mass lesion in the upper presacral space: A sacral transneuroforaminal approach. *Cardiovasc. Intervent. Radiol.* **35**, 1255–1257 (2012).

2. Siewerdsen, J. H. Cone-beam CT with a flat-panel detector: From image science to image-guided surgery. *Nucl. Instrum. Methods Phys. Res. A.* **648**, 241–250 (2011).

3. Zhao, B. et al. Evaluating variability in tumor measurements from same-day repeat CT scans of patients with non-small cell lung cancer. *Radiology.* **252**, 263–272 (2009).

4. Yu, L. et al. Radiation dose reduction in computed tomography: Techniques and future perspective. *Imag. Med.* **1**, 65–84 (2009).

5. McCollough, C. H., Bruesewitz, M. R. & Kofler, J. M. Jr. CT dose reduction and dose management tools: Overview of available options. *Radiographics.* **26**, 503–512 (2006).

6. Kalra, M. K. et al. Techniques and applications of automatic tube current modulation for CT. *Radiology.* **233**, 649–657 (2004).

7. Feldkamp, L., Davis, L. & Kress, J. Practical cone-beam algorithm. *JOSA A.* **1**, 612–619 (1984).

8. Donoho, D. L. Compressed sensing. *IEEE Trans. Inf. Theory.* **52**, 1289–1306 (2006).

9. Yu, H. & Wang, G. Compressed sensing based interior tomography. *Phys. Med. Biol.* **54**, 2791–2805 (2009).

10. Sidky, E., Kao, C. & Pan, X. Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT. *J X-Ray Sci Tech.* **14**, 119–139 (2006).

11. Sidky, E. & Pan, X. Image reconstruction in circular cone-beam computed tomography by constrained total-variation minimization. *Phys. Med. Biol.* **53**, 4777–4807 (2008).

12. Ritschl, L. et al. Improved total variation-based CT image reconstruction applied to clinical data. *Phys. Med. Biol.* **56**, 1545–1561 (2011).

13. Tian, Z. et al. Low-dose CT reconstruction via edge preserving total variation regularization. *Phys. Med. Biol.* **56**, 5949–5967 (2011).

14. Jia, X. et al. Gpu-based iterative cone-beam CT reconstruction using tight frame regularization. *Phys. Med. Biol.* **56**, 3787–3807 (2011).

15. Garduño, E., Herman, G. & Davidi, R. Reconstruction from a few projections by L1-minimization of the Haar transform. *Inverse problems.* **27**, 055006 (2011).

16. Xu, Q. et al. Low-dose X-ray CT reconstruction via dictionary learning. *IEEE Trans. Med. Imag.* **31**, 1682–1697 (2012).

17. Bai, T. et al. 3D dictionary learning based iterative cone beam CT reconstruction. *Int. J. Cancer Ther. Oncol.* **2**, 40–41 (2014).

18. Chen, Y. et al. Artifact suppressed dictionary learning for low dose CT image processing. *IEEE Trans. Med. Imaging.* **33**, 2271–2292 (2014).

19. Liu, J. et al. 3D Feature Constrained Reconstruction for Low Dose CT Imaging. *IEEE Trans. Circuits Syst. Video Technol.* (2016).

20. Chen, Y. et al. Discriminative feature representation: an effective post-processing solution to low dose CT imaging. *Phys. Med. Biol.* **62**, 2103–2132 (2017).

21. Caudevilla, O. & Brankov, J. G. Learning based prior for analyzer-based phase contrast image reconstruction”. *Biomedical Imaging (ISBI), 2015 IEEE 12th International Symposium on*. IEEE. **2015**, 1612–1615 (2015).

22. Mou, X. et al. Dictionary learning based low-dose x-ray CT reconstruction using a balancing principle. *Proc. SPIE 9212, Developments in X-Ray Tomography IX.* **9212**, 921207–921215 (2014).

23. Pfister, L. & Bresler, Y. Tomographic reconstruction with adaptive sparsifying transforms. *IEEE Int. Conf. Acoustics, Speech and Signal Processing(ICASSP)*, 6914–6918 (2014).

24. Chen, G. H., Tang, J. & Leng, S. Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets. *Med. Phys.* **35**, 660–663 (2008).

25. Chen, G. H. et al. Time-resolved interventional cardiac C-arm cone-beam CT: An application of the PICCS algorithm. *IEEE Trans. Med. Imag.* **31**, 907–923 (2012).

26. Chen, G. H., Tang, J. & Hsieh, J. Temporal resolution improvement using PICCS in MDCT cardiac imaging. *Med. Phys.* **36**, 2130–2135 (2009).

27. Lauzier, P. T. & Chen, G. H. Characterization of statistical prior image constrained compressed sensing (PICCS)—Part II. Application to dose reduction. *Med. Phys.* **40**, 021902 (2013).

28. Elbakri, I. & Fessler, J. Statistical image reconstruction for polyenergetic X-ray computed tomography. *IEEE Trans. Med. Imaging.* **21**, 89–99 (2002).

29. Elad, M. Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing. 121–350 (Springer, 2010).

30. Elad, M. & Aharon, M. Image denoising via sparse and redundant representations over learned dictionaries. *IEEE Trans. Image Process.* **15**, 3736–3745 (2006).

31. Rubinstein, R., Bruckstein, A. M. & Elad, M. Dictionaries for sparse representation modeling. *Proc. IEEE.* **98**, 1045–1057 (2010).

32. Candes, E. J., Eldar, Y. C., Needell, D. & Randall, P. Compressed sensing with coherent and redundant dictionaries. *Appl. Comput. Harmon. Anal.* **31**, 59–73 (2011).

33. Goldstein, T. & Osher, S. The split Bregman method for L1-regularized problems. *SIAM J. Imag. Sci.* **2**, 323–343 (2009).

34. Chen, S. S., Donoho, D. L. & Saunders, M. A. Saunders. Atomic decomposition by basis pursuit. *SIAM J. Sci. Comput.* **20**, 33–61 (1998).

35. Rubinstein, R., Zibulevsky, M. & Elad, M. Efficient implementation of the K-SVD algorithm using batch orthogonal matching pursuit. *CS Technion.* **40**, 1–15 (2008).

36. Wang, Z. & Bovik, A. C. A universal image quality index. *IEEE Signal Pro-cess Lett.* **9**, 81–84 (2002).

37. Ravishankar, S. & Bresler, Y. MR image reconstruction from highly undersampled k-space data by dictionary learning. *IEEE Trans. Med. Imaging.* **30**, 1028–1041 (2011).

38. McCollough, C. The 2016 NIH-AAPM-Mayo Clinic Low Dose CT Grand Challenge. http://www.aapm.org/GrandChallenge/LowDoseCT/# (2016).

39. Park, J. C. et al. Fast compressed sensing based CBCT reconstruction using Barzilai-Borwein formulation for application to on-line IGRT. *Med. Phys.* **39**, 1207–1217 (2012).

40. Zhang, C. et al. A Model of Regularization Parameter Determination in Low-Dose X-Ray CT Reconstruction Based on Dictionary Learning. *Computational and mathematical methods in medicine.* **2015**, 1–13 (2015).
**Acknowledgements**

The authors would also like to thank Dr. Cynthia McCollough (the Mayo Clinic, USA) for providing clinical projection data, which was arranged under the American Association of Physicists in Medicine, and grants EB017095 and EB017185 from the National Institute of Biomedical Imaging and Bioengineering.

**Author Contributions**

Y.C. and J.L. designed this project and wrote the code. Y.C, J.L. and L.X. prepared the data and performed all the experiments. Y.C. and J.L. wrote the main manuscript text. Y.H., H.S. and L.L. contributed to the method formation. L.Z., and G.C. checked and revised the whole paper language. All authors reviewed the manuscript.

**Additional Information**

**Competing Interests:** The authors declare that they have no competing interests.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_18.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2017