]Multi-Channel Convolutional Analysis Operator Learning for Dual-Energy CT Reconstruction

Alessandro Perelli^∗1, 2, Suxer Alfonso Garcia^∗1, Alexandre Bousse^1, Jean-Pierre Tasu^3, Nikolaos Efthimiadis^3, Dimitris Visvikis^1

$^1$ LaTIM, INSERM UMR 1101, Université de Bretagne Occidentale, 29238 Brest, France
$^2$ School of Science and Engineering, University of Dundee, Scotland DD1 4HN, UK
$^3$ Department of Radiology, University Hospital Poitiers, Poitiers, France
$^\ast$ These authors equally contributed to this work.

Objective. DECT has the potential to improve contrast, reduce artifacts and the ability to perform material decomposition in advanced imaging applications. The increased number or measurements results with a higher radiation dose and it is therefore essential to reduce either number of projections per energy or the source X-ray intensity, but this makes tomographic reconstruction more ill-posed.

Approach. We developed the MCAOL method to exploit common spatial features within attenuation images at different energies and we propose an optimization method which jointly reconstructs the attenuation images at low and high energies with a mixed norm regularization on the sparse features obtained by pre-trained convolutional filters through the CAOL algorithm.

Main results. Extensive experiments with simulated and real CT data were performed to validate the effectiveness of the proposed methods and we reported increased reconstruction accuracy compared to CAOL and iterative methods with single and joint TV regularization.

Significance. Qualitative and quantitative results on sparse-views and low-dose DECT demonstrate that the proposed MCAOL method outperforms both CAOL applied on each energy independently and several existing state-of-the-art MBIR techniques, thus paving the way for  dose reduction.

§ INTRODUCTION

DECT is an energy-resolved X-ray imaging technique that offers two sets of attenuation  measurements  utilizing  two  different  energy spectra acquired over the same anatomical region <cit.>. Since its invention, the technique has progressively been used for material decomposition and energy-selective imaging <cit.>.
Applications of DECT include differentiation and quantification of materials, tissue characterization, virtual monoenergetic imaging, automated bone removal in CT angiography, cardiovascular imaging, multiple contrast agent imaging and mapping of effective atomic number, and radiotherapy  <cit.>.

DECT acquisition techniques require two helical scans at two different tube voltages; therefore, two sets of projection data at different energy levels are collected and further reconstructed.
However, as the number of incident photons increases when irradiating with two sources the same anatomical region, the radiation dose increases proportionally <cit.>.

A reduction in radiation exposure as well as in acquisition time can be achieved by decreasing the number of projection angles and/or reducing the X-ray source dose.
However, aliasing artifacts can appear in the reconstructed images if the number of projection angles does not follow the Nyquist sampling theorem. Moreover, it is more challenging to achieve high-resolution, high-contrast image reconstruction due to the low SNR <cit.>.

In the literature, most of the development on low-dose CT reconstruction has focused on single image. Among the main techniques, MBIR methods are the most popular. These techniques exploit models of the imaging system's physics (forward models) along with statistical models of the measurements and noise and often simple object priors. They iteratively optimize model-based cost functions to estimate the underlying unknown image <cit.>. Typically, such cost functions consist of a data-fidelity term, e.g., least squares or NLL, capturing the imaging forward model and the measurement/noise statistical model and a regularizer term promoting smoothness, low-rank or sparsity <cit.>.
TV <cit.> has been proposed to solve incomplete projection data reconstruction problems and achieved good performance. However, TV assumes that the signal is piecewise constant, resulting in undesired patchy effects on the reconstructed images <cit.>.
Recent developments and interests within spectral CT refer to self-supervised methods using only noisy data <cit.> and dynamic DECT, which refers to randomly changing the energy threshold of the detectors to obtain the spectral information <cit.>.

Data-driven and learning-based approaches have gained much interest in recent years for biomedical image reconstruction. These methods learn representations of images and are used in combination with MBIR techniques to perform complex mappings between limited or corrupted measurements and high-quality images. Among those algorithms, data-driven sparse transforms such as DL <cit.> use a training dataset of high-resolution and denoised images to learn features, in an unsupervised manner, that can be used to reconstruct new images. These features take the form of “atoms”, which are regrouped into dictionaries and are used to sparsely represent the image <cit.>. DL-based image reconstruction integrates the learned atoms with the raw scanner data within a regularized MBIR context <cit.>. Other closely related methods include sparsifying transform learning <cit.> and the connection between data-adaptive models and convolutional deep learning algorithms <cit.> with an increase of interest in methods that leverage both learning-based and MBIR tools.

However, most DL methods are patch-based, and the learned features often contain shifted versions of the same features. The resulting learned dictionary may be redundant and therefore are memory demanding, which makes it difficult to utilize in 3D multimodal imaging.
To address these problems, CDL techniques utilize shift-invariant filters, providing a convenient and memory-efficient alternative to conventional DL techniques <cit.>. CDL approaches can be combined with MBIR by providing unsupervised prior knowledge of the target image.
The CDL approach can also be formulated from an analysis point of view <cit.> (sparse convolution) and is known as CAOL.

Despite the rapidly expanding research, the application of CDL to multi-channel images has received little attention <cit.>.
However, image reconstruction from DECT sparse-views or low-dose requires algorithms more advanced than the standard approach where attenuation at each measured energy is reconstructed independently. Notable models in the literature designed to promote structural similarity of images are JTV <cit.>, spectral patch-based penalty for the maximum-likelihood method <cit.>, tensor-based and coupled dictionary learning <cit.>, parallel level sets <cit.> and the prior rank, intensity and sparsity model (PRISM) <cit.>.

§.§ Main Contribution

We extend the CAOL approach to multichannel settings and we develop a MCAOL framework that can exploit direct joint reconstruction, given the low-dose DECT measurements,  where all the unknown images are reconstructed simultaneously by solving one combined optimization problem. This is the first time that MCAOL is applied to DECT image reconstruction and we demonstrate its superiority respect to CAOL. Furthermore, MCAOL requires considerably less memory compared to alternative DL approaches.
The joint reconstruction approach is developed for a low-dose data acquisition protocol which consists of collecting data using a sparse angular sampling, using a different X-ray energy in consecutive steps and low X-ray photon counts.

In DECT, a reasonable prior assumption is that attenuation images at different energies can be expected to be structurally similar in the sense that an edge (e.g., an organ boundary) that is present at one energy, is likely to be at same location and alignment with the other energies as well, even though the contrast between materials will be different at each energy.

MCAOL technique reconstructs attenuation images from the projection data combined with multi-channel filters trained on a dataset of reconstructed images.

The central idea of MCAOL is to learn unsupervised DECT muti-channel convolutional dictionaries that can provide a joint sparse representation of the underlined images by jointly learning filters for the different energies: each atom not only carries individual information for each energy individually but also inter-energy information.

By reconstructing DECT images using MBIR techniques in conjunction with MCAOL, the multi-energy information can be optimally used by allowing the images to “talk to each other” during the reconstruction process through the learned joint dictionaries, reducing noise while preserving image resolution. In order to deal with the extreme low-dose scenario, we model the Poisson model exactly and we solve the image optimization problem by using approximated quasi-Newton method with constrained memory to achieve accurate joint reconstruction with limited computational complexity.

§.§ Notation and Paper Organization

The remainder of this paper is structured as follows: Section <ref> briefly reviews the single-channel CAOL algorithm as initially proposed in <cit.> while Section <ref> presents the proposed approach to extend CAOL to multi-channel problems.  In Section <ref> we first describe the physical model of X-ray DECT from the continuous to discrete domain and the Poisson noise model; furthermore, we illustrate the proposed method to minimize the MCAOL reconstruction objective function using the exact low-dose DECT model. Finally, in Section <ref> comprehensive results of MCAOL on a numerical XCAT phantom <cit.> and clinical acquisition of a patient's full body are shown together with a comparison of its performance with state-of-the-art algorithms for both the sparse-views and low-dose DECT reconstruction.

Throughout the paper we adopt the following notations: matrices or discrete operators and column vectors are written respectively in capital and normal boldface type, i.e. $\*A$ and $\*a$, to distinguish from scalars and continuous variables written in normal weight; $[\*a]_j$ denotes the $j$-th entry of $\*a$; `$(\cdot)\transp$' denotes the transposition; `$\circledast$' denotes the 2D convolution operator; given an image  $\*x\in\R^J$ and a filter ${\*d} \in\R^P$, $P<J$, $\*d \circledast \*x$ represents the 2D convolution obtained by zero padding and then truncating the boundary such that the output dimension is the same as the input dimension $J$; an image $\*x\in\R^J$ is represented by a column vector for algebraic operations and by a $\sqrt{J} \times \sqrt{J}$ matrix for 2D convolutions. Finally, we indicate the set of positive real numbers with the notation $]0, +\infty[$.

§ LEARNING CONVOLUTIONAL REGULARIZERS FOR IMAGE RECONSTRUCTION: CAOL

In this Section we review the foundation of CAOL for MBIR which is achieved by solving an optimization problem of the form

\begin{equation}\label{eq:PML}
\min_{\*x \in \R^J} L(\*x,\*y) + \beta R(\*x)
\end{equation}

where $\*x \in \R^J$ is the 2D
image to reconstruct, $\*y\in\R^I$ is the observed measurement, $L$ is a data-fidelity term that incorporates the measurement model—generally taking the form of a NLL function—and $R$ is a regularizer weighted by $\beta>0$; $I$ and $J$ are respectfully the dimension of the measurement (number of detectors) and dimension of the image (number of pixels). The minimization is carried out with the help of iterative algorithms such as modified EM for ET <cit.> or PWLS combined with SPS for CT <cit.>.

The regularizer $R$ is designed such that the reconstructed image $\hat{\*x}(\*y)$ has desired properties, such as smoothness and sparsity of the gradient. It can be also trained so that $\hat{\*x}(\*y)$ can be sparsely represented as a linear combination of basic elements, or atoms, regrouped in a dictionary.

We consider the CAOL approach <cit.> where the image is sparsely represented with convolutional kernels (filters). In the analysis model, the image is represented with “sparsifying” filters ${\*d}_k \in\R^P$ by the analysis operator $\-A_{\*D}: \*x \mapsto \{\*d_k \circledast \*x\}$, such that

\begin{equation}
\*d_k \circledast \*x \approx \*z_k, \quad  \forall k = 1, \dots, K .
\end{equation}

where $\*z_k\in \R^J$ is a sparse feature image vector of the same dimension as the image $\*x$. The filters $\*d_k \in\R^P$ are vectorized images of dimension $P \ll J$ that are regrouped in a dictionary $\*D = \{\*d_k\}\in\R^{P\times K}$.

Learning the dictionary $\*D$ from a dataset of training images $\{\*x_l\in\R^J: l = 1,\ldots ,L\}$ corresponds to finding a collection of filters $\*D^{\star} = \{\*d^{\star}_k\}$ obtained by the following non-convex optimization problem

\begin{equation}\label{eq:d}
\*D^{\star} = \argmin_{\*D \in C} \, \min_{\{\*z_{l,k}\}} F_{\+a} (\*D, \{\*z_{l,k}\})
\end{equation}

where $C$ is the constraints on $\*D = \{\*d_k\}$ and with the training analysis objective function $F_\+a$ defined as

\begin{equation}\label{eq:cost_tr_CAOL}
F_{\+a}(\*D, \{\*z_{l,k}\}) = \sum_{l=1}^{L} \sum_{k=1}^{K} \frac{1}{2} \left\| \*d_k \circledast \*x_l  - \*z_{l,k}\right\|_2^2 + \alpha\left\| \*z_{l,k} \right\|_0
\end{equation}

where $\*z_{l,k}\in \R^J$ is the feature image associated to the training image $\*x_l$ and the filter $\*d_k$, $\left\|\cdot\right\|_0$ is the sparsity-promoting $l_0$ semi-norm defined for all $\*z=[z_1,\dots,z_J]\transp \in \R^J$ as

\begin{equation}
\|\*z\|_0 = \sum_{j=1}^J \bm{1}_{]0,+\infty[}(|z_j|)
\end{equation}

where $\bm{1}_A\colon\mathbb{R}\rightarrow\{0,1\}$ denotes the indicator function of a set $A\subset \mathbb{R}$, which is defined as $\*1_A(\xi)=1$ if $\xi\in A$ and $\*1_A(\xi)=0$ if $\xi\notin A$, and $\alpha > 0$ is a weight balancing between accuracy and sparsity. In Chun2019 the filters are enforced to satisfy the tight-frame conditions, i.e.,

\begin{equation}
C = \left\{ \{\*d_k\}: [\*d_1, \ldots , \*d_K] [\*d_1, \ldots , \*d_K]\transp = \frac{1}{P}\*I_K\right\}
\end{equation}

where $\*I_K$ is the $K\times K$ identity matrix, to promote filters diversity. The optimization problem eq:d is solved by alternating between minimization in $\*z$  and in $\*D$. The joint minimization in $\*D$, $\*z$ can be achieved with a BPEGM algorithm, which can be implemented using the CONVOLT <cit.>. BPEG-M alternates between minimization in $\*z$ and $\*D$. The minimization in $\*z$ is obtained with a hard-thresholding operator $\-T\colon\R^J\times\R_+^*\to\R^J$ defined at each row $j$ as

\begin{equation}\label{eq:thres}
[\-T(\*a,\beta)]_j = \cases{a_j  &for $\frac{1}{2} a_j^2  \ge \beta $\\
0 &otherwise\\}
\end{equation}

for all $\*a =  [a_1,\dots,a_J]\transp \in \R^J$ and for all $\beta >0$, which provides a global minimizer for $\*z\mapsto \frac{1}{2} \| \*a - \*z  \|_2^2  +   \beta \|\*z\|_0$, in such a way that

\begin{equation}
\mathcal{T}(\*d_k \circledast \*x_l , \alpha) = \argmin_{\*z_{l,k}} \frac{1}{2} \| \*d_k \circledast \*x_l  - \*z_{l,k} \|_2^2 + \alpha\| \*z_{l,k} \|_0 \, .
\end{equation}

The minimization in $\*D$ is achieved with a PMOC, which can be implemented with a singular value decomposition <cit.>. Finally the regularizer $R$ in the minimization problem eq:PML is derived from the learned filters $\*D^\star$ as

\begin{equation}\label{eq:pen_CAOL}
R(\*x) = \min_{\{\*z_k\}}\, \sum_{k=1}^{K} \frac{1}{2} \left\| \*d_k^\star \circledast \*x  - \*z_k\right\|_2^2 + \alpha\left\| \*z_k \right\|_0  \, .
\end{equation}

max totalsize=,center

[text_n, inner sep=0pt, node distance=0cm,  text width=12em] at (0,0) (tr_text) Training DECT dataset $\*x_{e,l} : e = 1, 2 ; l = 1, \ldots , L$;

[fig_n, above= 0.1 of  tr_text, inner sep=0pt] (tr_im1)

[fig_n, below= 0.1 of  tr_text, inner sep=0pt] (tr_im2)

[punkt,right= of  tr_text, label=above:MCAOL from $\{\*x_{e,l}\}$] (Tr_module)  ;

[->] (tr_text.east) edge node  (Tr_module.west);

[text_n, right= 0cm and 1.5cm of  Tr_module] (filt_text) Trained filters
$\*D^* = \{\*D^*_1, \*D^*_2\}$;

[->] (Tr_module.east) edge node  (filt_text.west);

[fig_n, below right= 0cm and -.4cm of  filt_text, inner sep=0pt] (acq_CT) 	;

[text_upn, below= of acq_CT, inner sep=0pt, node distance=0cm] (acq text) DECT acquisition;

[fig_n, above= 0.1cm and 0cmof  filt_text, inner sep=0pt] (filt_im) 	;

[punkt,right= 0cm and 7cm of filt_text] (Rec_module)  ;
[->] (filt_text.east) edge node  (Rec_module.west);

[text_upn, above= of  Rec_module] (recn_text) Reconstruction: apply $\*D^*$
to MBIR DECT;

[fig_n, below right= 0.1cm and -4.5cm of  Rec_module, inner sep=0pt] (sino_im) 	;

[text_upn, below= of  sino_im, inner sep=0pt, node distance=0cm] (sino text) Sinograms $\*y_e, e=1,2$;

[text_n, right= of  Rec_module] (rec_text) Reconstructed images
$\hat{\*x}_e(\*y_e), e = 1,2$;
[->] (Rec_module.east) edge node  (rec_text.west);

Diagram of MCAOL consisting of the unsupervised filter learning phase and the model-based iterative DECT reconstruction module.

§ MULTI-CHANNEL CONVOLUTIONAL ANALYSIS OPERATOR LEARNING

MBIR can be generalized to multi-channel imaging. Assuming we wish to reconstruct two images  $\*x_1,\*x_2\in\R^J$ of the same “object”  from two independent measurements $\*y_1\in\R^{I_1}$ and $\*y_2\in\R^{I_2}$ corresponding to two modalities, multi-channel MBIR can be achieved by using an iterative algorithm to solve

\begin{equation}\label{eq:MPML}
\min_{\*x_1,\*x_2\in\R^J} \rho_1 L_1(\*x_1,\*y_1) + \rho_2 L_2(\*x_2,\*y_2) + R_\mathrm{mc}(\*x_1,\*x_2)
\end{equation}

where $L_1$ and $L_2$ are the data-fidelity terms for $\*x_1$ and $\*x_2$, $R_\mathrm{mc}$ is a multi-channel regularizer and $\rho_1,\rho_2>0$ are weights. $R_\mathrm{mc}$ is designed to exploit the inference between the 2 channels $\*x_1$ and $\*x_2$, for example to promote structural similarities as proposed in ehrhardt2014.

MCAOL is a generalization of CAOL where the training is performed jointly on a set of images obtained from imaging modalities as depicted in Fig. <ref> for DECT. Let $\{(\*x_{1,l}, \*x_{2,l}) \in\R^J \times\R^J : l= 1,\ldots ,L\}$ be a training dataset consisting of $L$ pairs of images.
MCAOL learns the sparsifying filter pairs

\begin{equation}\label{eq:dual_d_def}
(\*d_{1,k}, \*d_{2,k}) \in\R^P \times\R^P : k= 1,\ldots ,K
\end{equation}

together with the extracted feature pairs

\begin{equation}\label{eq:dual_z}
(\*z_{1,l,k}, \*z_{2,l,k}) \in\R^J \times\R^J : k =1,\ldots , K,\;  l = 1,\ldots ,L \, .
\end{equation}

MCAOL is achieved by solving the following optimization problem, given the training image set $(\*x_{1,l}, \*x_{2,l})$

\begin{equation}\label{eq:dual_d_opt}
( \*D_1^{\star}, \*D_2^{\star} )
= \argmin_{\*D_1, \*D_2, \in C} \min_{\{\*z_{1,l,k}\} \atop \{\*z_{2,l,k}\}} F_{\+{mc}} \left(\*D_1, \*D_2, \{\*z_{1,l,k}\}, \{\*z_{2,l,k}\}\right)
\end{equation}

with the training cost function $F_{\+{mc}}$ defined as

\begin{eqnarray}
F_{\+{mc}} (\*D_1, \*D_2, \{\*z_{1,l,k}\}, \{\*z_{2,l,k}\})   = \sum_{l=1}^{L} \sum_{k=1}^{K} \frac{\gamma_1}{2} \left\| \*d_{1,k} \circledast \*x_{1,l}  - \*z_{1,l,k}\right\|_2^2 \nonumber  \\
+ \frac{\gamma_2}{2} \left\| \*d_{2,k} \circledast \*x_{2,l}  - \*z_{2,l,k}\right\|^2 + \left\| (\*z_{1,l,k}, \*z_{2,l,k}) \right\|_{1,0} \label{eq:Fma}
\end{eqnarray}

where $\gamma_1, \gamma_2 > 0$ are weights and the semi-norm $\|\cdot \|_{1,0}$ on $\R^J \times \R^J$ is defined for all $\*z_1 = [z_{1,1},\dots,z_{1,J}]\transp\in\R^J$ and for all $\*z_2 = [z_{2,1},\dots,z_{2,J}]\transp\in\R^J$ as

\begin{equation}\label{eq:l_1_0}
\| (\*z_1, \*z_2) \|_{1,0} = \sum_{j=1}^J\*1_{]0, +\infty[}\left(|z_{1,j}| + |z_{2,j}|\right)
\end{equation}

$\|\cdot \|_{1,0}$ denotes the $l_{1,0}$ norm. It promotes joint sparsity, i.e., with zero and non-zero values at the same locations, of image features in all the modalities, that are encoded by the multi-channel dictionary $\*D_1, \*D_2$.

To solve (<ref>) we utilize an iterative alternating minimization scheme where the update of the filters $( \*D_1, \*D_2 )$ is solving a block non-convex problem using the PMOC algorithm <cit.> while for the update of the sparse codes $(\*z_{1,l,k}, \*z_{2,l,k})$ we deploy a multi-channel hard-thresholding operator  $\-T_{\+{mc}}\colon\R^J\times \R^J\times (\R_+^*)^2 \to \R^J\times \R^J$ defined at each row $j$ as

\begin{equation}\label{eq:mthres}
[\-T_{\+{mc}}(\*a_1,\*a_2,\*\gamma)]_j = \cases{(a_{1,j},a_{2,j})  & for $\frac{1}{2}\gamma_1 a_{1,j}^2 + \frac{1}{2}\gamma_2 a_{2,j}^2 \ge 1$\\
(0,0) &otherwise\\}
\end{equation}

for all $\*a_1  =  [a_{1,1},\dots,a_{1,J}]\transp \in \R^J$, $\*a_2  =  [a_{2,1},\dots,a_{2,J}]\transp \in \R^J$  and  for all $\*\gamma = (\gamma_1,\gamma_2)\in(\R_+^*)^2$, which provides a global minimizer for $(\*z_1,\*z_2) \mapsto \frac{\gamma_1}{2} \| \*a_1 - \*z_1 \|_2^2  + \frac{\gamma_2}{2} \| \*a_2 - \*z_2 \|_2^2 + \|\*z_1,\*z_2\|_{1,0} $ <cit.>, in such a way that

\begin{eqnarray}
\-T_{\+{mc}}(\*d_{1,k} \circledast \*x_{1,l}, \*d_{2,k} \circledast \*x_{2,l}, \*\gamma  ) =  \argmin_{\*z_{1,l,k},\*z_{2,l,k}}  \, \Bigg\{ \frac{\gamma_1}{2} \| \*d_{1,k} \circledast \*x_{1,l}  - \*z_{1,l,k}\|_2^2 \nonumber  \\
+ \frac{\gamma_2}{2} \| \*d_{2,k} \circledast \*x_{2,l}  - \*z_{2,l,k}\|^2 + \| (\*z_{1,l,k}, \*z_{2,l,k}) \|_{1,0}   \Bigg\}   \, .
\end{eqnarray}

Finally the regularizer $R_{\+{mc}}$ in the minimization problem eq:MPML is derived from the learned filters $(\*D_1^\star,\*D_2^\star)$ as

\begin{eqnarray}
R_{\+{mc}}(\*x_1,\*x_2)  =   \min_{\{\*z_{1,k}\} \atop \{\*z_{2,k}\}} \sum_{k=1}^{K} \frac{\gamma_1}{2} \left\| \*d^\star_{1,k} \circledast \*x_1  - \*z_{1,k}\right\|_2^2 \nonumber \\
\quad + \frac{\gamma_2}{2} \left\| \*d^\star_{2,k} \circledast \*x_{2,l}  - \*z_{2,k}\right\|_2^2 + \left\| (\*z_{1,k}, \*z_{2,k}) \right\|_{1,0} \label{eq:pen_MCAOL}
\end{eqnarray}

The pseudo-code for MCAOL training procedure is summarized in Algorithm <ref>.

DE Training Dataset $\*x_{e,l}$, $l=1,\ldots ,L$, $e = 1,2$, joint sparsity weights $\bm\gamma = (\gamma_1,\gamma_2)$, #outer iterations $N_\mathrm{outer}$
Learned filters $( \*D_1^{\star}, \*D_2^{\star} )$
$(\*D_1^0, \*D_2^0) \leftarrow$ $\+{Normalized}$ $\+{random}$ $\+{initialization}$

$t = 0, \ldots, N_\mathrm{outer}-1$
Update sparse codes (in parallel)
$(\*z^{t+1}_{1,l,k}, \*z^{t+1}_{2,l,k}) \leftarrow$ $\mathcal{T}_{\+{mc}}(\*d^{t+1}_{1,k} \circledast \*x_{1,l}, \*d^{t+1}_{2,k} \circledast \*x_{2,l}, \bm\gamma)$

Update Filters
$ \*D^{t+1}_{1} \leftarrow \mathrm{PMOC}(\*x_{1,l},\*z^{t+1}_{1,k})$
$ \*D^{t+1}_{2} \leftarrow \mathrm{PMOC}(\*x_{2,l},\*z^{t+1}_{2,k})$

$\*D_1^{\star} \leftarrow  \*D_1^{N_\mathrm{outer}} $
$\*D_2^{\star} \leftarrow  \*D_2^{N_\mathrm{outer}} $
MCAOL Training Algorithm

§ DUAL-ENERGY CT RECONSTRUCTION WITH MULTI-CHANNEL CAOL

§.§ X-ray CT Discrete Model

In this section, we describe the CT discrete physical measurement process with the spectrum of the X-ray source beams composed of two different energies. We consider the case of 2D slice-by-slice imaging systems.

For image reconstruction we assume that the continuous attenuation image $\mu_e(\*r)$ which denotes the linear attenuation coefficient at position $\*r\in \R^2$ and the energy level $e=1,2$, can be represented by a linear combination of basis functions $\{b_j\}$  associated to a discrete sampling on a $\sqrt{J} \times \sqrt{J}$ Cartesian grid,

\begin{equation}\label{eq:paramdiscrete}
\mu_e(\*r) = \sum_{j=1}^{J} x_{e,j} b_j(\*r) \, ,
\end{equation}

where $x_{e,j} > 0$ for all $j=1,\dots,J$ and all $e=1,2$.
The line integral becomes a summation:

\begin{equation}\label{eq:linCT}
\int_{\R}\mu_e(\*\nu_i(l)) \, \ud l = \sum_{j=1}^{J} x_{e,j} \int_{\R}b_j(\*\nu_i(l))\, \ud l = \sum_{j=1}^{J} a_{i,j} x_{e,j}   % L_i
\end{equation}

where $\*\nu_i(l)  = \*s_i + l  \vec{\*\epsilon}_i \in \R^2$ is a parametrisation of the $i$-th ray emitted from the source $\*s_i$ with direction $\vec{\*\epsilon}_i$, $a_{i,j} \triangleq \int_{\R}b_j(\*\nu_i(l))\,\ud l$ is the contribution of the $j$-th pixel to the $i$-th ray.
The system matrix $\*A$ is constructed as an under-determined matrix of dimensions $I\times J$ where $I = N_\mathrm{d}\times N_{\theta}$ with $N_\mathrm{d}$ and $N_{\theta}$ being respectively the number of detectors and $N_{\theta}$ and the number of angles (projections), and is defined as $[\*A]_{i,j} = a_{i,j}$, $\forall\; i=1,\dots,I,\;\forall\; j=1,\dots,J$.
The spectral X-ray mathematical discrete model is based on the Beer's law which provides the X-ray intensity after transmission. The expected number of detected photons $\bar{y}_{i,e}$ is then redefined as a function of the discrete image $\*x_e$ as

\begin{equation}\label{eq:mean_Poisson_DECT}
\bar{y}_{i,e} (\*x_e) = S_e \+e^{ - [\*A \*x_e]_i } + \eta_{e,i}
\end{equation}

where $\*x_e = [x_{e,1},\ldots, x_{e,J}]\transp\in\R^J$ is the vector of attenuation coefficients at source energy $e$, $S_e$ is the mean photons flux at the $e$-th energy bin, as we assume a mono-energetic intensity, and $\eta_{e,i}\in\R^+$ is a known additive term representing the expected number of background events (primarily from
scatter). In the case of normal exposure, the number of detected photons follows a Poisson distribution, i.e.,

\begin{equation}\label{eq:poisson_measurment}
y_{i,e} \sim \+{Poisson}(\bar{y}_{i,e} (\*x_e))
\end{equation}

and the measurements at each energy bin $e=1,2$ are stored in a vector $\*y_e = [y_{e,1},\dots,y_{e,N_\mathrm{d} \cdot N_{\theta}}]\transp$.

Although monochromatic X-ray source does not usually hold for scanners in clinical practice, a common effective strategy consists of applying a polychromatic-to-monochromatic source correction pre-processing step <cit.>, and in the rest of the paper we will therefore assume that we have a monoenergetic source or that it has already been appropriately corrected.

§.§ Low-Dose CT Reconstruction

In case of low X-ray dose, since the photons counts can be very limited, the Gaussian approximation is no longer applicable as the  logarithm of the data cannot be computed. We therefore chose to perform sparse view CT reconstruction from the raw measurements $(\*y_1,\*y_2)$ by solving the minimization problem eq:MPML, with positivity constraints on $(\*x_1,\*x_2)$,  using the Poisson NLL functions $L_1$ and $L_2$  defined as

\begin{equation}\label{eq:NLL}
-L_e(\*x_e,\*y_e) = \sum_{i=1}^{I} y_{e,i}\log \bar{y}_{i,e} (\*x_e) - \bar{y}_{i,e} (\*x_e) , \quad e=1,2
\end{equation}

and the trained regularizer $R_\+{mc}$ derived from the learned filters $(\*D_1^\star,\*D_2^\star)$ as in eq:pen_MCAOL.
Therefore, substituting eq:NLL and eq:pen_MCAOL into the minimization eq:MPML,
we obtain the following explicit expression for the MCAOL DECT reconstruction problem:

\begin{eqnarray}\label{eq:NLL_MCAOLPoisson}
(\*x^{\star}_1,\*x^{\star}_2) = \argmin_{\*x_1,\*x_2 \geq 0} \, \sum_{e=1}^2\rho_e\underbrace{ \sum_{i=1}^{I} y_{e,i}\log \bar{y}_{i,e} (\*x_e) - \bar{y}_{i,e} (\*x_e)}_{ L_e(\*x_e,\*y_e)}  \nonumber \\
+ \underbrace{\min_{\{\*z_{1,k}\} \atop \{\*z_{2,k}\}} \sum_{k=1}^{K}  \Bigg\{\sum_{e=1}^2 \frac{\gamma_e}{2} \left\| \*d^{\star}_{e,k} \circledast \*x_{e}  - \*z_{e,k}\right\|_2^2 \Bigg\} + \left\| (\*z_{1,k}, \*z_{2,k}) \right\|_{1,0}}_{R_\mathrm{mc}(\*x_1,\*x_2)}
\end{eqnarray}

We solve the minimization problem (<ref>) by the alternating estimation of the sparse feature images and the linear attenuation images $\{\*x_e: e = 1,2\}$. Given the current estimates of the sparse coefficients $\{\*z^t_k: k=1,\ldots, K\}$, the image update $\*x_e^t$ at iteration $t$ is obtained through the following  minimization problem

\begin{eqnarray}
\*x_e^{t}  =   \argmin_{\*x_e \in (\R^+)^J}  \Phi^t_e(\*x_e) \label{eq:MCAOLoptimization} \\
\mathrm{with} \;\; \Phi_e(\*x_e) =   \rho_e L_e( \*x_e,\*y_e  )  + \frac{\gamma_e}{2} \sum_{k=1}^{K}  \left\| \*d^{\star}_{e,k} \circledast \*x_e  - \*z^{t}_{e,k}\right\|_2^2   \, .\nonumber
\end{eqnarray}

In this work, we utilized a LBFGS algorithm <cit.> to solve eq:MCAOLoptimization. The LBFGS iterative solver estimates $\*x_e^{t}$ starting the previous iterate $\*x_e^{t-1}$. We define the first estimate as $\*x_e^{t,(0)} = \*x_e^{t}$. At inner iteration $q$, given a current estimate $\*x_e^{t,(q)}$, the new estimate $\*x_e^{t,(q+1)}$ is obtained as

\begin{eqnarray}\label{eq:LBFGS}
\*x_e^{t,(q+1)} & = & \*x_e^{t,(q)} - s^{\star}\*B^{t,(q)}\nabla\Phi_e^t(\*x_e^{t,(q)}) \nonumber \\
\mathrm{with}\quad s^{\star} & = & \argmax_{s\in[0,1]}\chi(s) \label{eq:line_search} \\
\mathrm{and}\quad \chi(s) & = & \Phi_e^t\left(\*x_e^{t,(q)} - s\*B^{t,(q)}\nabla\Phi_e^t(\*x_e^{t,(q)})\right) \nonumber
\end{eqnarray}

where $\*B^{t,(q)}$ is an approximate inverse Hessian of $\Phi_e^t$ evaluated at $\*x_e^{t,(q)}$. The matrix/vector product $\*B^{t,(q)}\nabla\Phi_e^t(\*x_e^{t,(q)})$ in (<ref>) is directly computed (without storing $\*B^{t,(q)}$) from the $m$ previous iterates $\*x_e^{t,(q-p)}$, $p = 0, \ldots, m-1$. An approximate solution of the line-search sub-problem is obtained by backtracking to match the Wolfe Conditions. The iterative scheme (i.e., w.r.t. $q$) is repeated until either a convergence criterion is met. A more detailed explanation can be found in bousse2020hypo. We utilized the implementation proposed in zhu1997algorithm. We also used the LBFGS algorithm to minimize $L_1(\cdot,\*y_1)$ and $L_2(\cdot,\*y_2)$ (without penalty) in order to obtain initial images $\*x_1^0$ and $\*x_2^0$.

The other part of the alternating scheme is to update the sparse features update $\*z^t_{e,k}$ given the current estimate of $\*x^t_e$ the sparse features update $\*z^t_{e,k}$. This step is achieved using the multi-channel thresholding operator defined in eq:mthres.

The pseudo-code for MCAOL reconstruction algorithm is detailed in Algorithm <ref>.

Initial images $(\*x_1^0, \*x_2^0)$, DECT learned filters $\*D^{\star} = ( \*D_1^{\star}, \*D_2^{\star} )$, joint sparsity weight $\bm\gamma = (\gamma_1,\gamma_2)$, penalty weights $\bm\rho = (\rho_1,\rho_2)$, DE sinogram $\*y = (\*y_1,\*y_2)$, system matrix $\*A$, intensities $(S_1,S_2)$ , #outer iterations $N_\+{outer}$.
Reconstructed images $( \*x^{\star}_1, \*x^{\star}_2 )$
$t = 0, \ldots, N_\+{outer}-1$
Update sparse codes (in parallel)
$(\*z^{t+1}_{1,k}, \*z^{t+1}_{2,k}) \leftarrow$ $\mathcal{T}_{\+{mc}}(\*d^\star_{1,k} \circledast \*x_1^t, \*d^\star_{2,k} \circledast \*x_2^t, \bm\gamma)$

Update linear attenuation images
$\*x^{t+1}_{1} \leftarrow \mathrm{L}\mbox{-}\mathrm{BFGS} (\Phi_1^t , \mathrm{init}= \*x_1^{t} \mid \*y_1, \*z^{t+1}_{1}, \*D_1^{\star}, \*A,S_1,\rho_1,\gamma_1 )  $
$\*x^{t+1}_{2} \leftarrow \mathrm{L}\mbox{-}\mathrm{BFGS} (\Phi_2^t , \mathrm{init}= \*x_2^{t} \mid \*y_2, \*z^{t+1}_{2}, \*D_1^{\star}, \*A,S_2,\rho_2, ,\gamma_2)  $

$\*x^\star_{1} \leftarrow \*x^{N_\+{outer}}_{1} $
$\*x^\star_{2} \leftarrow \*x^{N_\+{outer}}_{2} $
MCAOL Reconstruction Algorithm

§ VALIDATION

We validated the proposed methods on two different DECT low-dose acquisition setup. In particular, we analyzed  the case of sparse-view DECT reconstruction with normal photon dose and the case of extreme low-photon counts with increased number of views. By approximating the dose as the product of the number of views and photon counts, the latter case represents a more challenging scenario since the overall dose considered is lower than the sparse-view case. Our implementation was based on CONVOLT <cit.>.

§.§ Methods Used for Comparison

The objective of the simulations with sparse views and normal X-ray source intensity is to demonstrate that MCAOL achieves improved accuracy compared reconstructing each energy separately by solving eq:PML with the CAOL regularizer defined in eq:pen_CAOL  and with the edge-preserving TV regularizers, as well as simultaneously by solving eq:MPML with the JTV regularizer, respectfully defined as

\begin{equation}\label{eq:tv}
R_\+{tv}  (\*x)  = \sum_{j=1}^J \sum_{k\in\-N_j}  \omega_{j,k} \sqrt{ (x_j - x_k)^2  +  \varepsilon }
\end{equation}

\begin{equation}	\label{eq:jtv}
R_\+{jtv}(\*x_1,\*x_2)   =
\sum_{j=1}^J \sum_{k\in\-N_j} \omega_{j,k}  \sqrt{(x_{1,j} - x_{1,k})^2  +  (x_{2,j} - x_{2,k})^2  + \varepsilon }
\end{equation}

where $\-N_j$ denotes the 8 nearest neighboring pixels of pixel $j$ and $\omega_{j,k}$ are weights ($\omega_{j,k}=1$ for axial neighbors and $\omega_{j,k}=1/\sqrt{2}$ for diagonal neighbors), and $\varepsilon>0$ is a small real value to ensure differentiability. For each method, we used the LBFGS solver to estimate $\*x_1$ and $\*x_2$.

The experiment with extreme low-counts aims at demonstrating that considering a weighted least-squares approximation of the log-likelihood function no longer guarantees effective reconstruction results, instead the exact Poisson statistics should be accounted. This results in a degradation of the performance of CAOL when optimized through the PWLS solver while using the quasi-Newton solver LBFGS leads to improved qualitative and quantitative results.

§.§ Methodology

All experiments were validated by generating the DECT measurements as in Eq. eq:mean_Poisson_DECT and then running $n=20$ Poisson noise instances as in Eq. eq:poisson_measurment from a GT image $\*x_e^\star = [x_{e,1}^*,\dots,x_{e,J}^*]\transp\in\R^J$, $e=1,2$. As performance metrics, we considered the mean absolute bias ($\+{AbsBias}$) error function defined as follows

\begin{equation}
\+{AbsBias} = \frac{1}{N_{\-R}}\frac{1}{N_{\+{noise}}} \sum_{j\in\-R} \sum_{n=1}^{N_{\mathrm{noise}}}\left| x_{e,j}^{[n]} - x_{e,j}^*  \right|
\end{equation}

where $x_{e,j}^{[n]}$ indicates the reconstructed linear attenuation coefficient at image pixel $j$ from the $n$-th Poisson noise replicate, $\-R$ is the spatial region of interest and $N_{\-R}$ is the number of pixels in the region $\-R$. Furthermore, we compute the STD defined as

\begin{equation}
\+{STD} = \frac{1}{N_{\-R}} \sum_{j\in\-R} \sqrt{\frac{1}{N_{\+{noise}}} \sum_{n=1}^{N_{\mathrm{noise}}}\left( x_{e,j}^{[n]} - \bar{x}_{e,j}  \right)^2}
\end{equation}

$\bar{x}_{e,j} = \frac{1}{N_{\+{noise}}}\sum_{n=1}^{N_{\mathrm{noise}}} x_{e,j}^{[n]}$. In this work $\-R$ corresponds to the nonnegative pixels region of $\*x_e^\star$ and is the same for both energy levels.

The simulations were repeated for all the methods, for different values of the regularization parameters in the objective functions eq:PML and eq:MPML in order to plot $\+{AbsBias}$ / $\+{STD}$ curves (Fig. <ref>, <ref> and <ref>). The quality of the reconstruction is assessed by the proximity of the curve to the origin. Training and reconstruction were performed according to the below-described settings.

The optimization problem (<ref>) is minimized using the BPEGM algorithm <cit.> with normalized input dataset.  To investigate the trade-off between accuracy and features sparsity, we tested (<ref>) with different values of $\gamma_1 = \gamma_2$ with filter $(\*d_{1,k}, \*d_{2,k})$ of dimension $P = 49$ and number of filters $K=49$. For each simulation, we tuned $\gamma_1 = \gamma_2$ by testing different values to investigate the effect. For all the datasets, we have used a training of $L=25$ images for each energy $e$. Regarding the BPEGM algorithm, we set the tolerance value equal to $10^{-4}$ and the maximum number of iterations to $3\times 10^3$. For the CAOL training algorithm, we have used the same settings as detailed for MCAOL except that we tuned a single regularization weight $\alpha$ in the optimization problem eq:cost_tr_CAOL for each separate energy channel.

Reconstruction MCAOL and JTV reconstructions were achieved by solving eq:MPML with $R_{\+{mc}}$ defined as eq:pen_MCAOL and eq:jtv respectively, while CAOL and TV reconstructions by solving eq:PML for each energy bin $e=1,2$ separately with $R$ defined as eq:pen_CAOL and eq:tv respectively. MCAOL and CAOL were achieved using $N_{\+{outer}}=300$ outer iterations while the inner image update is obtained using the LBFGS algorithm with $300$ iterations. The $(\gamma_1,\gamma_2)$-values and $\beta$-values were the same as for training. TV and JTV reconstructions were achieved with the LBFGS algorithm with $300$ iterations. The measurements were obtained from the GT images $\*x_e^\star$ outside of the training set and the reconstructions were repeated for each noise instance $n$, for a range of  $(\rho_1,\rho_2)$-values with $\rho_1 = \rho_2$ and for a range of $\beta$-values, in order to obtain $\+{AbsBias}$-versus-$\+{STD}$ curves.

We performed sparse-views and low-dose experiments on a simulated XCAT phantom and clinical data to assess the potential of the method for medical practice as detailed below.
The experiments were conducted with fixed X-ray dose amount, i.e., by selecting the number of angles and the X-ray source intensity, and we evaluated the quality of the linear attenuation  images reconstructed with different methods, both qualitatively and quantitatively.

§.§ Results on XCAT Phantom

For the unsupervised MCAOL and CAOL training, the numerical data consists of 1-mm pixel-width $512\times 512$ torso axial slice images generated from the XCAT phantom for $60$ keV and $120$ keV energies.

We utilized 20 slice pairs from the XCAT phantom, each pair consisting of a slice at $E_1=60$ KeV and a slice at $E_2=120$ KeV, to train the filters. An additional slice pair—not part of the training dataset—was used to generate the projection data  (i.e, following eq:mean_Poisson_DECT and eq:poisson_measurment) for reconstruction, as detailed below. We used the MCAOL weights parameters $\gamma_1 = \gamma_2 = 800$ and the CAOL parameter $\alpha = 0.01$.

Learned $l_{1,0}$ filters via MCAOL

($E_1$), $\gamma_1 = 800$     ($E_2$), $\gamma_2 = 800$

Learned $l_0$ filters via CAOL

($E_1$), $\alpha = 0.01$     ($E_2$), $\alpha = 0.01$

Learned filters $\{(\*d_{1,k}, \*d_{2,k})\}$ with $P=K=49$ using the XCAT training dataset, for subfig:MCAOL:filters:XCAT_l0_l10 MCAOL and subfig:MCAOL:filters:XCAT_l0_l0 CAOL.

$\bm{z}_{1,k}$, CAOL

$\bm{z}_{2,k}$, CAOL

$\bm{z}_{1,k}$, MCAOL

$\bm{z}_{2,k}$, MCAOL

XCAT Phantom: estimated sparse feature maps $\bm{z}_{e,k}$ for $e=1,2$ and $k=1,\dots,49$ using CAOL (subfig:z1_sep_xcat and subfig:z2_sep_xcat) and MCAOL (subfig:z1_joint_xcat and subfig:z2_joint_xcat); color scale: red for positive values, blue for negative values.

= 8.2pt

[spy using outlines=rectangle,yellow,magnification=1.25,size=12mm,connect spies]
on (0.3,-0.65) in node [left] at (-0.3,-1.2);

[spy using outlines=rectangle,yellow,magnification=1.25,size=12mm,connect spies]
on (0.3,-0.65) in node [left] at (-0.3,-1.2);

No prior

[spy using outlines=rectangle,yellow,magnification=1.25,size=12mm,connect spies]
on (0.3,-0.65) in node [left] at (-0.3,-1.2);

[spy using outlines=rectangle,yellow,magnification=1.25,size=12mm,connect spies]
on (0.3,-0.65) in node [left] at (-0.3,-1.2);

[spy using outlines=rectangle,yellow,magnification=1.25,size=12mm,connect spies]
on (0.3,-0.65) in node [left] at (-0.3,-1.2);

[spy using outlines=rectangle,yellow,magnification=1.25,size=12mm,connect spies]
on (0.3,-0.65) in node [left] at (-0.3,-1.2);

[spy using outlines=rectangle,yellow,magnification=1.25,size=12mm,connect spies]
on (0.3,-0.65) in node [left] at (-0.3,-1.2);

[spy using outlines=rectangle,yellow,magnification=1.25,size=12mm,connect spies]
on (0.3,-0.65) in node [left] at (-0.3,-1.2);

[spy using outlines=rectangle,yellow,magnification=1.25,size=12mm,connect spies]
on (0.3,-0.65) in node [left] at (-0.3,-1.2);

[spy using outlines=rectangle,yellow,magnification=1.25,size=12mm,connect spies]
on (0.3,-0.65) in node [left] at (-0.3,-1.2);

[spy using outlines=rectangle,yellow,magnification=1.25,size=12mm,connect spies]
on (0.3,-0.65) in node [left] at (-0.3,-1.2);

[spy using outlines=rectangle,yellow,magnification=1.25,size=12mm,connect spies]
on (0.3,-0.65) in node [left] at (-0.3,-1.2);

Comparison of reconstructed XCAT phantom from different reconstruction methods for sparse-view CT with top row corresponding to high energy $E_1 = 120$ keV and bottom row to low energy $E_2 = 60$ keV: sugfig:xcat_gt GT XCAT test image, sugfig:xcat_lbfgs minimization of the NLL function without prior, sugfig:xcat_caol CAOL reconstruction, sugfig:xcat_mcaol MCAOL joint reconstruction, sugfig:xcat_tv separate reconstruction using TV prior and sugfig:xcat_jtv joint reconstruction using JTV prior.

Fig. <ref> shows the pairs $(\*d_{1,k}, \*d_{2,k})$ of learned convolutional filters obtained by MCAOL (Fig. <ref>) and separate learning with CAOL (Fig. <ref>). From a qualitative point of view, we observe that MCAOL filter pairs $(\*d_{1}, \*d_{2})$ look similar as the edges are identical in the $2$ energy images,  as opposed to the CAOL filters.

$60$ keV

grid = major,
legend columns=2,
legend cell align=left,
legend entries=MCAOL, CAOL, TV, JTV,
legend style=at=(0.7,0.2),anchor=north
[color=blue,mark=*] table[x=xMCAOL, y=yMCAOL] Fig_5a.txt;
[color=black,mark=square*] table[x=xCAOL, y=yCAOL] Fig_5a.txt;
[color=red,mark=diamond*] table[x=xTV, y=yTV] Fig_5a.txt;
[color=green,mark=pentagon*] table[x=xJTV, y=yJTV] Fig_5a.txt;

$120$ keV

grid = major,
legend columns=2,
legend cell align=left,
legend entries=MCAOL, CAOL, TV, JTV,
legend style=at=(0.7,0.19),anchor=north
[color=blue,mark=*] table[x=xMCAOL, y=yMCAOL] Fig_5b.txt;
[color=black,mark=square*] table[x=xCAOL, y=yCAOL] Fig_5b.txt;
[color=red,mark=diamond*] table[x=xTV, y=yTV] Fig_5b.txt;
[color=green,mark=pentagon*] table[x=xJTV, y=yJTV] Fig_5b.txt;

Plot of the mean absolute bias ($\+{AbsBias}$) versus the standard deviation ($\+{STD}$) for the XCAT phantom at subfig:Bias_vs_std_xcat_e1 low X-ray source energy ($60$ keV) and subfig:Bias_vs_std_xcat_e2 high X-ray source energy ($120$ keV).

In order to generate the  sparse-view DECT projection measurements eq:mean_Poisson_DECT, we modeled the projector $\*A$ with a 2-mm FWHM resolution parallel beam system and we used 1-mm pixel-width $406\times 406$ GT torso axial-slice images with attenuation coefficients $\*x^\star_1, \*x^\star_2$ at energies $120$ keV (high) and $60$ keV (low) which differs from the training examples. The simulation consisted on generating sparse-view sinograms with $406$ detector pixels and $60$ regularly spaced projection angles, where $360^{\circ}$ is the full view rotation. A monochromatic source with $\bar{S}_e = 10^5$ incident photons and $100$ background events was used to generate each sinogram.

To support the statement that joint sparsity of allows both images to inform each other, which makes the estimation of $\*z_1$ and $\*z_2$ more robust, we show the estimated feature maps in Fig. <ref> obtained using the XCAT data. By comparing the estimated sparse feature maps $\bm{z}_{e,k}$ for $e=1,2$ and $k=1,\dots,49$ for separate reconstructions (Fig. <ref> and Fig. <ref>) and joint reconstruction (Fig. <ref> and Fig. <ref>), there are no similarities between the feature maps obtained from separate reconstructions while the feature maps obtained from joint reconstruction have similar structures.

Fig. <ref> shows the XCAT GT and the reconstruction images for both 60 keV and 120 keV energies obtained by MCAOL and the other algorithms used for comparison. The images are obtained using the parameters which corresponds to the minimum $\+{AbsBias}$ shown in Figs. <ref> and <ref>. It is worth noting that MCAOL (Fig. <ref>) manages to substantially reduce the noise as compared with CAOL (Fig. <ref>).

Fig. <ref> and Fig. <ref> show the $\+{AbsBias}$ against the $\+{STD}$ results respectively for low and high X-ray source energy. Among the methods used for comparison, TV promotes sparsity of the gradient, while JTV promotes joint sparsity of the $2$ gradients and therefore are particularly well-suited for XCAT. Despite this observation, it is possible to show that the minimum $\+{AbsBias}$ obtained by MCAOL outperforms all other algorithms, or in other words by fixing the STD, the $\+{AbsBias}$ achieved by MCAOL is always lower while it is possible to claim that by fixing the $\+{AbsBias}$, the $\+{STD}$ of MCAOL is reduced.

§.§ Results on Clinical Data

We utilized images reconstructed from data acquired on Philips IQon Spectral CT and reconstructed with a MBIR technique <cit.>.
All patients provided signed permission for the use of their clinical data for scientific purposes and anonymous publication of data.

The experiment was conducted in a similar fashion as for the XCAT simulation. We selected 22 slice pairs from a full body patient scan with $0.902$-mm pixel-width and  $512 \times 512$ image size for the training dataset corresponding to thorax. The energies used in this study are $70$ keV and $140$ keV. An additional slice pair was used to generate the projection data for reconstruction, as detailed below.

The pair of trained filters $(\*d_{1,k}, \*d_{2,k})$ obtained by both MCAOL and  CAOL unsupervised learning is shown in Fig. <ref>; we used the parameters $\gamma_1 = \gamma_2 = 10^4$ and the CAOL parameter $\alpha = 10^{-4}$.

Learned $l_{1,0}$ filters via MCAOL

($E_1$), $\gamma_1 = 10^4$     ($E_2$), $\gamma_2 = 10^4$

Learned $l_0$ filters via CAOL

($E_1$), $\alpha = 10^{-4}$     ($E_2$), $\alpha = 10^{-4}$

Learned filters $\{(\*d_{1,k}, \*d_{2,k})\}$ with $P=K=49$ using the clinical training dataset, for subfig:MCAOL:filters:XCLIN_l0_l10 MCAOL and subfig:MCAOL:filters:XCLIN_l0_l0 CAOL.

To generate the  sparse-view DECT measurements eq:mean_Poisson_DECT, we used the same geometrical and noise settings as for the XCAT simulation except that we used $451$ detector pixels and $451\times 451$ GT thorax images with attenuation coefficients $\*x^\star_1, \*x^\star_2$ at energies $140$ keV (high) and $70$ keV (low) which differs from the training examples.

In Fig. <ref> we show the estimated feature maps obtained using the clinical data. As already noted previously with the XCAT data, by comparing the estimated sparse feature maps $\bm{z}_{e,k}$ for $e=1,2$ and $k=1,\dots,49$ for separate reconstructions (Fig. <ref> and Fig. <ref>) and joint reconstruction (Fig. <ref> and Fig. <ref>), there are no similarities between the feature maps obtained from separate reconstructions while the feature maps obtained from joint reconstruction have similar structures.

$\bm{z}_{1,k}$, CAOL

$\bm{z}_{2,k}$, CAOL

$\bm{z}_{1,k}$, MCAOL

$\bm{z}_{2,k}$, MCAOL

Clinical data: estimated sparse feature maps $\bm{z}_{e,k}$ for $e=1,2$ and $k=1,\dots,49$ using CAOL (subfig:z1_sep_clinical and subfig:z2_sep_clinical) and MCAOL (subfig:z1_joint_clinical and subfig:z2_joint_clinical); color scale: red for positive values, blue for negative values.

= 8pt

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

No prior

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

Comparison of reconstructed clinical data from different reconstruction methods for sparse-view CT with top row corresponding to high energy $E_1 = 140$ keV and bottom row to low energy $E_2 = 70$ keV: sugfig:clin_gt GT image, sugfig:clin_lbfgs minimization of the NLL function without prior, sugfig:clin_caol CAOL reconstruction, sugfig:clin_mcaol MCAOL joint reconstruction, sugfig:clin_tv separate reconstruction using TV prior and sugfig:clin_jtv joint reconstruction using JTV prior.

$70$ keV.

grid = major,
legend columns=2,
legend cell align=left,
legend entries=MCAOL,CAOL, TV, JTV,
legend style=legend pos=north west
[color=blue,mark=*] table[x=xMCAOL, y=yMCAOL] Fig_9a.txt;
[color=black,mark=square*] table[x=xCAOL, y=yCAOL] Fig_9a.txt;
[color=red,mark=diamond*] table[x=xTV, y=yTV] Fig_9a.txt;
[color=green,mark=pentagon*] table[x=xJTV, y=yJTV] Fig_9a.txt;

$140$ keV.

grid = major,
legend columns=2,
legend cell align=left,
legend entries=MCAOL, CAOL, TV, JTV,
[color=blue,mark=*] table[x=xMCAOL, y=yMCAOL] Fig_9b.txt;
[color=black,mark=square*] table[x=xCAOL, y=yCAOL] Fig_9b.txt;
[color=red,mark=diamond*] table[x=xTV, y=yTV] Fig_9b.txt;
[color=green,mark=pentagon*] table[x=xJTV, y=yJTV] Fig_9b.txt;

Plot of the mean absolute bias ($\+{AbsBias}$) versus the standard deviation ($\+{STD}$) for the  clinical data at subfig:Bias_vs_std_clin_e1 low X-ray source energy ($70$ keV) and subfig:Bias_vs_std_clin_e2 high X-ray source energy ($140$ keV).

In Fig. <ref> the GT image and the reconstruction images for both energies and the different methods are shown; it is worth noting that the MCAOL reconstruction (Fig. <ref>) is less noisy than the CAOL reconstruction (Fig. <ref>).

Figs. <ref> and <ref> report the $\+{AbsBias}$ versus the $\+{STD}$ plots and we obtain a similar behavior compared to the XCAT simulations; although the relative distance of the $\+{AbsBias}$ among the simulated algorithms is reduced, MCAOL is still outperforming all other methods constantly either by fixing $\+{AbsBias}$ or $\+{STD}$, while the performance of CAOL is improving and it is close to the JTV solution accuracy.

§.§ Results for Low-Dose DECT

We conducted DECT reconstruction on the same set of data as in Section <ref> but with a different CT acquisition setup; we substantially decreased the initial photon counts to $\bar{S}_e = 10^3$ (reduction of 2 orders of magnitude compared to the previous experiments) and we doubled the number of views to 120. By approximating the total delivered X-ray dose as the product of the photons intensity times the number of views, it turns out that this scenario is considerably more challenging in terms of ill-posed problem with a total dose reduction of $50$ times.

= 8pt

No prior

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

[spy using outlines=rectangle,yellow,magnification=1.25,size=10mm,connect spies]
on (0.6,-0.60) in node [left] at (-0.3,-1);

Comparison of reconstructed clinical data from different reconstruction methods for low-dose CT with top row corresponding to high energy $E_1 = 140$ keV and bottom row to low energy $E_2 = 70$ keV: sugfig:clin_lbfgs_ld minimization of the NLL cost function without prior, sugfig:clin_caol_ld CAOL reconstruction, sugfig:clin_mcaol_ld MCAOL joint reconstruction, sugfig:clin_tv_ld separate reconstruction using TV prior, sugfig:clin_jtv_ld JTV prior and sugfig:clin_caol_pwls_ld CAOL-PWLS reconstruction.

$70$ keV.

grid = major,
legend columns=2,
legend cell align=left,
legend entries=MCAOL, CAOL, TV, JTV, CAOL-PWLS,
legend style=at=(0.8,1.05),anchor=north
[color=blue,mark=*] table[x=xMCAOL, y=yMCAOL] Fig_11a.txt;
[color=black,mark=square*] table[x=xCAOL, y=yCAOL] Fig_11a.txt;
[color=red,mark=triangle*] table[x=xTV, y=yTV] Fig_11a.txt;
[color=green,mark=square*] table[x=xJTV, y=yJTV] Fig_11a.txt;
[color=magenta,mark=pentagon*] table[x=xCAOLPWLS, y=yCAOLPWLS] Fig_11a.txt;

$140$ keV.

grid = major,
legend columns=2,
legend cell align=left,
legend entries=MCAOL, CAOL, TV, JTV, CAOL-PWLS,
legend style=at=(0.8,1.05),anchor=north
[color=blue,mark=*] table[x=xMCAOL, y=yMCAOL] Fig_11b.txt;
[color=black,mark=square*] table[x=xCAOL, y=yCAOL] Fig_11b.txt;
[color=red,mark=triangle*] table[x=xTV, y=yTV] Fig_11b.txt;
[color=green,mark=square*] table[x=xJTV, y=yJTV] Fig_11b.txt;
[color=magenta,mark=pentagon*] table[x=xCAOLPWLS, y=yCAOLPWLS] Fig_11b.txt;

Plot of the mean absolute bias ($\+{AbsBias}$) versus the standard deviation ($\+{STD}$) for the low-dose ($I_0 = 10^3$) reconstruction with clinical data at subfig:Bias_vs_std_clin_ld_e1 low X-ray source energy ($70$ keV) and subfig:Bias_vs_std_clin_ld_e2 high X-ray source energy ($140$ keV).

We used this simulation to prove that MCAOL returns a more accurate solution compared to other priors. Furthermore, we prove that despite the higher computational complexity to minimize the exact Poisson NLL in eq:NLL_MCAOLPoisson compared to solving the problem with a weighted least-squares approximated NLL, i.e., PWLS data-fit cost function, MCAOL achieves substantial improved bias accuracy compared to the PWLS solution.
To perform these experiments, we used the same optimal learned convolutional filters as obtained by the MCAOL training procedure detailed in Section <ref> and the GT images in Fig. <ref>(a).

Fig. <ref> show the reconstruction images for both energies and different methods; MCAOL (Fig. <ref>) accurately reconstruct the image features compared to all other methods and it is confirmed that the PWLS model performs poorly (Fig. <ref>).

Figs. <ref> and <ref> show either that MCAOL is consistently outperforming the other methods in terms of accuracy and variance and that the Poisson NLL formulation leads to a noticeable improvement compared to the PWLS formulation as it is indicated by comparing CAOL and CAOL-PWLS.

§ DISCUSSION

The bias-variance trade-off analysis of the estimation results over the regularization parameters confirms that MCAOL allows to achieve the minimum absolute bias compared to CAOL and other MBIR state-of-the-art methods and also reduce standard deviation.
Furthermore, MCAOL has the benefit of requiring less memory respect to DL methods because of the convolutional structure of the trained filters.

The MCAOL algorithm can be potentially extended to multi-energies $e=1,\ldots, E$ since the $l_{1,0}$ semi-norm in Eq. (<ref>) can be defined for a set of $E$-vectorised feature maps $\*z_1, \ldots, \*z_E$. Each $\*z_e, e=1, \ldots,E$ is a column vector of dimension $J\times 1$. Then the joint $l_{1,0}$ semi-norm is defined as

\begin{equation}
\| (\*z_1, \ldots, \*z_E) \|_{1,0} = \sum_{j=1}^J\*1_{]0, +\infty[}\left(|z_{1,j}| + \ldots + |z_{E,j}|\right)
\end{equation}

While the statistical noise tends to be higher in the multi-energy case, on each sub-band the contribute of the noise is reduced since the noise is split on more energy bands. Therefore, we believe that evaluating the joint norm, i.e., non-zeros elements in the feature vectors in overlapping positions for all energies at the same time, will reduce the degradation due to the increased overall noise.

The MCAOL framework allows to utilize any mixed norms for the jointly sparse regularization and other norms, such as the $l_{2,1}$-norm which as proposed by degraux2017online for convolutional synthesis operator learning, may also be considered.

In our experiments we have considered the product between the X-ray source intensity and the number of projection angles as an empirical measure for the total transmitted X-ray dose. While this metric gives a good approximation of the dose, we consider the analysis of the standardized measure of radiation dose, i.e., CT dose index (CTDI), as well as the absorbed dose as a follow-up study.

We account the open problems of how to optimally select both the regularization norm and regularization parameter according to the dataset for future algorithm development.

Although this work focuses on the multi-channel imaging reconstruction problem, we believe that our proposed method can be utilized in conjunction to DECT to task-oriented material decomposition problems. In particular, while an approach would be to design a material decomposition module in the image space which takes as input the MCAOL reconstructed images, a more compelling strategy would be designing a direct approach from sinograms to material images through MCAOL.

Furthermore, MCAOL method can be exploited for other multimodal imaging application such as PET/CT and PET/MRI. In the multimodal case, given the different intensity range on each channel, a further analysis on how to choose the NLL weights $\gamma_1\neq\gamma_2$ in (<ref>) should be conducted to properly balancing the information coming from the different modalities.

Finally, from a learning point of view, MCAOL training can be seen as a multi-channel single layer unsupervised convolutional autoencoder <cit.> which paves the way to extend this approach to deeper autoencoder architectures to capture more complex features such as textures.
The analysis and comparison of the proposed MCAOL approach with other supervised deep learning approaches is planned as a follow-up study. It is important to stress that MCAOL inherits a precise mathematical derivation and therefore it should not be susceptible of instabilities in the reconstruction which have been proven to occur with deep learning methods <cit.>.

We consider these problems as future development of the proposed algorithm.

§ CONCLUSION

In this work, we have extended the convolutional analysis operator framework to multi-channel imaging and we have applied and extensively analyzed the proposed method to the DECT application. The presented results show that by using the information coming from both energies and allowing the channels to “talk to each other” a more accurate solution of the reconstruction problem can be achieved together with a reduction of the noise in the estimate. The coupling between energies is encapsulated by using an $l_{1,0}$ sparse mixed norm in the MCAOL optimization problems both for training and reconstruction. We obtain consistently better performances across different DECT acquisition scenarios from sparse-views to low-dose photon counts.

This work was funded by the French National Research Agency (ANR) under grant .

§ REFERENCES

[Aharon et al.]Aharon, Elad   Bruckstein2006aharon2006k
Aharon, M., Elad, M.  Bruckstein, A.  2006.
K-SVD: An algorithm for designing overcomplete dictionaries for
sparse representation, IEEE Transactions on signal processing 54(11): 4311–4322.

[Antun et al.]Antun, Renna, Poon, Adcock   Hansen2020antun2020instabilities
Antun, V., Renna, F., Poon, C., Adcock, B.  Hansen, A. C.
On instabilities of deep learning in image reconstruction and the
potential costs of AI, Proceedings of the National Academy of
Sciences 117(48): 30088–30095.

[Block et al.]Block, Uecker   Frahm2007block2007undersampled
Block, K. T., Uecker, M.  Frahm, J.  2007.
Undersampled radial mri with multiple coils. iterative image
reconstruction using a total variation constraint, Magnetic Resonance in
Medicine: An Official Journal of the International Society for Magnetic
Resonance in Medicine 57(6): 1086–1098.

[Bousse et al.]Bousse, Courdurier, Émond, Thielemans, Hutton,
Irarrazaval  Visvikis2020bousse2020hypo
Bousse, A., Courdurier, M., Émond, E. C., Thielemans, K., Hutton, B. F.,
Irarrazaval, P.  Visvikis, D.  2020.
PET reconstruction with non-negativity constraints in projection
space: Optimization through hypo-convergence, IEEE Transactions on
Medical Imaging 39(1): 75–86.

Chun, I. Y.  2019.
Convolt: CONVolutional Operator Learning Toolbox.

Chun  Fessler2017Chun2017
Chun, I. Y.  Fessler, J. A.  2017.
Convolutional dictionary learning: Acceleration and convergence, IEEE Transactions on Image Processing 27(4): 1697–1712.

Chun  Fessler2019Chun2019
Chun, I. Y.  Fessler, J. A.  2019.
Convolutional analysis operator learning: Acceleration and
convergence, IEEE Transactions on Image Processing 29: 2108–2122.

[Cueva et al.]Cueva, Meaney, Siltanen   Ehrhardt2021cueva2021synergistic
Cueva, E., Meaney, A., Siltanen, S.  Ehrhardt, M. J.
Synergistic multi-spectral ct reconstruction with directional total
variation, arXiv preprint arXiv:2101.01834 .

De Pierro1995DePierro1995
De Pierro, A. R.  1995.
A modified expectation maximization algorithm for penalized
likelihood estimation in emission tomography, IEEE Transactions on
Medical Imaging 14(1): 132–137.

[Degraux et al.]Degraux, Kamilov, Boufounos   Liu2017degraux2017online
Degraux, K., Kamilov, U. S., Boufounos, P. T.  Liu, D.
Online convolutional dictionary learning for multimodal imaging, 2017 IEEE International Conference on Image Processing (ICIP),
pp. 1617–1621.

[Dong et al.]Dong, Niu  Zhu2014dong2014combined
Dong, X., Niu, T.  Zhu, L.  2014.
Combined iterative reconstruction and image-domain decomposition for
dual energy CT using total-variation regularization, Medical physics
41(5): 051909.

[Ehrhardt et al.]Ehrhardt, Thielemans, Pizarro, Atkinson,
Ourselin, Hutton  Arridge2014ehrhardt2014
Ehrhardt, M. J., Thielemans, K., Pizarro, L., Atkinson, D., Ourselin, S.,
Hutton, B. F.  Arridge, S. R.  2014.
Joint reconstruction of pet-mri by exploiting structural similarity,
Inverse Problems 31(1): 015001.

Elbakri  Fessler2002Elbakri2002
Elbakri, I. A.  Fessler, J. A.  2002.
Statistical image reconstruction for polyenergetic X-ray computed
tomography, IEEE Transactions on Medical Imaging 21(2): 89–99.

[Fang et al.]Fang, Wu, Kim, Kalra, Singh, Li   Li2021fang2021
Fang, W., Wu, D., Kim, K., Kalra, M. K., Singh, R., Li, L.  Li, Q.
Iterative material decomposition for spectral ct using
self-supervised noise2noise prior, Physics in Medicine & Biology .

Forghani   Mukherji2018forghani2018characterization
Forghani, R.  Mukherji, S. K.  2018.
Advanced dual-energy CT applications for the evaluation of the soft
tissues of the neck, Clinical Radiology 73(1): 70–80.

Garcia-Cardona   Wohlberg2018garcia2018convolutional_p
Garcia-Cardona, C.  Wohlberg, B.  2018.
Convolutional dictionary learning for multi-channel signals, 2018 52nd Asilomar Conference on Signals, Systems, and Computers, IEEE,
pp. 335–342.

[Kazantsev et al.]Kazantsev, Jørgensen, Andersen, Lionheart,
Lee  Withers2018kazantsev2018joint
Kazantsev, D., Jørgensen, J. S., Andersen, M. S., Lionheart, W. R. B., Lee,
P. D.  Withers, P. J.  2018.
Joint image reconstruction method with correlative multi-channel
prior for x-ray spectral computed tomography, Inverse Problems 34(6): 064001.

[Kim et al.]Kim, Ye, Worstell, Ouyang, Rakvongthai, El Fakhri
Kim, K., Ye, J. C., Worstell, W., Ouyang, J., Rakvongthai, Y., El Fakhri, G.
Li, Q.  2014.
Sparse-view spectral CT reconstruction using spectral patch-based
low-rank penalty, IEEE Transactions on Medical Imaging 34(3): 748–760.

[Kim et al.]Kim, Ye, Worstell, Ouyang, Rakvongthai, Fakhri
Kim, K., Ye, J. C., Worstell, W., Ouyang, J., Rakvongthai, Y., Fakhri, G. E.
Li., Q.  2015.
Sparse-view spectral CT reconstruction using spectral patch-based
low-rank penalty, IEEE Transactions on Medical Imaging 34(3): 748–760.

[McCollough et al.]McCollough, Leng, Yu   Fletcher2015mccollough2015dual
McCollough, C. H., Leng, S., Yu, L.  Fletcher, J. G.
Dual-and multi-energy CT: principles, technical approaches, and
clinical applications, Radiology 276(3): 637–653.

Nocedal  Wright2006NumOpti2006
Nocedal, J.  Wright, S. J.  2006.
Numerical Optimization, 2nd edn, Springer Science & Business

Philips IQon Elite Spectral CT product
Philips IQon Elite Spectral CT product specifications  2018.

Ravishankar  Bresler2012ravishankar2012
Ravishankar, S.  Bresler, Y.  2012.
Learning sparsifying transforms, IEEE Transactions on Signal
Processing 61(5): 1072–1086.

[Ravishankar et al.]Ravishankar, Nadakuditi   Fessler2017ravishankar2017efficient
Ravishankar, S., Nadakuditi, R. R.  Fessler, J. A.
Efficient sum of outer products dictionary learning (SOUP-DIL)
and its application to inverse problems, IEEE Transactions on
Computational Imaging 3(4): 694–709.

[Ravishankar et al.]Ravishankar, Ye   Fessler2019ravishankar2019
Ravishankar, S., Ye, J. C.  Fessler, J. A.  2019.
Image reconstruction: From sparsity to data-adaptive methods and
machine learning, Proceedings of the IEEE 108(1): 86–109.

[Sajja et al.]Sajja, Lee, Eriksson, Nordström, Sahgal,
Hashemi, Mainprize  Ruschin2020sajja2020technical
Sajja, S., Lee, Y., Eriksson, M., Nordström, H., Sahgal, A., Hashemi, M.,
Mainprize, J. G.  Ruschin, M.  2020.
Technical principles of dual-energy cone beam computed tomography and
clinical applications for radiation therapy, Advances in radiation
oncology 5(1): 1–16.

[Segars et al.]Segars, Mahesh, Beck, Frey   Tsui2008segars2008xcat
Segars, W. P., Mahesh, M., Beck, T. J., Frey, E. C.  Tsui, B. M. W.
Realistic CT simulation using the 4D XCAT phantom, Medical
Physics 35(8): 3800–3808.

[Sidky et al.]Sidky, Kao   Pan2006Sidky2006Accurate
Sidky, E. Y., Kao, C.  Pan, X.  2006.
Accurate image reconstruction from few-views and limited-angle data
in divergent-beam CT, J. X-ray Sci. Technol 14(2): 119–139.

Sidky  Pan2008sidky2008image
Sidky, E. Y.  Pan, X.  2008.
Image reconstruction in circular cone-beam computed tomography by
constrained, total-variation minimization, Physics in Medicine &
Biology 53(17): 4777.

[Song et al.]Song, Weizman, Mota, Eldar   Rodrigues2019song2019coupled
Song, P., Weizman, L., Mota, J. F., Eldar, Y. C.  Rodrigues, M. R.
Coupled dictionary learning for multi-contrast mri reconstruction,
IEEE Transactions on Medical Imaging 39(3): 621–633.

[van Elmpt et al.]van Elmpt, Landry, Das   Verhaegen2016van2016dual
van Elmpt, W., Landry, G., Das, M.  Verhaegen, F.  2016.
Dual energy ct in radiotherapy: current applications and future
outlook, Radiotherapy and Oncology 119(1): 137–144.

[Whiting et al.]Whiting, Massoumzadeh, Earl, O'Sullivan, Snyder
Whiting, B., Massoumzadeh, P., Earl, O., O'Sullivan, J., Snyder, D.
Williamson, J.  2006.
Properties of preprocessed sinogram data in X-ray CT, Medical physics 33(3): 3290–3303.

[Wu et al.]Wu, Zhang, Wang, Liu, Chen   Yu.2018wu2018dictionary
Wu, W., Zhang, Y., Wang, Q., Liu, F., Chen, P.  Yu., H.
Low-dose spectral CT reconstruction using image gradient
$\ell$0–norm and tensor dictionary, Applied Mathematical Modelling
63: 538–557.

[Xu et al.]Xu, Lu, Xu  Jia2011xu2011image
Xu, L., Lu, C., Xu, Y.  Jia, J.  2011.
Image smoothing via $l_0$ gradient minimization, 2011 SIGGRAPH
Asia Conference, pp. 1–12.

[Xu et al.]Xu, Yu, Mou, Zhang, Hsieh   Wang2012xu2012
Xu, Q., Yu, H., Mou, X., Zhang, L., Hsieh, J.  Wang, G.
Low-dose X-ray CT reconstruction via dictionary learning, IEEE Transactions on Medical Imaging 31(9): 1682–1697.

[Yang et al.]Yang, Cong  Wang2017yang2017prism
Yang, Q., Cong, W.  Wang, G.  2017.
Superiorization-based multi-energy CT image reconstruction, Inverse Problems 33(4).

[Yao et al.]Yao, Li  Chen2019yao2019
Yao, Y., Li, L.  Chen, Z.  2019.
Dynamic-dual-energy spectral ct for improving multi-material
decomposition in image-domain, Physics in Medicine & Biology 64(13): 135006.

[Yu et al.]Yu, Leng  McCollough2012yu2012mono
Yu, L., Leng, S.  McCollough, C. H.  2012.
Dual-energy CT–based monochromatic imaging, American Journal
of Roentgenology 199(5): S9–S15.

[Zhang et al.]Zhang, Liang, Wang, Cai, Wang, Tang, Zheng, Li, Yan
Zhang, W., Liang, N., Wang, Z., Cai, A., Wang, L., Tang, C., Zheng, Z., Li, L.,
Yan, B.  Hu, G.  2020.
Multi-energy CT reconstruction using tensor nonlocal similarity and
spatial sparsity regularization, Quantitative Imaging in Medicine and
Surgery 10(10): 1940.

[Zheng et al.]Zheng, Ravishankar, Long   Fessler2018zheng2018pwls
Zheng, X., Ravishankar, S., Long, Y.  Fessler, J. A.
PWLS-ULTRA: An efficient clustering and learning-based approach for
low-dose 3D CT image reconstruction, IEEE Transactions on Medical
Imaging 37(6): 1498–1510.

[Zhu et al.]Zhu, Byrd, Lu   Nocedal1997zhu1997algorithm
Zhu, C., Byrd, R. H., Lu, P.  Nocedal, J.  1997.
Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale
bound-constrained optimization, ACM Transactions on Mathematical
Software (TOMS) 23(4): 550–560.