# Regularization of Nonlinear Decomposition of Spectral X-ray Projection Images

Nicolas Ducros, J F P J Abascal, B Sixou, S Rit, Françoise Peyrin

► To cite this version:

Nicolas Ducros, J F P J Abascal, B Sixou, S Rit, Françoise Peyrin. Regularization of Nonlinear Decomposition of Spectral X-ray Projection Images. Medical Physics, 2017, 44 (9), pp.e174-e187. 10.1002/mp.12283 . hal-01391538v3

HAL Id: hal-01391538

https://hal.science/hal-01391538v3

Submitted on 18 May 2017

**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
# Regularization of Nonlinear Decomposition of Spectral X-ray Projection Images

N. Ducros,\* J.F.P.J. Abascal, B. Sixou, S. Rit, and F. Peyrin

Univ Lyon, INSA-Lyon, Université Claude Bernard Lyon 1, UJM-Saint Etienne,
CNRS, Inserm, CREATIS UMR 5220, U1206, F69621, LYON, France

(Dated: May 11, 2017)

**Purpose:** Exploiting the X-ray measurements obtained in different energy bins, spectral computed tomography (CT) has the ability to recover the 3D description of a patient in a material basis. This may be achieved solving two subproblems, namely the material decomposition and the tomographic reconstruction problems. In this work, we address the material decomposition of spectral x-ray projection images, which is a non-linear ill-posed problem.

**Methods:** Our main contribution is to introduce a material-dependent spatial regularization in the projection domain. The decomposition problem is solved iteratively using a Gauss-Newton algorithm that can benefit from fast linear solvers. A Matlab implementation is available online. The proposed regularized weighted least squares Gauss-Newton algorithm (RWLS-GN) is validated on numerical simulations of a thorax phantom made of up to five materials (soft tissue, bone, lung, adipose tissue and gadolinium), which is scanned with a 120 kV source and imaged by a 4-bin photon counting detector. To evaluate the method performance of our algorithm, different scenarios are created by varying the number of incident photons, the concentration of the marker and the configuration of the phantom. The RWLS-GN method is compared to the reference maximum likelihood Nelder-Mead algorithm (ML-NM). The convergence of the proposed method and its dependence on the regularization parameter are also studied.

**Results:** We show that material decomposition is feasible with the proposed method and that it converges in few iterations. Material decomposition with ML-NM was very sensitive to noise, leading to decomposed images highly affected by noise and artifacts even for the best case scenario. The proposed method was less sensitive to noise and improved contrast-to-noise ratio of the gadolinium image. Results were superior to those provided by ML-NM in terms of image quality and decomposition was 70 times faster. For the assessed experiments, material decomposition was possible with the proposed method when the number of incident photons was equal or larger than $10^5$ and when the marker concentration was equal or larger than $0.03 \text{ g.cm}^{-3}$.

**Conclusions:** The proposed method efficiently solves the nonlinear decomposition problem for spectral CT, which opens up new possibilities such as material-specific regularization in the projection domain and a parallelization framework, in which projections are solved in parallel.

## I. INTRODUCTION

Spectral computed tomography (CT) is gaining increasing attention. Exploiting X-ray measurements acquired at multiple energies, this new imaging modality has the ability to recover the concentration maps of the constituents of the tissues in a quantitative manner. While the principle of dual-energy CT has been known for more than 30 years, recent developments in energy selective detectors have boosted the research in this area [1-9]. In particular, it has been shown that photon counting detectors can be used to image high Z contrast agents by exploiting the K-edge energy discontinuity of their linear attenuation coefficient (LAC) [10-12]. Therefore, coupling spectral CT acquisition to the injection of contrast agent is promising not only to increase the sensitivity of standard CT but also to open new clinical applications for X-ray imaging as a functional imaging tool, e.g., for the characterization of the atherosclerotic plaque [13-15].

As in standard CT, spectral CT reconstruction requires to solve an inverse problem. However, spectral CT reconstruction can exploit measurements in different

energy channels, which allows to recover the concentration of different materials [16]. A wide variety of approaches for material reconstruction have been investigated. Spectral CT material reconstruction can be divided into image-based, direct inversion, and projection-based approaches.

The image-based approach is a pragmatic approach that consists in the tomographic reconstruction of each of the energy sinogram followed by a decomposition step. Although it provides valuable results for multi-material decomposition [17-19], it makes substantial assumptions concerning the forward model that may results in image artifacts and quantification inaccuracy. Development of sophisticated correction schemes to alleviate this problem is a topic of ongoing research [20-22].

Direct inversion and projection-based approaches both rely on a mathematical description of the physics of image formation. Direct inversion methods recover the material density directly from the measured energy-resolved sinograms while the projection-based methods first decompose the energy-resolved sinograms into material-specific sinograms and then perform a tomographic reconstruction of each of the decomposed sinograms. Direct inversion is probably the most elegant approach [23-26]. However, the resulting inverse problem

<footnote>\* nicolas.ducros@creastis.insa-lyon.fr</footnote>
is computationally demanding since data from all views have to be processed at the same time. On the contrary, the projection-based approach allows independent decomposition of each projection view, which can be parallelized and is computationally more tractable.

The decomposition step has led to various approaches, some rely on parametric models while the others are based on a physical model. On the one hand, a parametric model avoids measuring the source spectrum and the energy responses of the detector, which may be difficult [27]. Despite its excellent decomposition quality [28–32], this approach is more difficult to implement when three materials or more are to be recovered. This can be a serious drawback for K-edge imaging applications. On the other hand, approaches relying on a physical model of the imaging setup can be more easily applied to arbitrary numbers of materials [11].

In this paper, we present a new framework for the decomposition of a single projection, assuming a physical image formation model is known. The main contribution of this work is to introduce material-specific regularization in the projection domain.

Regularization is not a new idea in spectral CT which is an ill-posed inverse problem. Different regularization strategies have been investigated in order to improve the signal-to-noise ratio of the reconstructions, e.g., statistical penalties and spatial priors. Direct inversion approaches have the ability to incorporate spatial prior in the image domain. Previous works include non-linear edge-preserving statistical reconstruction [23], Gaussian Markov random field regularization with positivity constraints [24], extension of the ART algorithm [25], and TV-regularized primal-dual algorithm [26]. However, it may be noted that these approaches may be computationally expensive in terms of memory requirements and computation time, and have mainly been limited to the reconstruction of 2D slices exploiting a low number of energy bins. Moreover, simplified forward models are often used to moderate the computational cost of the reconstruction. Similarly, projection-based approaches focus on penalties in the image domain which are imposed regularizing the tomographic problem. Edge preserving regularization based on the Huber potential function has been considered in [12, 33] while total variation/ℓ₁ sparsity constrains were used in [34].

However, to the best of our knowledge, regularization of the decomposition step, which is nonlinear and perhaps more difficult to solve than the tomographic step, was only recently addressed [35, 36].

As proposed by Schlomka et al. [11], our forward model is nonlinear and takes into account the spectral response of the detector. In this paper, we also regularize each material projection image. Instead of considering the processing of a single projection line, we consider the decomposition of a 2D projection which is of interest for multi-line and 2D energy resolved detectors. In order to take into account the noise statistics, we have used a weighted least square data term. The inverse

<table><caption>TABLE I. Main mathematical notations</caption><thead><tr><th>Variable</th><th>Meaning</th><th>Dimension</th></tr></thead><tbody><tr><td>a</td><td>Projected mass density (our unknown)</td><td>R(<i>S</i>)</td></tr><tr><td>d</td><td>Detector response function</td><td>R(R, R)</td></tr><tr><td>i</td><td>Output energy bin index</td><td>N</td></tr><tr><td>j</td><td>input energy bin index</td><td>N</td></tr><tr><td>k</td><td>iteration index</td><td>N</td></tr><tr><td>m</td><td>Basis functions index</td><td>N</td></tr><tr><td>n</td><td>Number of input photons</td><td>R(R, <i>S</i>)</td></tr><tr><td>s</td><td>Number of output photons</td><td>R(R, <i>S</i>)</td></tr><tr><td>E</td><td>Input energy</td><td>R</td></tr><tr><td>I</td><td>Number of output energy bins</td><td>1</td></tr><tr><td>J</td><td>Number of input energy bins</td><td>1</td></tr><tr><td>M</td><td>Number of basis functions</td><td>1</td></tr><tr><td>P</td><td>Number of pixels</td><td>1</td></tr><tr><td>μ</td><td>Linear attenuation</td><td>R(R, Ω)</td></tr><tr><td>ρ</td><td>Material mass density</td><td>R(Ω)</td></tr><tr><td>τ</td><td>Material basis function</td><td>R(R)</td></tr><tr><td>Ω</td><td>Object volume</td><td>R<sup>3</sup></td></tr><tr><td>u</td><td>Detector pixel location</td><td>R<sup>2</sup></td></tr><tr><td>x</td><td>Object voxel location</td><td>R<sup>3</sup></td></tr><tr><td>a</td><td>Projected mass density image</td><td>R<sup>PM</sup></td></tr><tr><td>s</td><td>Photon counts image</td><td>R<sup>IP</sup></td></tr><tr><td>H</td><td>Hessian of <i>R</i></td><td>R<sup>MP×MP</sup></td></tr><tr><td>J</td><td>Jacobian of <i>F</i></td><td>R<sup>IP×MP</sup></td></tr><tr><td>ε</td><td>Output energy</td><td>R</td></tr><tr><td>R</td><td>Regularizer</td><td>R(R<sup>PM</sup>)</td></tr><tr><td>S</td><td>Detection surface</td><td>R<sup>2</sup></td></tr><tr><td>F</td><td>Forward model</td><td>R<sup>IP</sup>(R<sup>PM</sup>)</td></tr></tbody></table>

problem is solved using a Gauss-Newton algorithm that exploits second-order derivatives, which converges faster than the most commonly used first- or zeroth-order algorithms. Gauss-Newton minimization was successfully applied to several medical imaging inverse problems [37–40]. In this paper, it is implemented using explicit computation of the Hessian, which is a sparse and structured matrix. Hence, our algorithm can benefit from efficient linear equation solvers. The formulation is presented in a general framework, providing analytical formulas for general regularization functionals. Our algorithm, which is referred to as regularized weighted least square Gauss-Newton (RWLS-GN), is validated on synthetic data and the projection of a voxelized thorax phantom based on a clinical CT volume. We provide results for different scenarios created by varying the number of incident photons and the amount of contrast agent. To promote reproducibility of this research, the computer code (MATLAB scripts) related to this study is available online [41].

## II. THEORY

As depicted in figure 1, we consider a 3-dimensional (3-D) object in a bounded domain Ω that is imaged by a 2-D detector with a surface S. A summary of the notation used throughout the paper is provided in table I.
<figure><img src="image_2.png" /><figcaption>FIG. 1. Principle of acquisition and notations. Here, a conical geometry is chosen as an illustration. Ω is the object volume, S is the detector surface, ρm is the mass density of the mth material, and am is projected mass density of the mth material.</figcaption></figure>

### A. X-ray/matter interaction

Let $n(E, \mathbf{u})$ denotes the number of photons of energy $E$ that reaches the detector at the pixel position $\mathbf{u} \in S$. Neglecting scattering within the object, $n$ follows the Beer-Lambert law. Mathematically,

$$n(E, \mathbf{u}) = n_0(E, \mathbf{u}) \exp \left[ - \int_{\mathcal{L}(\mathbf{u})} \mu(E, \mathbf{x}) \, d\ell \right] \quad (1)$$

where $n_0(E, \mathbf{u})$ is the source spectrum, $\mathcal{L}(\mathbf{u})$ is the acquisition line that depends on the source emission geometry (parallel, cone-beam, etc), and $\mu(E, \mathbf{x})$ is the local LAC of the object at energy $E$ and point $\mathbf{x} \in \Omega$. For simplicity, the source spectrum is assumed to be the same at each pixel, i.e.,

$$n_0(E, \mathbf{u}) = n_0(E), \quad \forall \mathbf{u} \in S \quad (2)$$

### B. Detection model

The mean signal recorded by a photon counting detector may be modelled by

$$\bar{s}(\mathcal{E}, \mathbf{u}) = \int_{\mathbb{R}} d(\mathcal{E}, E) n(E, \mathbf{u}) \, dE \quad (3)$$

where $d(\mathcal{E}, E)$ accounts for the detector response function, $E$ being the incident photon energy and $\mathcal{E}$ the detected energy. The photons detected within the ith energy bin [$\mathcal{E}_i, \mathcal{E}_{i+1}$] are accumulated electronically

thanks to a counting circuit. The mean number of photons detected within the ith energy bin is given by

$$\bar{s}_i(\mathbf{u}) = \int_{\mathbb{R}} d_i(E) n(E, \mathbf{u}) \, dE \quad (4)$$

where

$$d_i(E) = \int_{\mathcal{E}_i}^{\mathcal{E}_{i+1}} d(\mathcal{E}, E) \, d\mathcal{E} \quad (5)$$

is the response function of the ith bin of the detector. Note that the detection model assumes that the value measured at a given pixel is not correlated to the value measured at neighbor pixels, i.e. charge sharing can be neglected or is corrected for [42].

### C. Object decomposition

It is assumed that the LAC can be decomposed as the sum of M basis functions that are separable in energy and space. We have

$$\mu(E, \mathbf{x}) = \sum_{m=1}^{M} \rho_m(\mathbf{x}) \tau_m(E), \quad \forall \mathbf{x} \in \Omega \quad (6)$$

where the $\tau_m$ are some well-chosen basis functions and the $\rho_m$ are the corresponding weights in the decomposition. Following the work of Alvarez and Macovski [16], two approaches have emerged for the choice of the basis functions: i) physics-based where $\tau_m$ models the physical effects, e.g. photoelectric, Compton scattering, and, ii) material-based where $\tau_m$ is the mass attenuation of the constituents of the objects (in cm$^2$.g$^{-1}$). In the latter approach, $\rho_m$ is the mass density (in g.cm$^{-3}$) of material m. For the sake of simplicity but with no loss of generality, $\rho_m$ is referred to as mass density in the following. K-edge materials may be added with the two types of basis [11].

### D. Forward model

Substituting (6) into (1) and the obtained result into (4), the mean measured signal may be written

$$\bar{s}_i(\mathbf{u}) = \int_{\mathbb{R}} n_0(E) d_i(E) \exp \left[ - \sum_{m=1}^{M} a_m(\mathbf{u}) \tau_m(E) \right] dE \quad (7)$$

where

$$a_m(\mathbf{u}) = \int_{\mathcal{L}(\mathbf{u})} \rho_m(\mathbf{x}) \, d\ell \quad (8)$$

is the projection of the mass density $\rho_m$ along the line $\mathcal{L}(\mathbf{u})$ (see figure 1) that is expressed in g.cm$^{-2}$ when a material-based decomposition is performed. In the following, $a_m$ is referred to as the projected mass density (PMD).
## E. Discretization

Let us assume that the detector consists of an array of $P$ pixels located at $\mathbf{u}_p$, $p \in \{1, \dots, P\}$. Each of the pixels counts photons in $I$ different energy bins. Let $\bar{s} \in \mathbb{R}^{PI}$ and $\mathbf{s} \in \mathbb{R}^{PI}$ be the mean and measured photon counts vector, which are defined by

$$ \bar{s} = [\bar{s}_{1,1} \dots \bar{s}_{I,1} \dots \dots \bar{s}_{1,P} \dots \bar{s}_{I,P}]^\top \quad (9) $$

$$ \mathbf{s} = [s_{1,1} \dots s_{I,1} \dots \dots s_{1,P} \dots s_{I,P}]^\top. \quad (10) $$

In a similar way, let $\mathbf{a} \in \mathbb{R}^{PM}$ be the (unknown) vector containing the PMD of each material in each pixel, which is defined by

$$ \mathbf{a} = [a_{1,1} \dots a_{M,1} \dots \dots a_{1,P} \dots a_{M,P}]^\top. \quad (11) $$

The discrete forward model can be formalized as

$$ \bar{s} = \mathcal{F}(\mathbf{a}) \quad (12) $$

where $\mathcal{F}(\mathbf{a})$ denotes the non-linear mapping induced by (7).

The material decomposition problem is to recover the PMD vector $\mathbf{a}$ from the measured data $\mathbf{s}$ given the known forward mapping $\mathcal{F}$. This is an ill-posed problem that requires prior knowledge about the solution to stabilize the inversion in the presence of noise.

## F. Noise model

The measurements are assumed to be corrupted by independent Poisson noise, i.e.,

$$ s_{i,p} \sim \mathbf{P}(\lambda = \bar{s}_{i,p}), \quad (13) $$

where $\mathbf{P}(\lambda)$ denotes the Poisson distribution of mean $\lambda$.

In order to investigate the influence of noise which is related to the value of $\bar{s}_{i,p}$, we introduce the total number of photons emitted by the source during the acquisition. It is defined by

$$ N_0 = \int n_0(E) dE \quad (14) $$

The larger $N_0$ is, the better the signal-to-noise ratio of the measurements.

# III. PROPOSED DECOMPOSITION METHOD

## A. Variational approach

In this manuscript, we propose to recover the projected mass vector within a variational framework, i.e., minimizing the cost function

$$ \mathcal{C}(\mathbf{a}) = \frac{1}{2} \|\mathbf{s} - \mathcal{F}(\mathbf{a})\|_{\mathbf{W}}^2 + \alpha \mathcal{R}(\mathbf{a}) \quad (15) $$

<figure><img src="image_3.png" /><figcaption>FIG. 2. Potential functions used in this manuscript. The quadratic potential function leads to the $\ell_2$, while the Huber function approximates the $\ell_1$-norm. The Huber function is depicted for $\epsilon = \{0.5, 0.2, 0.05, 0.02\}$. Note that $\psi^{\text{huber}}(b_q)$ behaves as $\frac{1}{2} \frac{b_q^2}{\epsilon}$ for $b_q \ll \epsilon$ and as $|b_q| - \epsilon$ for $b_q \gg \epsilon$</figcaption></figure>

where $\mathbf{W}$ is a weighting matrix, $\mathcal{R}$ is the regularization functional, and $\alpha$ is a global regularization parameter. The weighting matrix is chosen as

$$ \mathbf{W} = \mathrm{diag}\left(\frac{1}{\mathbf{s}}\right) \quad (16) $$

such that $\|\mathbf{s} - \mathcal{F}(\mathbf{a})\|_{\mathbf{W}}^2$ is a weighted least squares fidelity term adapted to data corrupted by Poisson noise.

One main difficulty of the inverse problem is the non linearity of the forward mapping, which is addressed by means of the Gauss-Newton algorithm described in Section III C.

## B. Regularization functional

Priors must be suited to the images of materials to be recovered. For instance, a soft tissue image is expected to be smoother than a bone image. We adopt the following regularizing functional, which allows the prior of the different materials to be tuned independently:

$$ \mathcal{R}(\mathbf{a}) = \sum_{m=1}^{M} \beta_m \mathcal{R}_m(\mathbf{b}_m), \quad \text{with} \quad \mathbf{b}_m = \mathbf{L}_m \mathbf{a}_m, \quad (17) $$

where $\beta_m$, $\mathcal{R}_m$, $\mathbf{a}_m = [a_{m,1} \dots a_{m,P}]^\top \in \mathbb{R}^P$, and $\mathbf{L}_m \in \mathbb{R}^{Q_m \times P}$ are the regularization parameter, the regularization functional, the PMD, and a linear transform of the $m$-th basis function (e.g., the derivative or the laplacian), respectively.

Many different regularizing functionals, together with many different transforms have been considered in the past. In this manuscript, we limit our study to the following cases.
a. Material regularization functional We consider the functionals that can be modeled by

$$ \mathcal{R}_m(\mathbf{b}_m) = \sum_{q=1}^{Q_m} \psi(b_q) \quad (18) $$

where $\psi: \mathbb{R} \to \mathbb{R}_+$ is a twice differentiable potential function chosen to promote some image prior and $b_q = (\mathbf{b}_m)_q$ is the q-th component of $\mathbf{b}_m$. In particular, the following two potential functions, illustrated in figure 2, are considered:

* Quadratic. The quadratic potential function

$$ \psi^{\text{quad}}(b_q) = b_q^2 \quad (19) $$

can be chosen to promote smooth solutions. It leads to the standard $\ell_2$ penalization also known as Tikhonov regularization:

$$ \mathcal{R}_m^{\text{quad}}(\mathbf{b}_m) = \|\mathbf{b}_m\|_2^2 \quad (20) $$

* Huber. The Huber penalty functional can be chosen to preserve edges in the solution. In particular, its smooth approximation is defined by

$$ \psi_m^{\text{huber}}(b_q) = \sqrt{b_q^2 + \epsilon^2} - \epsilon. \quad (21) $$

where $\epsilon$ is a parameter that can be tuned knowing the order of magnitude of the $b_q$'s. This choice for $\psi$ results in a $\ell_1$-like penalization term:

$$ \mathcal{R}_m^{\text{huber}}(\mathbf{b}_m) \approx \|\mathbf{b}_m\|_1 \quad (22) $$

b. Linear transform The choice for the linear operator $\mathbf{L}_m$ is restricted to discrete differential operators up to second order, i.e.,

$$ \mathbf{L}_m \in \{\mathbf{D}^{(0)}, \mathbf{D}^{(1)}, \mathbf{D}^{(2)}\} \quad (23) $$

where $\mathbf{D}^{(0)} \in \mathbb{R}^{P \times P}$ is the identity matrix, $\mathbf{D}^{(1)} \in \mathbb{R}^{2P \times P}$ is the 2-D discrete gradient, and $\mathbf{D}^{(2)} \in \mathbb{R}^{P \times P}$ is the 2-D discrete Laplacian. Any other linear operators could be used as well. However, these choices lead to sparse $\mathbf{L}_m$ matrices, which constitute a computational advantage as detailed in Section VI.

### C. Optimisation algorithm

We propose to minimize (15) using Gauss-Newton (GN)'s method, which is a classical iterative method for non-linear minimization. It starts with an initial guess $\mathbf{a}^{(0)}$ and builds new estimates with the following update rule:

$$ \mathbf{a}^{(k+1)} = \mathbf{a}^{(k)} + \lambda^{(k)} \delta\mathbf{a}^{(k)}. \quad (24) $$

where $\delta\mathbf{a}^{(k)} \in \mathbb{R}^{PM}$ is the so-called GN step and $\lambda^{(k)} \in \mathbb{R}$ is the step length.

The GN step is obtained solving

$$ (\mathbf{J}^{(k)})^{\top} \mathbf{W} \mathbf{J}^{(k)} + \alpha \mathbf{H}^{(k)} \delta\mathbf{a}^{(k)} = -\mathbf{g}^{(k)}, \quad (25) $$

where $\mathbf{J}^{(k)}$ is the Jacobian matrix of $\mathcal{F}$ at the iterate $\mathbf{a}^{(k)}$, $\mathbf{H}^{(k)}$ is the Hessian matrix of $\mathcal{R}$ at $\mathbf{a}^{(k)}$, and $\mathbf{g}^{(k)}$ is the gradient of $\mathcal{C}$ at $\mathbf{a}^{(k)}$ that is given by:

$$ \mathbf{g}^{(k)} = -\mathbf{J}^{(k)}^{\top} \mathbf{W}^{\frac{1}{2}} \left[ \mathbf{s} - \mathcal{F}(\mathbf{a}^{(k)}) \right] + \alpha \nabla \mathcal{R}(\mathbf{a}^{(k)}). \quad (26) $$

The step length is obtained minimizing the cost function along the direction provided by the GN direction. Mathematically, it is defined as

$$ \lambda^{(k)} = \arg \min \mathcal{C}(\mathbf{a}^{(k)} + \lambda \delta\mathbf{a}^{(k)}) \quad (27) $$

### D. Computation of the Jacobian

The Jacobian matrix of the forward model introduced in (25) is the $IP \times MP$ real-valued matrix that maps variations in the unknown vector $\mathbf{a}$ onto variations in the forward model $\bar{\mathbf{s}} = \mathcal{F}(\mathbf{a})$. It is defined by

$$ \mathbf{J} = \frac{\partial \bar{\mathbf{s}}}{\partial \mathbf{a}} \quad (28) $$

Differentiating the forward model given in (7), it can be shown that

$$ \frac{\partial \bar{s}_{i,p}}{\partial a_{m,p'}} = \begin{cases} 0, & \text{for } p \neq p' \\ j_{i,m,p}, & \text{for } p = p' \end{cases} \quad (29) $$

with

$$ j_{i,m,p} = - \int_{\mathbb{R}} \tau_m(E)n(E, \mathbf{u}_p)d_i(E) dE. \quad (30) $$

As a result, the Jacobian can be expressed as the block diagonal matrix

$$ \mathbf{J} = \text{diag}([\mathbf{J}_1 \dots \mathbf{J}_P]) \quad (31) $$

where $\mathbf{J}_p \in \mathbb{R}^{I \times M}$ is such that $(\mathbf{J}_p)_{i,m} = j_{i,m,p}$. It is noteworthy that $\mathbf{J}$ is a sparse matrix having only $IMP$ non-zero entries out of $IMP^2$. Physically, the many zeros of $\mathbf{J}$ indicate that a change of the PMD in a particular pixel does not affect the value measured at another pixel.

### E. Computation of the gradient of the regularizer

Computation of the gradient of the regularizer $\mathcal{R}$ is required by (26). By linearity of the gradient, we obtain from the definition of the regularizer (17) that

$$ \nabla \mathcal{R}(\mathbf{a}) = \sum_{m=1}^{M} \beta_m \nabla \mathcal{R}_m(\mathbf{L}_m \mathbf{a}_m), \quad (32) $$
where $\nabla$ is the gradient with respect to $\mathbf{a}$, which is defined by

$$ \nabla = \left[ \frac{\partial}{\partial a_{1,1}} \cdots \frac{\partial}{\partial a_{M,1}} \cdots \frac{\partial}{\partial a_{1,P}} \cdots \frac{\partial}{\partial a_{M,P}} \right]^{\top}. \quad (33) $$

Since regularization is performed on each material independently, i.e. $\mathcal{R}_m$ only depends on $\mathbf{a}_m$, we have $\frac{\partial \mathcal{R}_m}{\partial a_{m',p}} = 0, \forall m' \neq m$. It is therefore convenient to introduce the partial gradient w.r.t. the PMD of the m-th material

$$ \nabla_m = \left[ \frac{\partial}{\partial a_{m,1}} \cdots \frac{\partial}{\partial a_{m,P}} \right]^{\top}. \quad (34) $$

and to remark that the full gradient $\nabla \mathcal{R}_m \in \mathbb{R}^{MP}$ is simply obtained by inserting zeros into the partial gradients $\nabla_m \mathcal{R}_m \in \mathbb{R}^P$. Precisely, $\nabla \mathcal{R}_m = \nabla_m \mathcal{R}_m \otimes \mathbf{e}_m$, where $\mathbf{e}_m \in \mathbb{R}^M$ is the natural basis vector given by $\mathbf{e}_m = [\delta_{1,m} \cdots \delta_{M,m}]^\top$ where $\delta_{m',m}$ is the Kronecker delta. Therefore, we have

$$ \nabla \mathcal{R}(\mathbf{a}) = \sum_{m=1}^{M} \beta_m \nabla_m \mathcal{R}_m(\mathbf{L}_m \mathbf{a}_m) \otimes \mathbf{e}_m. \quad (35) $$

The latter equation is of particular interest when closed-form formulas for $\nabla_m \mathcal{R}_m$ are available. Applying the appendix formula (A2) to the potential functions considered in (19) and (21), the following simple identities are obtained:

$$ \nabla_m \mathcal{R}_m^{\text{quad}}(\mathbf{b}_m) = 2 \mathbf{L}_m^{\top} \mathbf{b}_m \quad (36) $$

and

$$ \nabla_m \mathcal{R}_m^{\text{huber}}(\mathbf{b}_m) = \mathbf{L}_m^{\top} \left[ \frac{b_1}{\sqrt{b_1^2 + \epsilon^2}} \cdots \frac{b_Q}{\sqrt{b_Q^2 + \epsilon^2}} \right]^{\top} \quad (37) $$

where $b_q = (\mathbf{L}_m \mathbf{a}_m)_q$.

## F. Computation of the Hessian of the regularizer

The Hessian of the regularizer is required to compute the GN update that solves (25). It is defined by:

$$ \mathbf{H} = \frac{\partial^2 \mathcal{R}}{\partial \mathbf{a} \partial \mathbf{a}^\top}. \quad (38) $$

With the same reasoning as in the previous section, we have

$$ \mathbf{H}(\mathbf{a}) = \sum_m \beta_m \mathbf{H}_m(\mathbf{b}_m) \otimes \mathbf{e}_m \otimes \mathbf{e}_m^\top. \quad (39) $$

where $\mathbf{H}_m \in \mathbb{R}^{P \times P}$ is a shorthand for the Hessian of $\mathcal{R}_m$ w.r.t. $\mathbf{b}_m$. At this point it is important to notice that $\mathbf{H}$ is built from the $M$ matrices $\mathbf{H}_m$, each of size $P \times P$.

Moreover, closed-form equations for $\mathbf{H}_m$ are obtained applying the appendix equation (B2) to the potential functions given in (19) and (21). We have:

$$ \mathbf{H}_m^{\text{quad}}(\mathbf{b}_m) = 2 \mathbf{L}_m^{\top} \mathbf{L}_m \quad (40) $$

<figure><img src="image_4.png" /><figcaption>FIG. 3. Top row: Spectrum of the x-ray source. Middle row: detector response for each energy bin. Bottom row: mass attenuation of the three materials considered in the numerical experiments (gadolinium in green, bone in grey, lung tissue in red, soft tissue in yellow, adipose tissue in cyan; lung and soft tissues are superimposed).</figcaption></figure>

and

$$ \mathbf{H}_m^{\text{huber}}(\mathbf{b}_m) = \epsilon^2 \mathbf{L}_m^{\top} \text{diag} \left( \left[ (b_1^2 + \epsilon^2)^{-\frac{3}{2}} \cdots (b_Q^2 + \epsilon^2)^{-\frac{3}{2}} \right] \right) \mathbf{L}_m \quad (41) $$

where $b_q = (\mathbf{L}_m \mathbf{a}_m)_q$.

# IV. NUMERICAL SIMULATIONS AND IMAGE ANALYSIS

## A. Acquisition parameters

The source spectrum $n_0(E)$ was simulated using the spekCalc software [43, 44], considering a tube voltage of 120 kV, a 12° anode angle, and a filter of 1.2 mm of Al. It is depicted in the top row of figure 3. Measurements were performed in $I = 4$ energy bins. The response function $d_i(E)$ of each of the 4 bins was computed according to (5), where $d(\mathcal{E}, E)$ was simulated using the parametric model provided by Schlomka et al. [11, (see appendix A.2)] and choosing the following energy thresholds: $\mathcal{E}_1 = 15$ keV, $\mathcal{E}_2 = 36$ keV, $\mathcal{E}_3 = 60$ keV, and $\mathcal{E}_4 = 91$ keV (see figure 3). In order to investigate the robustness of our decomposition method to noise, different total numbers of photons were considered in the range $N_0 \in [10^4, 10^7]$ photons.

The projections were simulated with RTK [45] using a point x-ray source at 1.5 m from a 41 × 41 cm² flat panel
with 256 × 256 pixels. The ray casting used bi-linear interpolation in each slice of the phantom following and 16 × 16 rays were calculated and averaged per pixel.

## B. Phantom generation

We used the thoracic-abdominal CT scan depicted in figure 4a), which is part of the 3D-IRCADb database [46]. It is a 480 × 370 × 167 volume with a 0.961 × 0.961 × 1.8 mm³ voxel size and an 8-bit grayscale. We also used the anatomical atlas represented in figure 4d). It was obtained and provided by Kéchichian et al. [47] who segmented the previously mentioned CT scan via vicinity prior. Given both the CT and atlas volumes, the two material mass density volumes depicted in figure 4b) and c) were created as follows. First, segments were grouped. For the three-material phantom, all the segments corresponding to soft tissue[48] were grouped into the same segment Ωsoft. For the five material phantom, Ωsoft did not include adipose tissue[49] and lung tissue[50] that were assigned to the segments Ωfat and Ωlung, respectively. For both phantoms, the bone segment Ωbone remained unchanged. Second, the mass density in each voxel was estimated from the CT images. A linear mapping between the 8-bit grayscale of the CT scan and the mass density was assumed, i.e., ρm = λμCT. The constant λ (in g.cm⁻³) was calibrated such that mean value of the mass density of visceral fat is set to 0.95 g.cm⁻³. Third, the portal vein segment Ωportal was marked a uniform concentration cgd of gadolinium, which is a promising spectral CT contrast agent [13, 33]. The previous three steps result in the computation of three volumes that were each projected onto the detector according to (8) and the acquisition geometry L(u) described in Section IV A. The resulting projected material mass density image is displayed in figure 4d).

The mass attenuations of soft tissue, bone, and gadolinium, as required in (7), were taken from the ICRU report 44 [51]. They are depicted in figure 3 (bottom row).

## C. Image quality metrics

The quality of the decomposed images has been assessed by means of two different objective performance metrics. First, the normalized euclidean distance has been considered. Given the ground-truth image a<sub>m</sub><sup>true</sup> and the decomposed image a<sub>m</sub><sup>dec</sup>, the normalized Euclidean error for the m-th material is defined as

$$ error_m = \frac{\|a_m^{\text{dec}} - a_m^{\text{true}}\|_2}{\|a_m^{\text{true}}\|_2}. \qquad (42) $$

The lower the normalized Euclidean error, the better. We also define the total normalized Euclidean error as

$$ error_{\text{tot}} = \frac{1}{M} \sum_{m=1}^{M} error_{m}. \qquad (43) $$

Second, the contrast-to-noise ratio (CNR) has been considered. Given a region of interest S<sub>roi</sub> ∈ S and a background region S<sub>back</sub> = S \ S<sub>roi</sub>, we define the CNR of the decomposed image a<sub>m</sub><sup>dec</sup> as

$$ cnr_m = \frac{\bar{a}_{\text{roi}} - \bar{a}_{\text{back}}}{(w_{\text{roi}}\sigma_{\text{roi}}^2 + w_{\text{back}}\sigma_{\text{back}}^2)^{1/2}}, \qquad (44) $$

where $\bar{a}_{\text{roi}}$ ($\bar{a}_{\text{back}}$) and $\sigma_{\text{roi}}$ ($\sigma_{\text{back}}$) are the average and standard deviation of a<sub>m</sub><sup>dec</sup> in the region-of-interest S<sub>roi</sub> (background S<sub>back</sub>) and w<sub>roi</sub> (w<sub>back</sub>) is the partial surface of S<sub>roi</sub> (S<sub>back</sub>) in the image. The larger the CNR, the better.

## D. Algorithm parameters

The quality of the reconstructed images depends dramatically on the choice of the regularization parameters α introduced in (15). In order to fairly compare decompositions obtained from different acquisition scenarios, each decomposition was performed for a set of regularization parameters in the range log(α) ∈ [-2, 2], with a 0.5 step. The best regularization parameter α* has been chosen to obtain the lowest reconstruction error, for the phantom with three materials, and then the same regularization parameter has been used to assess the method for different phantom configurations. For the selection of the material specific regularization parameters β<sub>soft</sub>, β<sub>bone</sub>, β<sub>gd</sub>, we evaluated the algorithm for values between 0.1 and 10 and found that the simple choice β<sub>soft</sub> = β<sub>bone</sub> = β<sub>gd</sub> = 1 provided results close to the best ones (see supplemental data in table V at the end of the manuscript). The hyperparameter ϵ was set to a small value ϵ = 0.01.

Our RWLS-GN algorithm has two stopping criteria. It is interrupted either if i) the length step is small enough, i.e., λ(k) < λ<sub>stop</sub> or ii) the cost function decrease is small enough, i.e., 1 - C(k)/C(k-1) < δ<sub>stop</sub>. We chose λ<sub>stop</sub> = 5 × 10⁻³ and δ<sub>stop</sub> = 10⁻³. The initial guess a<sub>0</sub><sup>0</sup> was chosen as constant images with values in the range of the expected values for each material, i.e., a<sub>0</sub><sup>0</sup><sub>soft</sub> = 10 g.cm⁻², a<sub>0</sub><sup>0</sup><sub>bone</sub> = 1 g.cm⁻², and a<sub>0</sub><sup>0</sup><sub>gd</sub> = 0 g.cm⁻².

# V. RESULTS

## A. Comparison between ML-NM and RWLS-GN methods

Our algorithm was compared to the reference method described in [11], which relies on the maximization of
<figure><img src="image_5.png" /><figcaption>FIG. 4. Numerical phantom. (a) Slice of the 3D-IRCADb [46] thoracic-abdominal CT scan (in Hounsfield units). (b) Slice of the 3-material phantom showing the mass density ρ<sub>m</sub> (g.cm<sup>-3</sup>) of soft tissue in yellow, bone in gray, and gadolinium in green. (c) Slice of the 5-material phantom showing the mass density ρ<sub>m</sub> (g.cm<sup>-3</sup>) of soft tissue in yellow, bone in gray, gadolinium in green, lung tissue in red, and adipose tissue in blue. (d) Anatomical atlas obtained segmenting the 3D-IRCADb volume via vicinity prior [47], different colors indicate different labels. (e) Projection of the 3-material phantom in the coronal plane. The projected material mass density ρ<sub>m</sub> (g.cm<sup>-2</sup>) is displayed with the color code used in (b). (f) Projection of the 5-material phantom in the coronal plane. Same color code as in (c).</figcaption></figure>

a likelihood data fidelity term by means of the Nelder-Mead method. It is referred to as ML-NM and implemented using the `fminsearch` Matlab function, setting the maximum number of iterations to 2000 × M.

Ground truth images of the PMD of soft tissue, bone and gadolinium are shown in the top row of figure 5. For the best case scenario (number of incident photons of 10⁷ and concentration of the marker 1 g.cm⁻³), decomposed images provided by both the RWLS-GN and the ML-NM methods are given in the middle row of figure 5. Difference between decomposition and ground truth images is shown in the bottom row of figure 5. While materials were well separated when using RWLS-GN, the images decomposed using ML-NM have poor quality. Large artifacts are observed in the soft and bone tissue images and the gadolinium image is very noisy. For ML-NM, the largest errors are related to the recovery of gadolinium, which results in a dramatic cross-talk with the soft and bone tissue images. For RWLS-GN, it can be seen that the gadolinium image presents the lowest error and that the largest error occurred at borders of the spine. Quantitatively, RWLS-GN led to error<sub>soft</sub> = 0.014, error<sub>bone</sub> = 0.271, error<sub>gd</sub> = 0.071, and cnr<sub>gd</sub> = 3.42 for α* = 10<sup>-0.5</sup>, while ML-NM led to error<sub>soft</sub> = 1430, error<sub>bone</sub> = 1558, error<sub>gd</sub> = 15, and cnr<sub>gd</sub> = 0.3627.

For the considered 256 × 256 images, a RWLS-GN ite-
ration requires about 2s on a standard laptop (2.6 GHz

i7 CPU and 16 GiB of RAM) which leads a decomposition time of about 20s assuming 10 iterations before convergence as observed in all the assessed experiments. The ML-NM is much slower with a decomposition time of about 1500s for the same images.

B. Convergence of the proposed RWLS-GN
algorithm

Convergence of the algorithm is shown in figure 6. The method converged in around ten iterations. The step length was above one in almost all iterations and the decrease of the cost function was close to one in all except for the last iterations. Number of iterations for convergence was between 10 and 15 for all scenarios considered in this work.

C. Influence of the regularization parameter

The influence of the regularization parameter α is il-
lustrated in figure 7. Figure 7a shows the solution error
for the different materials as a function of α. The lowest
error for all materials was attained for similar values of α
in the range [0.1, 1]. Figure 7b and 7c shows underregu-
larized and overregularized images, respectively. When
<figure><img src="image_6.png" /><figcaption>FIG. 5. Typical material decompositions. The projected mass densities (g.cm⁻²) are obtained considering the 3-material phantom shown in figure 4b/4e. The gadolinium concentration $c_{gd}$ was set to 1 g.cm⁻³ and the number of incident photons $N_0$ to $10^7$ photons. (a) Ground truth images. (b) Images recovered by our GN-RWLS method for $\alpha^* = 10^{-0.5}$ (best regularization parameter). (c) Images recovered by the maximum likelihood Nelder-Mead (ML-NM) method described in [11]. (d) GN-RWLS error. (e) ML-NM error.</figcaption></figure>

regularization is not sufficient, images are noisy and the method fails to decompose materials. On the contrary, excessive regularization leads to oversmoothened images, failing to properly decompose soft and bone tissue. Overall, the gadolinium image presented the best image quality and was more robust to the selection of the regularization parameter, and separation of bone and soft tissue were more sensitive to the selection of $\alpha$.

The optimal regularization parameter for all scenarios is shown in table II. For a sufficiently large number of photons ($N_0 \ge 10^5$) the optimal value $\log(\alpha^*)$ was found to be almost independent of the number of photons and concentration of the marker ($\log(\alpha^*) = 0$ or $-0.5$, see table II). For a low number of photons ($N_0 \le 10^5$) the method converged to unsatisfying solutions. In addition, it was found that values of $\alpha$ in the range 0.1 and 1 led to similar results in terms of image decomposition and recovery of the concentration of the marker.

For the rest of the analysis in this work, the regulari-

<table><caption>TABLE II. Best regularization parameter. The table displays $\log(\alpha^*)$</caption><thead><tr><th>$N_0$</th><th>1</th><th>0.3</th><th>0.1</th><th>0.03</th><th>0.01</th></tr></thead><tbody><tr><td>$10^7$</td><td>-0.5</td><td>-0.5</td><td>-0.5</td><td>0.0</td><td>0.0</td></tr><tr><td>$10^6$</td><td>0.0</td><td>-0.5</td><td>0.0</td><td>0.0</td><td>0.0</td></tr><tr><td>$10^5$</td><td>-0.5</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td></tr></tbody></table>

zation parameter $\alpha$ has been set to $10^{-0.5}$.

### D. Influence of noise

The influence of noise has been assessed by varying the number of incident photons from $N_0 = 10^7$ to $N_0 = 10^4$ for a concentration of the marker of $c_{gd} = 1$ g.cm⁻³.
<figure><img src="image_7.png" /><figcaption>FIG. 6. Convergence of the algorithm for $\alpha = 10^{-0.5}$, $c_{gd} = 1$ g.cm$^{-3}$, and $N_0 = 10^7$ photons. (Top left) cost function $C$ as a function of the iteration number $k$. (Top right) total normalized Euclidean error $error_{tot}$ as a function of the iteration number $k$. (Bottom left) step length $\lambda^{(k)}$ vs iteration number $k$. (Bottom right) cost function decrease $\frac{c^{(k-1)}-c^{(k)}}{c^{(k-1)}}$ vs iteration number. The red dashed lines indicate the values of the chosen stopping criteria.</figcaption></figure>

<table><caption>TABLE III. Total normalized Euclidean distance $error_{tot}$ for different scenarios</caption><thead><tr><th>$N_0$</th><th>1</th><th>0.3</th><th>0.1</th><th>0.03</th><th>0.01</th></tr></thead><tbody><tr><td>10<sup>7</sup></td><td>0.12</td><td>0.13</td><td>0.21</td><td>0.49</td><td>1.26</td></tr><tr><td>10<sup>6</sup></td><td>0.19</td><td>0.21</td><td>0.33</td><td>0.75</td><td>1.97</td></tr><tr><td>10<sup>5</sup></td><td>0.33</td><td>0.39</td><td>0.61</td><td>1.45</td><td>4.00</td></tr></tbody></table>

For $N_0 = 10^7$ materials were well separated and quantification of the marker concentration was correct (see figure 5). Lowering photons from $N_0 = 10^6$ to $N_0 = 10^5$ decreased image quality for bone and soft tissue while the gadolinium image was less affected (see top and middle rows in figure 8). For a number of photons $N_0 = 10^4$, material decomposition was not feasible.

These results are confirmed by the Euclidean distance error (first column in table III) and contrast-to-noise ratio of the marker (first column in table IV). Lowering the number of photons from $N_0 = 10^7$ to $N_0 = 10^5$ increased the error between two and three times but has limited influence on the contrast-to-noise ratio of the marker.

<figure><img src="image_8.png" /><figcaption>FIG. 7. Influence of regularization. Measurements were simulated for $N_0 = 10^7$ photons and $c_{gd} = 1$ g.cm$^{-3}$. (a) Normalized euclidean distance as a function of the regularization parameter $\alpha$. (b) Underregularized images recovered for $\alpha = 10^{-4}$. (c) Overregularized images recovered for $\alpha = 10^2$.</figcaption></figure>

### E. Influence of the contrast agent marker

The influence of the contrast agent marker was studied by varying the marker concentration $c_{gd}$ from 1 to 0.01. Decomposed gadolinium images for different marker concentration and number of photons are shown in figure 9. Marker concentrations of 1 g.cm$^{-3}$ were well recovered even for a low number of photons ($N_0 = 10^5$). Lowering the marker concentration to 0.1 had a large effect on image quality even for a large number of photons, which is indicated by an increased decomposition error. For concentrations below 0.01 g.cm$^{-3}$, recovery of the
<figure><img src="image_9.png" /><figcaption>FIG. 8. Influence of noise. Decomposed images for decreasing number of incident photons: $N_0 = 10^6$ photons (top), $N_0 = 10^5$ photons (middle), and $N_0 = 10^4$ photons (bottom). The marker concentration was fixed to $1 \text{ g.cm}^{-3}$.</figcaption></figure>

<table><caption>TABLE IV. Contrast-to-noise ratio of the decomposed marker image.</caption><thead><tr><th>N<sub>0</sub></th><th>1</th><th>0.3</th><th>0.1</th><th>0.03</th><th>0.01</th></tr></thead><tbody><tr><td>10<sup>7</sup></td><td>3.42</td><td>3.41</td><td>3.05</td><td>1.86</td><td>0.82</td></tr><tr><td>10<sup>6</sup></td><td>3.48</td><td>3.34</td><td>2.73</td><td>1.31</td><td>0.56</td></tr><tr><td>10<sup>5</sup></td><td>3.66</td><td>2.91</td><td>1.82</td><td>0.65</td><td>0.24</td></tr></tbody></table>

marker was not feasible.

Contrast-to-noise of the decomposed marker shows that the recovery of the marker was more affected by marker concentration than by noise (see table IV). Contrast-to-noise was similar to the best case scenario for $c_{gd} = 1$ when $N_0 \ge 10^5$, $c_{gd} = 0.3$ when $N_0 \ge 10^6$ and for $c_{gd} = 0.1$ when $N_0 \ge 10^7$.

## F. Influence of the phantom

Robustness of the algorithm with respect to the phantom is shown in figure 10. Results correspond to different phantom compositions (soft tissue, bone tissue, and gadolinium only or together with lung and adipose tissue) and projection angles. In all cases, three materials (soft tissue, bone tissue and gadolinium) are considered for decomposition. Decomposition of the sagittal view is found to be feasible considering the algorithm parameters used to decompose the coronal view. For a phantom with five materials (soft tissue, bone, lung tissue, adipose

<figure><img src="image_10.png" /><figcaption>FIG. 9. Decomposed gadolinium images obtained for decreasing marker concentrations $c_{gd} = \{1, 0.1, 0.01\} \text{ g.cm}^{-3}$ (from left to right) and decreasing number of incident photons $N_0 \in \{10^7, 10^6, 10^5\}$ from top to bottom. In all cases $\alpha = \alpha^*$.</figcaption></figure>

tissue and gadolinium) decomposed images are visually similar to those from a phantom with only three materials (compare figure 5b and 10a and figure 10b and 10c). While the bone error is found to be higher in a five-material phantom, it is interesting to note that the gadolinium error is of the same order in both phantom (about 0.07 for the coronal view and 0.2 for the sagittal view).

## VI. DISCUSSION AND CONCLUSION

This work presented a new and efficient method to solve the nonlinear decomposition problem in spectral X-ray imaging. The problem was solved using a Gaussian Newton step that exploited sparsity of Jacobian and Hessian operators.

It has been validated on numerical experiments of a thorax phantom and has been compared to the ML-NM algorithm. We showed that the proposed RWLS-GN method solves the nonlinear decomposition problem in an efficient manner. The algorithm converged in around ten iterations for all scenarios. It was also shown to improve the accuracy of the projected material maps. Decomposed material images presented high image quality, low error and accurate concentration of the marker. On the contrary, ML-NM decomposition was highly affected by noise, leading to artifacts in bone and soft tissue images and to a very noisy gadolinium image. For the assessed experiments, RWLS-GN decomposition was feasible for a
<figure><img src="image_11.png" /><figcaption>FIG. 10. Robustness of the decomposition setting $\alpha = 10^{-0.5}$. From left to right, recovered soft tissue, bone, and gadolinium images. (a) Decomposition of the 5-material phantom (coronal view) leading to error_bone = 0.697, error_gd = 0.066, and $\alpha^* = 10^{-0.5}$. (b) Decomposition of the 3-material phantom (sagittal view) leading to error_soft = 0.019, error_bone = 0.330, error_gd = 0.212, and $\alpha^* = 10^{-1}$. (c) Decomposition of the 5-material phantom (sagittal view) leading to error_bone = 0.701, error_gd = 0.201, and $\alpha^* = 10^{-1}$. In all cases, measurements were simulated for $N_0 = 10^7$ photons and $c_{gd} = 1 \text{ g.cm}^{-3}$.</figcaption></figure>

number of incident photons larger or equal than $N_0 = 10^5$ and for a marker concentration above $0.03 \text{ g.cm}^{-3}$. In the context of spectral CT, our method could exploit data redundancy in the projection domain and be parallelized across projections.

Most previous studies in spectral CT included regularization in the tomographic problem [12, 33, 34]. These methods decompose the data in a previous step and then perform the tomographic problem by regularizing in the image domain, which is a simpler but computationally more expensive problem as it requires to project and backproject the data at each iteration. In contrast, in this work we tackle regularizing the nonlinear decomposition problem, which has not been widely addressed. This problem is nonlinear, harder to solve and it presents

some differences but it also offers some advantages that can be exploited. The two main advantages are that the problem is smaller (i.e. one projection at a time) and that the forward operator (and also gradient and Hessian operators) is very sparse. For instance the number of nonzero elements in the Hessian is of the order of the number of pixels instead of the square of the number of pixels. The proposed algorithm exploits these features by using sparse matrices to explicitly represent gradient and Hessian operator and a Gauss-Newton step, which has a higher convergence rate than the most widely used linearized methods. At each iteration, the linear system (25) can be efficiently solved using a Cholesky decomposition. Fast solvers for the factorization of sparse matrix are available. Their complexity depends on the sparsity (percentage of non-zero entries) of the matrix and also the sparsity pattern (location of the non zero entries). These solvers are perfectly adapted to the resolution of (25) since the matrix $\hat{\mathbf{H}} = 2\mathbf{J}^\top\mathbf{J} + \alpha\mathbf{H}$ is both sparse and structured. Indeed, it follows from (31) that $\mathbf{J}^\top\mathbf{J} = \text{diag}([\mathbf{J}_1^\top\mathbf{J}_1 \dots \mathbf{J}_P^\top\mathbf{J}_P])$, i.e., $\mathbf{J}^\top\mathbf{J}$ is block diagonal with a sparsity of $M^2P/M^2P^2 = 1/P$. Considering $P = 256 \times 256$ pixels, the sparsity (portion of non-zero values) of $\mathbf{J}^\top\mathbf{J}$ is 0.0015 %. Moreover, the matrix $\mathbf{H}$ is built from the products of sparse banded matrices $\mathbf{H}_m$, hence it is also sparse banded.

This work has focused on the decomposition of spectral X-ray images acquired under a particular projection angle and may find direct applications in interventional radiology. In this context, it could avoid the need of image subtraction that is subject to motion artifacts. In the context of spectral CT, its main advantage is that it can be parallelized, processing each projection in parallel. This has the potential to provide a fast method for sinogram decomposition. Although it has been shown that RWLS-GN decomposition is feasible for different projection views and that it improves standard ML-NM decomposition, its extension to spectral CT requires further investigation.

One of the main differences with previous studies in the literature is that most of the approaches regularize in the spatial domain and very few works have addressed regularization in the projection domain. The proposed algorithm included regularization on the projection domain and allowed the use of specific regularization functional for each material. The problem was solved using a Gauss-Newton step with efficient computation of operators. The selection of the regularization parameter was studied and the improvement brought by regularization was highlighted for different experimental scenarios. The gain of including regularization in the material decomposition step was shown by comparing the proposed algorithm to ML-NM, which led to unsatisfactory decomposed images highly affected by noise and artifacts.

The proposed method allows the use of different regularization functionals for each material and it is formulated in a general framework for any choice of the operator. In this study, we selected Total Variation for the
contrast material and first-order Tikhonov and second-order Tikhonov for bone and soft-tissue, respectively. These choices are based on assumptions that gadolinium projection images should be piecewise constant and projection image for soft tissue is smooth. Nevertheless, other choices of regularization functionals could be more suitable and provide better results, and should be investigated in future work.

## ACKNOWLEDGMENTS

This work was supported by the grant ANR-11-INBS-0006 of the French National Research Agency (ANR) and was performed within the framework of the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program ANR-11-IDEX-0007. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement number 701915.

The authors have no relevant conflicts of interest to disclose.

## Appendix A: Gradient of the regularizer

We seek an expression for the gradient of the regularizer $\mathcal{R}(\mathbf{b})$, $\mathbf{b} = \mathbf{La}$, with respect to $\mathbf{a}$. Applying the chain rule leads to:

$$ \nabla_{\mathbf{a}} \mathcal{R}(\mathbf{b}) = \mathbf{L}^{\top} \nabla_{\mathbf{b}} \mathcal{R}(\mathbf{b}), \quad (A1) $$

where $\nabla_{\mathbf{a}}$ is the gradient w.r.t. the variable $\mathbf{a}$ and $\nabla_{\mathbf{b}}$ is the gradient w.r.t. $\mathbf{b}$. Differentiating $\mathcal{R}(\mathbf{b}) = \sum_q \psi(b_q)$ with respect to $\mathbf{b}$, we obtain

$$ \nabla_{\mathbf{a}} \mathcal{R}(\mathbf{b}) = \mathbf{L}^{\top} [\psi'(b_1) \dots \psi'(b_Q)]^{\top}. \quad (A2) $$

## Appendix B: Hessian of the regularizer

Applying the chain rule twice, it may be shown that the Hessian of $\mathcal{R}(\mathbf{b})$ with respect to the variable $\mathbf{a}$, for a linear change of variable $\mathbf{b} = \mathbf{La}$, is given by (see equation 2.1 of Lebl [52])

$$ \mathbf{H}_{\mathbf{a}}[\mathcal{R}](\mathbf{b}) = \mathbf{L}^{\top} \mathbf{H}_{\mathbf{b}}[\mathcal{R}](\mathbf{b}) \mathbf{L} \quad (B1) $$

where $\mathbf{H}_{\mathbf{b}}[\mathcal{R}]$ is the Hessian of $\mathcal{R}(\mathbf{b})$ w.r.t. the variable $\mathbf{b}$. Since $\mathcal{R}(\mathbf{b}) = \sum_q \psi(b_q)$ we have

$$ \mathbf{H}_{\mathbf{a}}[\mathcal{R}](\mathbf{b}) = \mathbf{L}^{\top} \mathrm{diag} ([\psi''(b_1) \dots \psi''(b_Q)]) \mathbf{L}. \quad (B2) $$

[1] Polad M Shikhaliev, “Energy-resolved computed tomography: first experimental results,” Physics in Medicine and Biology 53, 5595 (2008).

[2] Michael Overdick, Christian Baumer, K.J. Engel, Johannes Fink, Christoph Herrmann, H. Kruger, Matthias Simon, R. Steadman, and G. Zeitler, “Status of direct conversion detectors for medical imaging with x-rays,” Nuclear Science, IEEE Transactions on 56, 1800–1809 (2009).

[3] F. Cassol Brunner, M. Dupont, C. Meessen, Y. Boursier, H. Ouamara, A. Bonissent, C. Kronland-Martinet, J. Clemens, F. Debarbieux, and C. Morel, “First k-edge imaging with a micro-CT based on the xpad3 hybrid pixel detector,” Nuclear Science, IEEE Transactions on 60, 103–108 (2013).

[4] R. Ballabriga, M. Campbell, E. Heijne, X. Llopart, L. Tlustos, and W. Wong, “Medipix3: A 64k pixel detector readout chip working in single photon counting mode with improved spectrometric performance,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 633, Supplement 1, S15–S18 (2011), 11th International Workshop on Radiation Imaging Detectors (IWORID).

[5] E.N. Gimenez, R. Ballabriga, M. Campbell, I. Horswell, X. Llopart, J. Marchal, Kawal J S Sawhney, N. Tartoni, and D. Turecek, “Characterization of medipix3 with synchrotron radiation,” Nuclear Science, IEEE Transactions on 58, 323–332 (2011).

[6] Cheng Xu, M. Persson, Han Chen, S. Karlsson, M. Danielsson, C. Svensson, and H. Bornefalk, “Evaluation of a second-generation ultra-fast energy-resolved asic for photon-counting spectral CT,” Nuclear Science, IEEE Transactions on 60, 437–445 (2013).

[7] X Wang, D Meier, S Mikkelsen, G E Maehlum, D J Wagenaar, B M W Tsui, B E Patt, and E C Frey, “MicroCT with energy-resolved photon-counting detectors,” Physics in Medicine and Biology 56, 2791 (2011).

[8] Xiaolan Wang, Dirk Meier, Katsuyuki Taguchi, Douglas J. Wagenaar, Bradley E. Patt, and Eric C. Frey, “Material separation in x-ray CT with energy resolved photon-counting detectors,” Medical Physics 38, 1534–1546 (2011).

[9] Katsuyuki Taguchi and Jan S. Iwanczyk, “Vision 20/20: Single photon counting x-ray detectors in medical imaging,” Medical Physics 40, 100901 (2013).

[10] E Roessl and R Proksa, “K-edge imaging in x-ray computed tomography using multi-bin photon counting detectors,” Physics in Medicine and Biology 52, 4679 (2007).

[11] J P Schlomka, E Roessl, R Dorscheid, S Dill, G Martens, T Istel, C Bumer, C Herrmann, R Steadman, G Zeitler, A Livne, and R Proksa, “Experimental feasibility of multi-energy photon-counting k-edge imaging in pre-clinical computed tomography,” Physics in Medicine and Biology 53, 4031 (2008).
[12] C.O. Schirra, E. Roessl, T. Koehler, B. Brendel, A. Thran, D. Pan, M.A. Anastasio, and R. Proksa, “Statistical reconstruction of material decomposed data in spectral CT,” Medical Imaging, IEEE Transactions on 32, 1249–1257 (2013).

[13] Sebastian Feuerlein, Ewald Roessl, Roland Proksa, Gerhard Martens, Oliver Klass, Martin Jeltsch, Volker Rasche, Hans-Juergen Brambs, Martin H. K. Hoffmann, and Jens-Peter Schlomka, “Multienergy photon-counting k-edge imaging: Potential for improved luminal depiction in vascular imaging,” Radiology 249, 1010–1016 (2008), http://radiology.rsna.org/content/249/3/1010.full.pdf+html.[27]

[14] David P. Cormode, Ewald Roessl, Axel Thran, Torjus Skajaa, Ronald E. Gordon, Jens-Peter Schlomka, Valentin Fuster, Edward A. Fisher, Willem J. M. Mulder, Roland Proksa, and Zahi A. Fayad, “Atherosclerotic plaque composition: Analysis with multicolor CT and targeted gold nanoparticles,” Radiology 256, 774–782 (2010), http://radiology.rsna.org/content/256/3/774.full.pdf+html.

[15] Dipanjan Pan, Ewald Roessl, Jens-Peter Schlomka, Shelton D. Caruthers, Angana Senpan, Mike J. Scott, John S. Allen, Huiying Zhang, Grace Hu, Patrick J. Gaffney, Eric T. Choi, Volker Rasche, Samuel A. Wickline, Roland Proksa, and Gregory M. Lanza, “Computed tomography in color: Nanok-enhanced spectral ct molecular imaging,” Angewandte Chemie International Edition 49, 9635–9639 (2010).

[16] R E Alvarez and A Macovski, “Energy-selective reconstructions in x-ray computerised tomography,” Physics in Medicine and Biology 21, 733 (1976).

[17] Huy Q. Le and Sabee Molloi, “Least squares parameter estimation methods for material decomposition with energy discriminating detectors,” Med. Phys. 38, 245 (2011).

[18] Darin P Clark, Ketan Ghaghada, Everett J Moding, David G Kirsch, and Cristian T Badea, “In vivo characterization of tumor vasculature using iodine and gold nanoparticles and dual energy micro-CT,” Physics in Medicine and Biology 58, 1683 (2013).

[19] Sebastian Faby, Stefan Kuchenbecker, Stefan Sawall, David Simons, Heinz-Peter Schlemmer, Michael Lell, and Marc Kachelrieß, “Performance of today’s dual energy CT and future multi energy CT in virtual non-contrast imaging and in iodine quantification: A simulation study,” Medical Physics 42, 4349–4366 (2015).

[20] Mats Persson, Ben Huber, Staffan Karlsson, Xuejin Liu, Han Chen, Cheng Xu, Moa Yveborg, Hans Bornefalk, and Mats Danielsson, “Energy-resolved CT imaging with a photon-counting silicon-strip detector,” Physics in Medicine and Biology 59, 6709 (2014).

[21] P.R.S. Mendonca, P. Lamb, and D.V. Sahani, “A flexible method for multi-material decomposition of dual-energy CT images,” Medical Imaging, IEEE Transactions on 33, 99–116 (2014).

[22] Mengheng Touch, Darin P Clark, William Barber, and Cristian T Badea, “A neural network-based method for spectral distortion correction in photon counting x-ray CT,” Physics in Medicine and Biology 61, 6132 (2016).

[23] Y. Long and J. A. Fessler, “Multi-material decomposition using statistical image reconstruction for spectral CT,” IEEE Transactions on Medical Imaging 33, 1614–1626 (2014).

[24] Ruoqiao Zhang, J.-B. Thibault, C.A. Bouman, K.D. Sauer, and Jiang Hsieh, “Model-based iterative recon-

struction for dual-energy x-ray CT using a joint quadratic likelihood model,” Medical Imaging, IEEE Transactions on 33, 117–134 (2014).

[25] Yunsong Zhao, Xing Zhao, and Peng Zhang, “An extended algebraic reconstruction technique (E-ART) for dual spectral CT,” IEEE transactions on medical imaging 34, 761–768 (2015).

[26] Rina Foygel Barber, Emil Y Sidky, Taly Gilat Schmidt, and Xiaochuan Pan, “An algorithm for constrained one-step inversion of spectral CT data,” Physics in Medicine and Biology 61, 3784 (2016).

[27] E.Y. Sidky, L. Yu, X. Pan, Y. Zou, and M. Vannier, “A robust method of x-ray source spectrum estimation from transmission measurements: Demonstrated on computer simulated, scatter-free transmission data,” Journal of applied physics 97, 124701 (2005).

[28] Robert E. Alvarez, “Estimator for photon counting energy selective x-ray imaging with multibin pulse height analysis,” Medical Physics 38, 2324–2334 (2011).

[29] N. Maaß, S. Sawall, M. Knaup, and M. Kachelrieß, “Empirical multiple energy calibration (emec) for material-selective CT,” in 2011 IEEE Nuclear Science Symposium Conference Record (2011) pp. 4222–4229.

[30] Kevin C Zimmerman and Taly Gilat Schmidt, “Experimental comparison of empirical material decomposition methods for spectral CT,” Physics in Medicine and Biology 60, 3175 (2015).

[31] R. E. Alvarez, “Efficient, non-iterative estimator for imaging contrast agents with spectral x-ray detectors,” IEEE Transactions on Medical Imaging 35, 1138–1146 (2016).

[32] Dufan Wu, Li Zhang, Xiaohua Zhu, Xiaofei Xu, and Sen Wang, “A weighted polynomial based material decomposition method for spectral x-ray CT imaging,” Physics in Medicine and Biology 61, 3749–3783 (2016).

[33] A. Sawatzky, Q. Xu, C. O. Schirra, and M. A. Anastasio, “Proximal admm for multi-channel image reconstruction in spectral x-ray CT,” (2014).

[34] Qiaofeng Xu, Alex Sawatzky, Mark A Anastasio, and Carsten O Schirra, “Sparsity-regularized image reconstruction of decomposed k-edge data in spectral CT,” Physics in Medicine and Biology 59, N65 (2014).

[35] Bernhard Brendel, Frank Bergner, Kevin Brown, and Thomas Koehle, “Penalized likelihood decomposition for dual layer spectral CT,” in *Proceeding of the fourth International Meeting in X-ray Computed Tomography* (2016).

[36] Nicolas Ducros, Simon Rit, Bruno Sixou, and Françoise Peyrin, “Non-linear regularized decomposition of spectral x-ray projection images,” in *Proceeding of the fourth International Meeting in X-ray Computed Tomography* (2016).

[37] M. Vauhkonen, D. Vadasz, P. A. Karjalainen, E. Somersalo, and J. P. Kaipio, “Tikhonov regularization and prior information in electrical impedance tomography,” IEEE Transactions on Medical Imaging 17, 285–293 (1998).

[38] Martin Schweiger, Simon R Arridge, and Ilkka Nissil, “Gauss-Newton method for image reconstruction in diffuse optical tomography,” Physics in Medicine and Biology 50, 2365 (2005).

[39] J. F. P.-J. Abascal, J. Chamorro-Servent, J. Aguirre, S. Arridge, T. Correia, J. Ripoll, J. J. Vaquero, and M. Desco, “Fluorescence diffuse optical tomography using the split bregman method,” Medical Physics 38, 6275–6284 (2011).
[40] Nicolas Ducros, Teresa Correia, Andrea Bassi, Gianluca Valentini, Simon Arridge, and Cosimo D'Andrea, "Reconstruction of an optical inhomogeneity map improves fluorescence diffuse optical tomography," Biomedical Physics & Engineering Express 2, 055020 (2016).

[41] "https://www.creatis.insa-lyon.fr/~ducros/WebPage/spectral_ct.html".

[42] Thomas Koenig, Marcus Zuber, Elias Hamann, Angelica Cecilia, Rafael Ballabriga, Michael Campbell, Marie Ruat, Lukas Tlustos, Alex Fauler, Michael Fiederle, and Tilo Baumbach, "How spectroscopic x-ray imaging benefits from inter-pixel communication," Physics in Medicine and Biology 59, 6195-6213 (2014).

[43] Gavin G. Poludniowski, "Calculation of x-ray spectra emerging from an x-ray tube. part ii. x-ray production and filtration in x-ray targets," Medical Physics 34, 2175-2186 (2007).

[44] G Poludniowski, G Landry, F DeBlois, P M Evans, and F Verhaegen, "Spekcalc : a program to calculate photon spectra from tungsten anode x-ray tubes," Physics in Medicine and Biology 54, N433 (2009).

[45] S. Rit, M. Vila Oliva, S. Brousmiche, R. Labarbe, D. Sarrut, and G.C. Sharp, "The reconstruction toolkit (RTK), an open-source cone-beam CT reconstruction toolkit based on the insight toolkit (ITK)," J. Phys.: Conf. Ser. 489, 012079 (2014).

[46] "3d image reconstruction for comparison of algorithm database," .

[47] R. Kéchichian, S. Valette, M. Desvignes, and R. Prost, "Shortest-path constraints for 3d multiobject semiautomatic segmentation via clustering and graph cut," IEEE Transactions on Image Processing 22, 4224-4236 (2013).

[48] Namely, os, heart, lung (R), lung (L), liver, spleen, gallbladder, pancreas, kidney (R), kidney (L), adrenal gland (R), adrenal gland (L), hyperplasia, aorta, vena cava, portal vein, skin, fat, visceral fat, muscle, marrow, renal pyramids (R), renal pyramids (L), digestive tract tissue, digestive tract cavity, pulmonary arteries (R), pulmonary arteries (L), bronchi(R), bronchi (L), intervertebral discs, hepatic arteries, hepatic vena cava, hepatic portal veins.

[49] Namely, fat and visceral fat.

[50] Namely, bronchi(R), bronchi (L).

[51] International Commission on Radiation Units and Measurements (ICRU), *Report 44 of the Tissue Substitutes in Radiation Dosimetry and Measurement*, Tech. Rep. (Bethesda, MD, 1989).

[52] Jiri Lebl, *Tasty Bits of Several Complex Variables*, second edition ed., 9781365095573 (Copyright Jiri Lebl, 2016).
<table><caption>TABLE V. (Supplemental data). Quality of the decomposition for different hyperparameters {$\beta_{soft}, \beta_{bone}, \beta_{gd}$} and the best global regularization parameter $\alpha^*$. In the simulations, we set $N_0 = 10^7$ photons and $c_{gd} = 1 \text{ g.cm}^{-3}$.</caption><thead><tr><th>β<sub>soft</sub></th><th>β<sub>bone</sub></th><th>β<sub>gd</sub></th><th>log(α<sup>*</sup>)</th><th>error</th><th>cnr<sub>gd</sub></th></tr></thead><tbody><tr><td>1</td><td>1</td><td>1</td><td>-0.5</td><td>0.12</td><td>3.42</td></tr><tr><td>1</td><td>0.1</td><td>0.1</td><td>0</td><td>0.12</td><td>3.42</td></tr><tr><td>1</td><td>0.1</td><td>10</td><td>-1</td><td>0.17</td><td>3.64</td></tr><tr><td>1</td><td>10</td><td>0.1</td><td>-1</td><td>0.14</td><td>3.33</td></tr><tr><td>1</td><td>10</td><td>10</td><td>-1.5</td><td>0.13</td><td>3.40</td></tr><tr><td>0.1</td><td>1</td><td>0.1</td><td>0</td><td>0.14</td><td>3.35</td></tr><tr><td>0.1</td><td>1</td><td>10</td><td>-1</td><td>0.15</td><td>3.68</td></tr><tr><td>10</td><td>1</td><td>0.1</td><td>-1</td><td>0.12</td><td>3.41</td></tr><tr><td>10</td><td>1</td><td>10</td><td>-1</td><td>0.12</td><td>3.52</td></tr><tr><td>0.1</td><td>0.1</td><td>1</td><td>0</td><td>0.13</td><td>3.57</td></tr><tr><td>0.1</td><td>10</td><td>1</td><td>-0.5</td><td>0.16</td><td>3.42</td></tr><tr><td>10</td><td>0.1</td><td>1</td><td>-0.5</td><td>0.13</td><td>3.45</td></tr><tr><td>10</td><td>10</td><td>1</td><td>-1.5</td><td>0.12</td><td>3.39</td></tr></tbody></table>