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

Alessandro Perelli, Suxer Alfonso Garcia, Alexandre Bousse, Jean-Pierre Tasu, Nikolaos Efthimiadis, Dimitris Visvikis

► To cite this version:

Alessandro Perelli, Suxer Alfonso Garcia, Alexandre Bousse, Jean-Pierre Tasu, Nikolaos Efthimiadis, et al.. Multi-Channel Convolutional Analysis Operator Learning for Dual-Energy CT Reconstruction. Physics in Medicine and Biology, 2022, 67 (6), pp.065001. 10.1088/1361-6560/ac4c32. hal-04994416

HAL Id: hal-04994416

https://hal.science/hal-04994416v1

Submitted on 17 Mar 2025

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>

Distributed under a Creative Commons Attribution 4.0 International License
# Multi-Channel Convolutional Analysis Operator Learning for Dual-Energy CT Reconstruction

Alessandro Perelli*¹,², Suxer Alfonso Garcia*¹,

Alexandre Bousse¹, Jean-Pierre Tasu³, Nikolaos Efthimiadis³,

Dimitris Visvikis¹

¹ LaTIM, INSERM UMR 1101, Université de Bretagne Occidentale, 29238 Brest, France

² School of Science and Engineering, University of Dundee, Scotland DD1 4HN, UK

³ Department of Radiology, University Hospital Poitiers, Poitiers, France

* These authors equally contributed to this work.

E-mail: aperelli001@dundee.ac.uk

E-mail: suxer-lazara.alfonso-garcia@inserm.fr

**Abstract.** *Objective.* Dual-energy computed tomography (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 multi-channel convolutional analysis operator learning (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 convolutional analysis operator learning (CAOL) algorithm.

*Main results.* Extensive experiments with simulated and real computed tomography (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 total-variation (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 model-based iterative reconstruction (MBIR) techniques, thus paving the way for dose reduction.

## 1. Introduction

Dual-energy computed tomography (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 (Zhang, Liang, Wang, Cai, Wang, Tang,
Zheng, Li, Yan & Hu 2020). Since its invention, the technique has progressively been used for material decomposition and energy-selective imaging (Dong, Niu & Zhu 2014). Applications of DECT include differentiation and quantification of materials, tissue characterization, virtual monoenergetic imaging, automated bone removal in computed tomography (CT) angiography, cardiovascular imaging, multiple contrast agent imaging and mapping of effective atomic number, and radiotherapy (McCollough, Leng, Yu & Fletcher 2015, Forghani & Mukherji 2018, Yu, Leng & McCollough 2012, van Elmpt, Landry, Das & Verhaegen 2016).

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 (Sajja, Lee, Eriksson, Nordström, Sahgal, Hashemi, Mainprize & Ruschin 2020).

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 signal-to-noise ratio (SNR) (Zhang et al. 2020).

In the literature, most of the development on low-dose CT reconstruction has focused on single image. Among the main techniques, model-based iterative reconstruction (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 (Elbakri & Fessler 2002). Typically, such cost functions consist of a data-fidelity term, e.g., least squares or negative log-likelihood (NLL), capturing the imaging forward model and the measurement/noise statistical model and a regularizer term promoting smoothness, low-rank or sparsity (Kim, Ye, Worstell, Ouyang, Rakvongthai, El Fakhri & Li 2014). Total-variation (TV) (Sidky, Kao & Pan 2006, Sidky & Pan 2008) 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 (Block, Uecker & Frahm 2007). Recent developments and interests within spectral CT refer to self-supervised methods using only noisy data (Fang, Wu, Kim, Kalra, Singh, Li & Li 2021) and dynamic DECT, which refers to randomly changing the energy threshold of the detectors to obtain the spectral information (Yao, Li & Chen 2019).

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 dictionary learning (DL) (Xu, Yu, Mou, Zhang,
Hsieh & Wang 2012) 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 (Aharon, Elad & Bruckstein 2006). DL-based image reconstruction integrates the learned atoms with the raw scanner data within a regularized MBIR context (Ravishankar, Nadakuditi & Fessler 2017, Zheng, Ravishankar, Long & Fessler 2018). Other closely related methods include sparsifying transform learning (Ravishankar & Bresler 2012) and the connection between data-adaptive models and convolutional deep learning algorithms (Ravishankar, Ye & Fessler 2019) 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 3-dimensional (3-D) multimodal imaging. To address these problems, convolutional dictionary learning (CDL) techniques utilize shift-invariant filters, providing a convenient and memory-efficient alternative to conventional DL techniques (Chun & Fessler 2017). 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 (Chun & Fessler 2019) (sparse convolution) and is known as convolutional analysis operator learning (CAOL).

Despite the rapidly expanding research, the application of CDL to multi-channel images has received little attention (Degraux, Kamilov, Boufounos & Liu 2017, Garcia-Cardona & Wohlberg 2018). 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 joint total variation (JTV) (Ehrhardt, Thielemans, Pizarro, Atkinson, Ourselin, Hutton & Arridge 2014, Cueva, Meaney, Siltanen & Ehrhardt 2021), spectral patch-based penalty for the maximum-likelihood method (Kim, Ye, Worstell, Ouyang, Rakvongthai, Fakhri & Li. 2015), tensor-based and coupled dictionary learning (Wu, Zhang, Wang, Liu, Chen & Yu. 2018, Song, Weizman, Mota, Eldar & Rodrigues 2019), parallel level sets (Kazantsev, Jørgensen, Andersen, Lionheart, Lee & Withers 2018) and the prior rank, intensity and sparsity model (PRISM) (Yang, Cong & Wang 2017).

## 1.1. Main Contribution

We extend the CAOL approach to multichannel settings and we develop a multi-channel convolutional analysis operator learning (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 multi-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.

## 1.2. Notation and Paper Organization

The remainder of this paper is structured as follows: Section 2 briefly reviews the single-channel CAOL algorithm as initially proposed in (Chun & Fessler 2019) while Section 3 presents the proposed approach to extend CAOL to multi-channel problems. In Section 4 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 5 comprehensive results of MCAOL on a numerical extended cardiac-torso (XCAT) phantom (Segars, Mahesh, Beck, Frey & Tsui 2008) 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; $[\mathbf{a}]_j$ denotes the j-th entry of **a**; ‘(·)ᵀ’ denotes the transposition; ‘⊗’ denotes the 2-dimensional (2-D) convolution operator; given an image $\mathbf{x} \in \mathbb{R}^J$ and a filter $\mathbf{d} \in \mathbb{R}^P$, $P < J$, $\mathbf{d} \otimes \mathbf{x}$ represents the 2-D 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 $\boldsymbol{x} \in \mathbb{R}^J$ is represented by a column vector for algebraic operations and by a $\sqrt{J} \times \sqrt{J}$ matrix for 2-D convolutions. Finally, we indicate the set of positive real numbers with the notation $]0, +\infty[$.

## 2. 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

$$ \min_{\boldsymbol{x} \in \mathbb{R}^J} L(\boldsymbol{x}, \boldsymbol{y}) + \beta R(\boldsymbol{x}) \quad (1) $$

where $\boldsymbol{x} \in \mathbb{R}^J$ is the 2-D image to reconstruct, $\boldsymbol{y} \in \mathbb{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 expectation-maximization (EM) for emission tomography (ET) (De Pierro 1995) or penalized weighted least-squares (PWLS) combined with separable paraboloidal surrogate (SPS) for CT (Elbakri & Fessler 2002).

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

We consider the CAOL approach (Chun & Fessler 2019) where the image is sparsely represented with convolutional kernels (filters). In the analysis model, the image is represented with “sparsifying” filters $\boldsymbol{d}_k \in \mathbb{R}^P$ by the analysis operator $\mathcal{A}_D: \boldsymbol{x} \mapsto \{\boldsymbol{d}_k \otimes \boldsymbol{x}\}$, such that

$$ \boldsymbol{d}_k \otimes \boldsymbol{x} \approx \boldsymbol{z}_k, \quad \forall k = 1, \dots, K. \quad (2) $$

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

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

$$ \mathcal{D}^* = \arg \min_{\mathcal{D} \in \mathcal{C}} \min_{\{\boldsymbol{z}_{l,k}\}} F_a(\mathcal{D}, \{\boldsymbol{z}_{l,k}\}) \quad (3) $$

where $\mathcal{C}$ is the constraints on $\mathcal{D} = \{\boldsymbol{d}_k\}$ and with the training analysis objective function $F_a$ defined as

$$ F_a(\mathcal{D}, \{\boldsymbol{z}_{l,k}\}) = \sum_{l=1}^{L} \sum_{k=1}^{K} \frac{1}{2} \| \boldsymbol{d}_k \otimes \boldsymbol{x}_l - \boldsymbol{z}_{l,k} \|_2^2 + \alpha \| \boldsymbol{z}_{l,k} \|_0 \quad (4) $$
where $z_{l,k} \in \mathbb{R}^J$ is the feature image associated to the training image $x_l$ and the filter $d_k$, $\|\cdot\|_0$ is the sparsity-promoting $l_0$ semi-norm defined for all $z = [z_1, \dots, z_J]^\top \in \mathbb{R}^J$ as

$$ \|z\|_0 = \sum_{j=1}^{J} \mathbf{1}_{[0, +\infty[}(|z_j|)} \quad (5) $$

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

$$ C = \left\{ \{d_k\} : [d_1, \dots, d_K][d_1, \dots, d_K]^\top = \frac{1}{P} I_K \right\} \quad (6) $$

where $I_K$ is the $K \times K$ identity matrix, to promote filters diversity. The optimization problem (3) is solved by alternating between minimization in $z$ and in $D$. The joint minimization in $D$, $z$ can be achieved with a *Block Proximal Extrapolated Gradient method using a Majorizer* (BPEG-M) algorithm, which can be implemented using the CONVolutional Operator Learning Toolbox (CONVOLT) (Chun & Fessler 2019, Chun 2019). BPEG-M alternates between minimization in $z$ and $D$. The minimization in $z$ is obtained with a hard-thresholding operator $\mathcal{T}: \mathbb{R}^J \times \mathbb{R}_+^* \to \mathbb{R}^J$ defined at each row $j$ as

$$ [\mathcal{T}(\mathbf{a}, \beta)]_j = \begin{cases} a_j & \text{for } \frac{1}{2}a_j^2 \ge \beta \\ 0 & \text{otherwise} \end{cases} \quad (7) $$

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

$$ \mathcal{T}(d_k \otimes x_l, \alpha) = \arg \min_{z_{l,k}} \frac{1}{2} \|d_k \otimes x_l - z_{l,k}\|_2^2 + \alpha \|z_{l,k}\|_0. \quad (8) $$

The minimization in $D$ is achieved with a *proximal mapping with orthogonality constraints* (PMOC), which can be implemented with a singular value decomposition (Chun & Fessler 2019). Finally the regularizer $R$ in the minimization problem (1) is derived from the learned filters $D^*$ as

$$ R(\mathbf{x}) = \min_{\{z_k\}} \sum_{k=1}^{K} \frac{1}{2} \|d_k^* \otimes \mathbf{x} - z_k\|_2^2 + \alpha \|z_k\|_0. \quad (9) $$

### 3. 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 \mathbb{R}^J$ of the same “object” from two independent measurements $y_1 \in \mathbb{R}^{I_1}$ and $y_2 \in \mathbb{R}^{I_2}$ corresponding to two modalities, multi-channel MBIR can be achieved by using an iterative algorithm to solve

$$ \min_{x_1, x_2 \in \mathbb{R}^J} \rho_1 L_1(x_1, y_1) + \rho_2 L_2(x_2, y_2) + R_{\text{mc}}(x_1, x_2) \quad (10) $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Diagram of MCAOL consisting of the unsupervised filter learning phase and the model-based iterative DECT reconstruction module.</figcaption></figure>

where $L_1$ and $L_2$ are the data-fidelity terms for $x_1$ and $x_2$, $R_{mc}$ is a multi-channel regularizer and $\rho_1, \rho_2 > 0$ are weights. $R_{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 Ehrhardt et al. (2014).

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. 1 for DECT. Let $\{(x_{1,l}, x_{2,l}) \in \mathbb{R}^J \times \mathbb{R}^J : l = 1, \dots, L\}$ be a training dataset consisting of $L$ pairs of images. MCAOL learns the sparsifying filter pairs

$$ (\mathbf{d}_{1,k}, \mathbf{d}_{2,k}) \in \mathbb{R}^P \times \mathbb{R}^P : k = 1, \dots, K \qquad (11) $$

together with the extracted feature pairs

$$ (\mathbf{z}_{1,l,k}, \mathbf{z}_{2,l,k}) \in \mathbb{R}^J \times \mathbb{R}^J : k = 1, \dots, K, \quad l = 1, \dots, L. \qquad (12) $$

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

$$ (\mathbf{D}_1^*, \mathbf{D}_2^*) = \arg\min_{\mathbf{D}_1, \mathbf{D}_2 \in \mathcal{C}} \min_{\{z_{1,l,k}\}, \{z_{2,l,k}\}} F_{mc}(\mathbf{D}_1, \mathbf{D}_2, \{z_{1,l,k}\}, \{z_{2,l,k}\}) \qquad (13) $$

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

$$ F_{mc}(\mathbf{D}_1, \mathbf{D}_2, \{z_{1,l,k}\}, \{z_{2,l,k}\}) = \sum_{l=1}^{L} \sum_{k=1}^{K} \frac{\gamma_1}{2} \|\mathbf{d}_{1,k} \oplus \mathbf{x}_{1,l} - \mathbf{z}_{1,l,k}\|_2^2 \\ + \frac{\gamma_2}{2} \|\mathbf{d}_{2,k} \oplus \mathbf{x}_{2,l} - \mathbf{z}_{2,l,k}\|_2^2 + \|(z_{1,l,k}, z_{2,l,k})\|_{1,0} \qquad (14) $$

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

$$ \|(z_1, z_2)\|_{1,0} = \sum_{j=1}^{J} \mathbf{1}_{[0, +\infty[} (|z_{1,j}| + |z_{2,j}|) \qquad (15) $$

$\| \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 $\mathbf{D}_1, \mathbf{D}_2$.
**Algorithm 1:** MCAOL Training Algorithm

**Input:** DE Training Dataset $\mathbf{x}_{e,l}$, $l = 1, \dots, L$, $e = 1, 2$, joint sparsity weights $\gamma = (\gamma_1, \gamma_2)$, #outer iterations $N_{\text{outer}}$

**Output:** Learned filters $(\mathbf{D}_1^*, \mathbf{D}_2^*)$ $(\mathbf{D}_1^0, \mathbf{D}_2^0) \leftarrow$ Normalized random initialization;

for $t = 0, \dots, N_{\text{outer}} - 1$ do
  Update sparse codes (in parallel) ;
  for $k,l=1,1,\dots,K,L$ do
    $(\mathbf{z}_{1,l,k}^{t+1}, \mathbf{z}_{2,l,k}^{t+1}) \leftarrow \mathcal{T}_{\text{mc}}(\mathbf{d}_{1,k}^{t+1} \circledast \mathbf{x}_{1,l}, \mathbf{d}_{2,k}^{t+1} \circledast \mathbf{x}_{2,l}, \boldsymbol{\gamma})$ ;
  end
  Update Filters ;
  $\mathbf{D}_1^{t+1} \leftarrow$ PMOC($\mathbf{x}_{1,l}, \mathbf{z}_{1,k}^{t+1}$) ;
  $\mathbf{D}_2^{t+1} \leftarrow$ PMOC($\mathbf{x}_{2,l}, \mathbf{z}_{2,k}^{t+1}$) ;

end
$\mathbf{D}_1^* \leftarrow \mathbf{D}_1^{N_{\text{outer}}}$ ;
$\mathbf{D}_2^* \leftarrow \mathbf{D}_2^{N_{\text{outer}}}$ ;

To solve (13) we utilize an iterative alternating minimization scheme where the update of the filters $(\mathbf{D}_1, \mathbf{D}_2)$ is solving a block non-convex problem using the PMOC algorithm (Chun & Fessler 2019, Chun 2019) while for the update of the sparse codes $(\mathbf{z}_{1,l,k}, \mathbf{z}_{2,l,k})$ we deploy a multi-channel hard-thresholding operator $\mathcal{T}_{\text{mc}}: \mathbb{R}^J \times \mathbb{R}^J \times (\mathbb{R}_+^*)^2 \to \mathbb{R}^J \times \mathbb{R}^J$ defined at each row $j$ as

$$[\mathcal{T}_{\text{mc}}(\mathbf{a}_1, \mathbf{a}_2, \boldsymbol{\gamma})]_j = \begin{cases} (a_{1,j}, a_{2,j}) & \text{for } \frac{1}{2}\gamma_1 a_{1,j}^2 + \frac{1}{2}\gamma_2 a_{2,j}^2 \ge 1 \\ (0, 0) & \text{otherwise} \end{cases} \quad (16)$$

for all $\mathbf{a}_1 = [a_{1,1}, \dots, a_{1,J}]^\top \in \mathbb{R}^J$, $\mathbf{a}_2 = [a_{2,1}, \dots, a_{2,J}]^\top \in \mathbb{R}^J$ and for all $\boldsymbol{\gamma} = (\gamma_1, \gamma_2) \in (\mathbb{R}_+^*)^2$, which provides a global minimizer for $(\mathbf{z}_1, \mathbf{z}_2) \mapsto \frac{\gamma_1}{2} \|\mathbf{a}_1 - \mathbf{z}_1\|_2^2 + \frac{\gamma_2}{2} \|\mathbf{a}_2 - \mathbf{z}_2\|_2^2 + \|\mathbf{z}_1, \mathbf{z}_2\|_{1,0}$ (Xu, Lu, Xu & Jia 2011, Section 3), in such a way that

$$\mathcal{T}_{\text{mc}}(\mathbf{d}_{1,k} \circledast \mathbf{x}_{1,l}, \mathbf{d}_{2,k} \circledast \mathbf{x}_{2,l}, \boldsymbol{\gamma}) = \arg \min_{\mathbf{z}_{1,l,k}, \mathbf{z}_{2,l,k}} \left\{ \frac{\gamma_1}{2} \|\mathbf{d}_{1,k} \circledast \mathbf{x}_{1,l} - \mathbf{z}_{1,l,k}\|_2^2 \right. \\ \left. + \frac{\gamma_2}{2} \|\mathbf{d}_{2,k} \circledast \mathbf{x}_{2,l} - \mathbf{z}_{2,l,k}\|_2^2 + \|(\mathbf{z}_{1,l,k}, \mathbf{z}_{2,l,k})\|_{1,0} \right\}. \quad (17)$$

Finally the regularizer $R_{\text{mc}}$ in the minimization problem (10) is derived from the learned filters $(\mathbf{D}_1^*, \mathbf{D}_2^*)$ as

$$R_{\text{mc}}(\mathbf{x}_1, \mathbf{x}_2) = \min_{\substack{\mathbf{z}_{1,k} \\ \mathbf{z}_{2,k}}} \sum_{k=1}^{K} \frac{\gamma_1}{2} \|\mathbf{d}_{1,k}^* \circledast \mathbf{x}_1 - \mathbf{z}_{1,k}\|_2^2 \\ + \frac{\gamma_2}{2} \|\mathbf{d}_{2,k}^* \circledast \mathbf{x}_2 - \mathbf{z}_{2,k}\|_2^2 + \|(\mathbf{z}_{1,k}, \mathbf{z}_{2,k})\|_{1,0} \quad (18)$$

The pseudo-code for MCAOL training procedure is summarized in Algorithm 1.
# 4. Dual-Energy CT Reconstruction with Multi-Channel CAOL

## 4.1. 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 2-D slice-by-slice imaging systems.

For image reconstruction we assume that the continuous attenuation image $\mu_e(\mathbf{r})$ which denotes the linear attenuation coefficient at position $\mathbf{r} \in \mathbb{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,

$$ \mu_e(\mathbf{r}) = \sum_{j=1}^{J} x_{e,j} b_j(\mathbf{r}), \qquad (19) $$

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

$$ \int_{\mathbb{R}} \mu_e(\boldsymbol{\nu}_i(l)) \, dl = \sum_{j=1}^{J} x_{e,j} \int_{\mathbb{R}} b_j(\boldsymbol{\nu}_i(l)) \, dl = \sum_{j=1}^{J} a_{i,j} x_{e,j} \qquad (20) $$

where $\boldsymbol{\nu}_i(l) = \mathbf{s}_i + l\boldsymbol{\epsilon}_i \in \mathbb{R}^2$ is a parametrisation of the $i$-th ray emitted from the source $\mathbf{s}_i$ with direction $\boldsymbol{\epsilon}_i$, $a_{i,j} \triangleq \int_{\mathbb{R}} b_j(\boldsymbol{\nu}_i(l)) \, dl$ is the contribution of the $j$-th pixel to the $i$-th ray. The system matrix $\mathbf{A}$ is constructed as an under-determined matrix of dimensions $I \times J$ where $I = N_d \times N_\theta$ with $N_d$ and $N_\theta$ being respectively the number of detectors and $N_\theta$ and the number of angles (projections), and is defined as $[\mathbf{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 $\mathbf{x}_e$ as

$$ \bar{y}_{i,e}(\mathbf{x}_e) = S_e e^{-[\mathbf{A}\mathbf{x}_e]_i} + \eta_{e,i} \qquad (21) $$

where $\mathbf{x}_e = [x_{e,1}, \dots, x_{e,J}]^\top \in \mathbb{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 \mathbb{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.,

$$ y_{i,e} \sim \text{Poisson}(\bar{y}_{i,e}(\mathbf{x}_e)) \qquad (22) $$

and the measurements at each energy bin $e = 1, 2$ are stored in a vector $\mathbf{y}_e = [y_{e,1}, \dots, y_{e,N_d \cdot N_\theta}]^\top$.

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 (Whiting, Massoumzadeh, Earl, O'Sullivan, Snyder & Williamson 2006), 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.
## 4.2. 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 $(\mathbf{y}_1, \mathbf{y}_2)$ by solving the minimization problem (10), with positivity constraints on $(\mathbf{x}_1, \mathbf{x}_2)$, using the Poisson NLL functions $L_1$ and $L_2$ defined as

$$ -L_e(\mathbf{x}_e, \mathbf{y}_e) = \sum_{i=1}^{I} y_{e,i} \log \bar{y}_{i,e}(\mathbf{x}_e) - \bar{y}_{i,e}(\mathbf{x}_e), \quad e = 1, 2 \quad (23) $$

and the trained regularizer $R_{\text{mc}}$ derived from the learned filters $(\mathbf{D}_1^*, \mathbf{D}_2^*)$ as in (18). Therefore, substituting (23) and (18) into the minimization (10), we obtain the following explicit expression for the MCAOL DECT reconstruction problem:

$$ \begin{aligned} (\mathbf{x}_1^*, \mathbf{x}_2^*) = & \arg \min_{\mathbf{x}_1, \mathbf{x}_2 \ge 0} \sum_{e=1}^{2} \rho_e \underbrace{\sum_{i=1}^{I} y_{e,i} \log \bar{y}_{i,e}(\mathbf{x}_e) - \bar{y}_{i,e}(\mathbf{x}_e)}_{L_e(\mathbf{x}_e, \mathbf{y}_e)} \\ & + \min_{\substack{\{\mathbf{z}_{1,k}\} \\ \{\mathbf{z}_{2,k}\}}} \sum_{k=1}^{K} \left\{ \sum_{e=1}^{2} \frac{\gamma_e}{2} \| \mathbf{d}_{e,k}^* \otimes \mathbf{x}_e - \mathbf{z}_{e,k} \|_2^2 \right\} + \|(\mathbf{z}_{1,k}, \mathbf{z}_{2,k})\|_{1,0} \end{aligned} \quad (24) $$

We solve the minimization problem (24) by the alternating estimation of the sparse feature images and the linear attenuation images $\{\mathbf{x}_e : e = 1, 2\}$. Given the current estimates of the sparse coefficients $\{\mathbf{z}_k^t : k = 1, \dots, K\}$, the image update $\mathbf{x}_e^t$ at iteration $t$ is obtained through the following minimization problem

$$ \mathbf{x}_e^t = \arg \min_{\mathbf{x}_e \in (\mathbb{R}^+)^J} \Phi_e^t(\mathbf{x}_e) \quad (25) $$

$$ \text{with } \Phi_e(\mathbf{x}_e) = \rho_e L_e(\mathbf{x}_e, \mathbf{y}_e) + \frac{\gamma_e}{2} \sum_{k=1}^{K} \| \mathbf{d}_{e,k}^* \otimes \mathbf{x}_e - \mathbf{z}_{e,k}^t \|_2^2. $$

In this work, we utilized a limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm (Nocedal & Wright 2006, Chapter 7) to solve (25). The L-BFGS iterative solver estimates $\mathbf{x}_e^t$ starting the previous iterate $\mathbf{x}_e^{t-1}$. We define the first estimate as $\mathbf{x}_e^{t,(0)} = \mathbf{x}_e^t$. At inner iteration $q$, given a current estimate $\mathbf{x}_e^{t,(q)}$, the new estimate $\mathbf{x}_e^{t,(q+1)}$ is obtained as

$$ \begin{aligned} \mathbf{x}_e^{t,(q+1)} &= \mathbf{x}_e^{t,(q)} - s^* \mathbf{B}^{t,(q)} \nabla \Phi_e^t(\mathbf{x}_e^{t,(q)}) \\ \text{with } s^* &= \arg \max_{s \in [0,1]} \chi(s) \end{aligned} \quad (26) $$

$$ \text{and } \chi(s) = \Phi_e^t(\mathbf{x}_e^{t,(q)} - s \mathbf{B}^{t,(q)} \nabla \Phi_e^t(\mathbf{x}_e^{t,(q)})) $$

where $\mathbf{B}^{t,(q)}$ is an approximate inverse Hessian of $\Phi_e^t$ evaluated at $\mathbf{x}_e^{t,(q)}$. The matrix/vector product $\mathbf{B}^{t,(q)} \nabla \Phi_e^t(\mathbf{x}_e^{t,(q)})$ in (26) is directly computed (without storing $\mathbf{B}^{t,(q)}$) from the $m$ previous iterates $\mathbf{x}_e^{t,(q-p)}$, $p = 0, \dots, m-1$. An approximate solution of the line-search
<caption>**Algorithm 2:** MCAOL Reconstruction Algorithm</caption>

**Input:** Initial images ($\mathbf{x}_1^0, \mathbf{x}_2^0$), DECT learned filters $\mathbf{D}^* = (\mathbf{D}_1^*, \mathbf{D}_2^*)$, joint sparsity weight $\gamma = (\gamma_1, \gamma_2)$, penalty weights $\rho = (\rho_1, \rho_2)$, DE sinogram $\mathbf{y} = (\mathbf{y}_1, \mathbf{y}_2)$, system matrix $\mathbf{A}$, intensities $(S_1, S_2)$, #outer iterations $N_{\text{outer}}$.

**Output:** Reconstructed images ($\mathbf{x}_1^*, \mathbf{x}_2^*$)

**for** $t = 0, \dots, N_{\text{outer}} - 1$ **do**

| **Update sparse codes (in parallel);**

| **for** $k, \dots, K$ **do**

| | $(\mathbf{z}_{1,k}^{t+1}, \mathbf{z}_{2,k}^{t+1}) \leftarrow \mathcal{T}_{\text{mc}}(\mathbf{d}_{1,k}^* \oplus \mathbf{x}_1^t, \mathbf{d}_{2,k}^* \oplus \mathbf{x}_2^t, \gamma)$;

| **end**

| **Update linear attenuation images;**

| | $\mathbf{x}_1^{t+1} \leftarrow \text{L-BFGS}(\Phi_1^t, \text{init} = \mathbf{x}_1^t | \mathbf{y}_1, \mathbf{z}_1^{t+1}, \mathbf{D}_1^*, \mathbf{A}, S_1, \rho_1, \gamma_1)$;

| | $\mathbf{x}_2^{t+1} \leftarrow \text{L-BFGS}(\Phi_2^t, \text{init} = \mathbf{x}_2^t | \mathbf{y}_2, \mathbf{z}_2^{t+1}, \mathbf{D}_1^*, \mathbf{A}, S_2, \rho_2, \gamma_2)$;

| **end**

| $\mathbf{x}_1^* \leftarrow \mathbf{x}_1^{N_{\text{outer}}}$;

| $\mathbf{x}_2^* \leftarrow \mathbf{x}_2^{N_{\text{outer}}}$;

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 Bousse, Courdurier, Émond, Thielemans, Hutton, Irarrazaval & Visvikis (2020). We utilized the implementation proposed in Zhu, Byrd, Lu & Nocedal (1997). We also used the L-BFGS algorithm to minimize $L_1(\cdot, \mathbf{y}_1)$ and $L_2(\cdot, \mathbf{y}_2)$ (without penalty) in order to obtain initial images $\mathbf{x}_1^0$ and $\mathbf{x}_2^0$.

The other part of the alternating scheme is to update the sparse features update $\mathbf{z}_{e,k}^t$ given the current estimate of $\mathbf{x}_e^t$ the sparse features update $\mathbf{z}_{e,k}^t$. This step is achieved using the multi-channel thresholding operator defined in (16).

The pseudo-code for MCAOL reconstruction algorithm is detailed in Algorithm 2.

## 5. 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 (Chun 2019).
## 5.1. 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 (1) with the CAOL regularizer defined in (9) and with the edge-preserving TV regularizers, as well as simultaneously by solving (10) with the JTV regularizer, respectfully defined as

$$R_{tv}(\boldsymbol{x}) = \sum_{j=1}^{J} \sum_{k \in \mathcal{N}_j} \omega_{j,k} \sqrt{(x_j - x_k)^2 + \varepsilon} \quad (27)$$

and

$$R_{jtv}(\boldsymbol{x}_1, \boldsymbol{x}_2) = \sum_{j=1}^{J} \sum_{k \in \mathcal{N}_j} \omega_{j,k} \sqrt{(x_{1,j} - x_{1,k})^2 + (x_{2,j} - x_{2,k})^2 + \varepsilon} \quad (28)$$

where $\mathcal{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 L-BFGS solver to estimate $\boldsymbol{x}_1$ and $\boldsymbol{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 L-BFGS leads to improved qualitative and quantitative results.

## 5.2. Methodology

All experiments were validated by generating the DECT measurements as in Eq. (21) and then running $n = 20$ Poisson noise instances as in Eq. (22) from a ground truth (GT) image $\boldsymbol{x}_e^* = [x_{e,1}^*, \dots, x_{e,J}^*]^\top \in \mathbb{R}^J$, $e = 1, 2$. As performance metrics, we considered the mean absolute bias (AbsBias) error function defined as follows

$$\text{AbsBias} = \frac{1}{N_{\mathcal{R}}} \frac{1}{N_{\text{noise}}} \sum_{j \in \mathcal{R}} \sum_{n=1}^{N_{\text{noise}}} \left| x_{e,j}^{[n]} - x_{e,j}^{*} \right| \quad (29)$$

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

$$\text{STD} = \frac{1}{N_{\mathcal{R}}} \sum_{j \in \mathcal{R}} \sqrt{\frac{1}{N_{\text{noise}}} \sum_{n=1}^{N_{\text{noise}}} (x_{e,j}^{[n]} - \bar{x}_{e,j})^2} \quad (30)$$

where $\bar{x}_{e,j} = \frac{1}{N_{\text{noise}}} \sum_{n=1}^{N_{\text{noise}}} x_{e,j}^{[n]}$. In this work $\mathcal{R}$ corresponds to the nonnegative pixels region of $\boldsymbol{x}_e^*$ 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 (1) and (10) in order to plot AbsBias / STD curves (Fig. 5, 9 and 11). 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.

*Training* The optimization problem (13) is minimized using the BPEG-M algorithm (Chun & Fessler 2019) with normalized input dataset. To investigate the trade-off between accuracy and features sparsity, we tested (13) with different values of $γ_1 = γ_2$ with filter ($d_{1,k}, d_{2,k}$) of dimension $P = 49$ and number of filters $K = 49$. For each simulation, we tuned $γ_1 = γ_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 BPEG-M 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 $α$ in the optimization problem (4) for each separate energy channel.

*Reconstruction* MCAOL and JTV reconstructions were achieved by solving (10) with $R_{mc}$ defined as (18) and (28) respectively, while CAOL and TV reconstructions by solving (1) for each energy bin $e = 1, 2$ separately with $R$ defined as (9) and (27) respectively. MCAOL and CAOL were achieved using $N_{outer} = 300$ outer iterations while the inner image update is obtained using the L-BFGS algorithm with 300 iterations. The $(γ_1, γ_2)$-values and $β$-values were the same as for training. TV and JTV reconstructions were achieved with the L-BFGS algorithm with 300 iterations. The measurements were obtained from the GT images $x_e^*$ outside of the training set and the reconstructions were repeated for each noise instance $n$, for a range of $(ρ_1, ρ_2)$-values with $ρ_1 = ρ_2$ and for a range of $β$-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.

### 5.3. Results on XCAT Phantom

For the unsupervised MCAOL and CAOL training, the numerical data consists of 1-mm pixel-width 512 × 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 (21) and (22)) for reconstruction, as detailed below. We used the MCAOL
weights parameters $\gamma_1 = \gamma_2 = 800$ and the CAOL parameter $\alpha = 0.01$.

<figure><img src="image_3.png" /><figcaption>Figure 2: Learned filters {($d_{1,k}$, $d_{2,k}$)} with $P = K = 49$ using the XCAT training dataset, for (a) MCAOL and (b) CAOL.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: XCAT Phantom: estimated sparse feature maps $z_{e,k}$ for $e = 1, 2$ and $k = 1, \dots, 49$ using CAOL ((a) and (b)) and MCAOL ((c) and (d)); color scale: red for positive values, blue for negative values.</figcaption></figure>

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

In order to generate the sparse-view DECT projection measurements (21), we modeled the projector **A** with a 2-mm full width at half maximum (FWHM) resolution parallel beam system and we used 1-mm pixel-width 406 × 406 GT torso axial-slice images with attenuation coefficients $\mathbf{x}_1^*$, $\mathbf{x}_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° 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. 3 obtained using the XCAT data. By comparing the estimated
<figure><img src="image_5.png" /><figcaption>Figure 4: 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: (a) GT XCAT test image, (b) minimization of the NLL function without prior, (c) CAOL reconstruction, (d) MCAOL joint reconstruction, (e) separate reconstruction using TV prior and (f) joint reconstruction using JTV prior.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Plot of the mean absolute bias (AbsBias) versus the standard deviation (STD) for the XCAT phantom at (a) low X-ray source energy (60 keV) and (b) high X-ray source energy (120 keV).</figcaption></figure>

sparse feature maps $z_{e,k}$ for $e = 1, 2$ and $k = 1, ..., 49$ for separate reconstructions (Fig. 3(a) and Fig. 3(b)) and joint reconstruction (Fig. 3(c) and Fig. 3(d)), there are no similarities between the feature maps obtained from separate reconstructions while the feature maps obtained from joint reconstruction have similar structures.

Fig. 4 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. 5(a) and 5(b). It is worth noting that MCAOL (Fig. 4(d)) manages to substantially reduce the noise as compared with CAOL (Fig. 4(c)).

Fig. 5(a) and Fig. 5(b) 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.

## 5.4. Results on Clinical Data

We utilized images reconstructed from data acquired on Philips IQon Spectral CT and reconstructed with a MBIR technique (Philips IQon Elite Spectral CT product specifications 2018). 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 × 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. 6; we used the parameters $\gamma_1 = \gamma_2 = 10^4$ and the CAOL parameter $\alpha = 10^{-4}$.

<figure><img src="image_7.png" /><figcaption>Figure 6: Learned filters {($d_{1,k}$, $d_{2,k}$)} with $P = K = 49$ using the clinical training dataset, for (a) MCAOL and (b) CAOL.</figcaption></figure>

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

In Fig. 7 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 $z_{e,k}$ for $e = 1, 2$ and $k = 1, ..., 49$ for separate reconstructions (Fig. 7(a) and Fig. 7(b)) and joint reconstruction (Fig. 7(c) and Fig. 7(d)), there are no similarities between the feature maps obtained from separate reconstructions while the feature maps obtained from joint reconstruction have similar structures.

<figure><img src="image_8.png" /><figcaption>Figure 7: Clinical data: estimated sparse feature maps $z_{e,k}$ for $e = 1, 2$ and $k = 1, ..., 49$ using CAOL ((a) and (b)) and MCAOL ((c) and (d)); color scale: red for positive values, blue for negative values.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: 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: (a) GT image, (b) minimization of the NLL function without prior, (c) CAOL reconstruction, (d) MCAOL joint reconstruction, (e) separate reconstruction using TV prior and (f) joint reconstruction using JTV prior.</figcaption></figure>

In Fig. 8 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. 8(d)) is less noisy than the CAOL reconstruction (Fig. 8(c)).

Figs. 9(a) and 9(b) 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
<figure><img src="image_10.png" /><figcaption>Figure 9: Plot of the mean absolute bias (AbsBias) versus the standard deviation (STD) for the clinical data at (a) low X-ray source energy (70 keV) and (b) high X-ray source energy (140 keV).</figcaption></figure>

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.

## 5.5. Results for Low-Dose DECT

We conducted DECT reconstruction on the same set of data as in Section 5.4 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.

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 (24) 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 5.4 and the GT images in Fig. 8(a).

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

Figs. 11(a) and 11(b) 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
<figure><img src="image_11.png" /><figcaption>Figure 10: 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: (a) minimization of the NLL cost function without prior, (b) CAOL reconstruction, (c) MCAOL joint reconstruction, (d) separate reconstruction using TV prior, (e) JTV prior and (f) CAOL-PWLS reconstruction.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Plot of the mean absolute bias (AbsBias) versus the standard deviation (STD) for the low-dose ($I_0$ = 10³) reconstruction with clinical data at (a) low X-ray source energy (70 keV) and (b) high X-ray source energy (140 keV).</figcaption></figure>

by comparing CAOL and CAOL-PWLS.

## 6. 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, \dots, E$ since the $l_{1,0}$ semi-norm in Eq. (15) can be defined for a set of $E$-vectorised feature maps $z_1, \dots, z_E$. Each $z_e, e = 1, \dots, E$ is a column vector of dimension $J \times 1$. Then the joint $l_{1,0}$ semi-norm is defined as

$$ \|(z_1, \dots, z_E)\|_{1,0} = \sum_{j=1}^{J} \mathbf{1}_{[0, +\infty[} (|z_{1,j}| + \dots + |z_{E,j}|) \quad (31) $$

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 Degraux et al. (2017) 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 (10) 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 (Chun & Fessler 2019, Appendix A) 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 (Antun, Renna, Poon, Adcock & Hansen 2020).

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

## 7. 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.

## Acknowledgments

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

## References

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, V., Renna, F., Poon, C., Adcock, B. & Hansen, A. C. (2020). 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, 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, 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.
URL: https://github.com/mechatoz/convolt

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

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

Cueva, E., Meaney, A., Siltanen, S. & Ehrhardt, M. J. (2021). Synergistic multi-spectral ct reconstruction with directional total variation, arXiv preprint arXiv:2101.01834.

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, K., Kamilov, U. S., Boufounos, P. T. & Liu, D. (2017). Online convolutional dictionary learning for multimodal imaging, *2017 IEEE International Conference on Image Processing (ICIP)*, pp. 1617–1621.
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, 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, 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, W., Wu, D., Kim, K., Kalra, M. K., Singh, R., Li, L. & Li, Q. (2021). Iterative material decomposition for spectral ct using self-supervised noise2noise prior, *Physics in Medicine & Biology*.

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, 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, 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, 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, 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, C. H., Leng, S., Yu, L. & Fletcher, J. G. (2015). Dual-and multi-energy CT: principles, technical approaches, and clinical applications, *Radiology* **276**(3): 637–653.

Nocedal, J. & Wright, S. J. (2006). *Numerical Optimization*, 2nd edn, Springer Science & Business Media.

*Philips IQon Elite Spectral CT product specifications* (2018).

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

Ravishankar, S., Nadakuditi, R. R. & Fessler, J. A. (2017). 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, 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, 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, W. P., Mahesh, M., Beck, T. J., Frey, E. C. & Tsui, B. M. W. (2008). Realistic CT simulation using the 4D XCAT phantom, *Medical Physics* **35**(8): 3800–3808.

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, 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, P., Weizman, L., Mota, J. F., Eldar, Y. C. & Rodrigues, M. R. (2019). Coupled dictionary learning for multi-contrast mri reconstruction, *IEEE Transactions on Medical Imaging* **39**(3): 621–633.

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, 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, W., Zhang, Y., Wang, Q., Liu, F., Chen, P. & Yu., H. (2018). Low-dose spectral CT reconstruction
using image gradient $l_0$-norm and tensor dictionary, *Applied Mathematical Modelling* **63**: 538–557.

Xu, L., Lu, C., Xu, Y. & Jia, J. (2011). Image smoothing via $l_0$ gradient minimization, *2011 SIGGRAPH Asia Conference*, pp. 1–12.

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

Yang, Q., Cong, W. & Wang, G. (2017). Superiorization-based multi-energy CT image reconstruction, *Inverse Problems* **33**(4).

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, L., Leng, S. & McCollough, C. H. (2012). Dual-energy CT-based monochromatic imaging, *American Journal of Roentgenology* **199**(5): S9–S15.

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, X., Ravishankar, S., Long, Y. & Fessler, J. A. (2018). 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, 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.