# A Kullback-Leibler approach for 3D reconstruction of spectral CT data corrupted by Poisson noise

Tom Hohweiller, Nicolas Ducros, Françoise Peyrin, Bruno Sixou

► To cite this version:

Tom Hohweiller, Nicolas Ducros, Françoise Peyrin, Bruno Sixou. A Kullback-Leibler approach for 3D reconstruction of spectral CT data corrupted by Poisson noise. SPIE: Optics + Photonics, Aug 2017, San Diego, United States. hal-01557381v2

HAL Id: hal-01557381

https://hal.science/hal-01557381v2

Submitted on 17 Jul 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
# A Kullback-Leibler approach for 3D reconstruction of spectral CT data corrupted by Poisson noise

Tom Hohweiller¹, Nicolas Ducros¹, Françoise Peyrin¹,², and Bruno Sixou¹

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

²The European Synchrotron Radiation Facility, Grenoble, France

## ABSTRACT

While standard computed tomography (CT) data do not depend on energy, spectral computed tomography (SPCT) acquire energy-resolved data, which allows material decomposition of the object of interest. Decompositions in the projection domain allow creating projection mass density (PMD) per materials. From decomposed projections, a tomographic reconstruction creates 3D material density volume. The decomposition is made possible by minimizing a cost function. The variational approach is preferred since this is an ill-posed non-linear inverse problem. Moreover, noise plays a critical role when decomposing data. That is why in this paper, a new data fidelity term is used to take into account of the photonic noise. In this work two data fidelity terms were investigated: a weighted least squares (WLS) term, adapted to Gaussian noise, and the Kullback-Leibler distance (KL), adapted to Poisson noise. A regularized Gauss-Newton algorithm minimizes the cost function iteratively. Both methods decompose materials from a numerical phantom of a mouse. Soft tissues and bones are decomposed in the projection domain; then a tomographic reconstruction creates a 3D material density volume for each material. Comparing relative errors, KL is shown to outperform WLS for low photon counts, in 2D and 3D. This new method could be of particular interest when low-dose acquisitions are performed.

## 1. INTRODUCTION

Spectral computed tomography (SPCT) is a new modality that is expected to be the future of computed tomography (CT). While conventional CT data are obtained integrating the whole energy spectrum received, new detectors, called photo-counting detectors, allow counting how many photons have hit the detectors in several energy bins. Since material attenuation depends on the incident energy of the photon the concentration of the object in a basis of different materials can be estimated from energy-resolved data. Thus SPCT produces more information from a single imaging exam than standard CT. This new modality is expected to have many applications in the medical field, such as imaging high Z contrast agents.¹,²

From 2D energy-resolved projections, three different approaches are available to estimate the 3D volumes of material density for a chosen basis of materials. A first approach consists in performing tomographic reconstruction for the different energy bins, thus creating 3D energy-resolved volumes. From those, materials decomposition can be achieved³ in the object domain. A second approach consists in solving the joint problem of decomposition and tomographic reconstruction, and can be referred as the one-step approach.⁴ The last possibility is to decompose materials in the projection domain,⁵⁻⁷ creating material density map (called projected mass density, PMD) and perform the tomographic reconstruction from the decomposed data to create the 3D volumes of material density. This decomposition approach is highly parallelizable since each can be processed independently of the others. Decomposition in projection domain also embeds all the non-linearities of the forward model. While tomographic reconstruction is linear. For all these reasons, we have selected this approach, thus performing material decomposition in the projection domain followed by the standard filtered-back-projection to create 3D volumes of material density.

Material decomposition in SPCT is a non-linear ill-posed inverse problem. It can be solved by minimizing a cost function, composed of a data fidelity term (representing the distance between the data and the forward

<footnote>Further author information: (Send correspondence to Tom Hohweiller)</footnote>

<footnote>Tom Hohweiller: E-mail: tom.hohweiller@creatis.insa-lyon.fr</footnote>
model) and a regularization term (which make the solution stable regarding the noise). Material decomposition was achieved by minimizing such cost function with different type of algorithms. Both zero-th and first-order algorithms succeeded to that task.⁵,⁷⁻⁹ However, recent work proposed a second-order iterative algorithm (Gauss-Newton scheme) and proved to be more efficient than lower order algorithms.⁶ Ducros et al. introduced a cost function with a weighted least squares (WLS) data fidelity term, that was minimized with a Gauss-Newton algorithm. However, the noise in SPCT is Poisson distributed, and WLS is adapted to Gaussian noise. In this work, we consider a data fidelity term based on the Kullback-Leibler distance to take into account the Poisson nature of the noise.¹⁰ From this change, one can expect to improve the quality of the decomposition.

This paper is organized as follow. In Section 2, the forward model, the inverse problem, and the minimization algorithm are presented. Section 3 presents the KL distance as new data fidelity. Numerical experiments are described in Section 4. Section 5 displays results of experiments and discusses them. Finally, Section 6 concludes this paper.

# 2. MATERIAL DECOMPOSITION IN SPECTRAL CT

## 2.1 Object model

Alvarez and Macovski showed that¹¹ the object of interest can be decomposed on a basis of M materials. This allows to write the linear attenuation coefficient of the object as the following linear combination:

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

where $\mu(E, \mathbf{x})$ is the linear attenuation coefficient (in cm⁻¹) at point $\mathbf{x}$ and energy E, $\rho_m(\mathbf{x})$ is the mass density (in g.cm⁻³) of the m-th material at point $\mathbf{x}$, and $\tau_m(E)$ is the mass attenuation (in cm².g⁻¹) of the m-th material.

## 2.2 Forward model

The standard SPCT forward model⁶,¹² is written as

$$ \mathbf{s}_i(\mathbf{u}) = \int_{\mathbb{R}} d_i(E)n^0(E) \exp \left( -\sum_{m=1}^{M} a_m(\mathbf{u})\tau_m(E) \right) dE \quad (2) $$

where $\mathbf{s}_i(\mathbf{u})$ is the number of photons hitting the detector in the i-th energy bin at pixel $\mathbf{u}$, $d_i(E)$ is the response function of the i-th energy bin of the detector,⁵ $n^0(E)$ is the emitted X-ray spectrum and $\mathbf{a}_m(\mathbf{u})$ is the projection of all the mass of the m-th material along the x-ray path $\mathcal{L}(\mathbf{u})$:

$$ \mathbf{a}_m(\mathbf{u}) = \int_{\mathcal{L}(\mathbf{u})} \rho_m(\mathbf{x}) d\mathbf{x} \quad (3) $$

## 2.3 Discretization and noise

The detector is assumed to have P pixels, with $\mathbf{u}_p$ the location of the p-th pixel, and I energy bins. We define the measurement vector $\mathbf{s} \in \mathbb{R}^{IP}$ as

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

where $s_{i,p}$ is the photon count measured in the i-th energy bin at pixel $\mathbf{u}_p$. Similarly, the projected mass density (PMD) vector $\mathbf{a} \in \mathbb{R}^{MP}$ is defined by

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

where $a_{m,p}$ is PMD of the m-th material at pixel $\mathbf{u}_p$.

Data is assumed to be corrupted by Poisson noise. Let $\mathbf{s}^*$ be the noiseless data vector, $\mathbf{s}$ the noisy data, and $\mathcal{P}(\gamma)$ the Poisson distribution with mean $\gamma$. We have:

$$ \mathbf{s} \sim \mathcal{P}(\gamma = \mathbf{s}^*). \quad (6) $$
## 2.4 Inverse problem

Material decomposition is an ill-posed problem that consists in recovering the PMD $\mathbf{a}$ from the data $\mathbf{s}$. Let $\mathbf{s} = \mathcal{F}(\mathbf{a})$ be the forward model defined in (2). The PMD can be recovered by minimizing the following cost function

$$ \mathcal{C}(\mathbf{s}, \mathcal{F}(\mathbf{a})) = \mathcal{D}(\mathbf{s}, \mathcal{F}(\mathbf{a})) + \alpha \mathcal{R}(\mathbf{a}) \quad (7) $$

where $\mathcal{D}(\mathbf{s}, \mathcal{F}(\mathbf{a}))$ is the fidelity term, $\alpha$ the regularization parameter and $\mathcal{R}(\mathbf{a})$ the regularization term used to stabilize the solution. The regularization terms used was chosen as$^6$

$$ \mathcal{R}(\mathbf{a}) = ||\Delta\mathbf{a}_{\text{soft}}||_2^2 + ||\nabla\mathbf{a}_{\text{bone}}||_2^2 \quad (8) $$

where $\Delta$ and $\nabla$ are first- and second-order differential operator. The usual WLS fidelity term is defined by

$$ \mathcal{D}_{\text{WLS}}(\mathbf{s}, \mathcal{F}(\mathbf{a})) = \frac{1}{2} ||\mathbf{W}(\mathbf{s} - \mathcal{F}(\mathbf{a}))||^2 \quad (9) $$

where $||.||^2$ is the $\ell_2$-norm and $\mathbf{W}$ is the weighting matrix defined by $\mathbf{W} = \text{diag}(\frac{1}{\sqrt{s+1}})$ where $\text{diag}(\mathbf{x})$ represents a diagonal matrix with the elements $\mathbf{x}$.

## 2.5 Minimization of the cost function

In this work, the cost function (7) is minimized by using a Gauss-Newton's algorithm. This second-order algorithm provides faster convergence than first order methods. Starting with an initial guess $\mathbf{a}^0$, the solution is computed iteratively with the following update formula

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

where $\mathbf{a}^k$ is the solution at the iteration $k$, $\delta\mathbf{a}^k$ is the Newton's step and $\lambda^k$ a scalar optimizing the minimization in the direction $\delta\mathbf{a}^k$. The Newton's step is computed by solving

$$ (\mathbf{J}^{\top}\mathbf{Z}_{\text{h}}\mathbf{J} + \alpha\mathbf{H}_{\text{reg}})\delta\mathbf{a}^k = -(\mathbf{J}^{\top}\mathbf{Z}_{\text{g}}(\mathcal{F}(\mathbf{a}^k) - \mathbf{s}) + \alpha\mathbf{g}_{\text{reg}}) \quad (11) $$

where $\mathbf{J}$ is the Jacobian of the forward problem, $\mathbf{H}_{\text{reg}}$ and $\mathbf{g}_{\text{reg}}$ is the Hessian and the gradient of the regularization term, respectively. $\mathbf{Z}_{\text{h}}$ and $\mathbf{Z}_{\text{g}}$ are diagonals matrices which depend on the fidelity term. The step length $\lambda^k$ is obtained by solving

$$ \lambda^k = \underset{\lambda}{\operatorname{argmin}} \mathcal{C}(\mathbf{s}, \mathcal{F}(\mathbf{a}^k + \lambda\delta\mathbf{a}^k)) \quad (12) $$

$\mathbf{Z}_{\text{g}}$ and $\mathbf{Z}_{\text{h}}$ are, in the WLS data fidelity term case, as follow$^{13}$

$$ \text{WLS:} \begin{cases} \mathbf{Z}_{\text{g}} = \mathbf{W}^{\top}\mathbf{W} \\ \mathbf{Z}_{\text{h}} = \mathbf{W}^{\top}\mathbf{W} \end{cases} $$

# 3. A NEW APPROACH: THE KULLBACK-LEIBLER DISTANCE

## 3.1 Noise in spectral CT

The standard data fidelity term WLS is adapted to Gaussian noise. However, in SPCT the noise is Poisson distributed since it is a comes from the detection process of the photons. The WLS is thus not well adapted to our problem for low signal-to-noise ratios. That is why in this work the cost function uses the Kullback-Leibler distance as data fidelity term. The KL term allows taking into account the statistics of the noise more faithfully.$^{10}$
## 3.2 A noise adapted fidelity term

The KL distance between two Poisson-distributed data s and $\mathcal{F}(a)$ is written¹⁴

$$ \mathcal{D}_{\text{KL}}(\mathbf{s}, \mathcal{F}(\mathbf{a})) = \sum_{i} (s_i + \zeta) \log \left( \frac{s_i + \zeta}{\mathcal{F}(a)_i + \zeta} \right) + \mathcal{F}(a)_i - s_i \quad (13) $$

where $\zeta$ is a small integer that avoids dividing by zero and having the log of zero. The Gauss-Newton algorithm requires the computation of the Jacobian and the approximation of the Hessian, as seen in Equation (11). For KL, we have:

$$ \text{KL: } \begin{cases} \mathbf{Z}_g = \mathbf{Y}^{-1} \\ \mathbf{Z}_h = \mathbf{Y}^{-2}\mathbf{S} \end{cases} $$

where $\mathbf{Y} = \text{diag}(\mathcal{F}(\mathbf{a}) + \zeta)$ and $\mathbf{S} = \text{diag}(\mathbf{s} + \zeta)$.

## 4. NUMERICAL SIMULATIONS

Both methods are tested on the numerical phantom DigiMouse.¹⁵ It is a volume of 104 × 190 × 496 voxels composed by $M = 2$ materials: soft tissues and bones. As mentioned earlier, PMD are first decomposed, then, a tomographic reconstruction is performed to retrieve a 3D volume per material. The projection data are obtained with the Radon function of Matlab (219 × 496 pixels). 360 projections angles are used over 360°.

We consider a X-ray spectrum generated by a 90KeV source using the spekCalc software.¹⁶ Three bins of energy ($I = 3$) are used by the photon-counting detectors to detect photons. We define $N^0$ as the number of photons send at each pixel detector:

$$ N^0 = \int_{\mathbb{R}} n^0(E)dE \quad (14) $$

$N^0$ controls the signal-to-noise ratio, thus a large range of $N^0$ values is used: $N^0 \in [10^2, 10^4]$.

The WLS decomposition is performed using the SPRAY Matlab toolbox.¹⁷ The initial guess of the iterative algorithm are uniform PMD maps with $\mathbf{a}_{\text{soft}}^0 = 2 \text{ g.cm}^{-2}$ and $\mathbf{a}_{\text{bones}}^0 = 1 \text{ g.cm}^{-2}$. The algorithm stops when the relative decrease of the cost function is lower than 0.1%, if $\lambda^k$ is lower than $10^{-3}$ or if the maximum number of iteration is reached. This maximum number of iterations is set at 50. A large range of regularization parameter is used in order to find the best decomposition, $\alpha \in [10^{-2}, 10^4]$.

We define the relative decomposition error $\xi_{\theta}$ of the $\theta$-th projection as:

$$ \xi_{\theta} = \sum_{m=1}^{M} \frac{\|\mathbf{a}_{m}^{\text{dec}}(\theta) - \mathbf{a}_{m}^{\text{true}}(\theta)\|}{\|\mathbf{a}_{m}^{\text{true}}(\theta)\|} \quad (15) $$

where $\mathbf{a}_{m}^{\text{dec}}(\theta)$ is the decomposed PMD of the $m^{th}$ material and $\mathbf{a}_{m}^{\text{true}}(\theta)$ is the ground truth for the material $m$. Similarly we define $\delta_Z$ as the relative reconstruction error in the $z$-th slice of the object as:

$$ \delta_Z = \sum_{m=1}^{M} \frac{\|\rho_m^{\text{rec}}(z) - \rho_m^{\text{true}}(z)\|}{\|\rho_m^{\text{true}}(z)\|} \quad (16) $$

where $\rho_m^{\text{true}}(z)$ is the $m$-th material density ground-truth and $\rho_m^{\text{rec}}(z)$ is the reconstructed $m$-th material density.

For each projection angle $\theta$, $\hat{\alpha}_{\theta}$ represents the regularization parameter that minimizes the relative error, i.e.,

$$ \hat{\alpha}_{\theta} = \underset{\alpha}{\operatorname{argmin}} \xi_{\theta}(\alpha) \quad (17) $$

The 3D volumes (154 × 154 × 496 voxels) are reconstructed from the Filtered Back Projection algorithm implemented by the iradon function of Matlab. We used the Hann filter and a v5cubic interpolation method (cubic interpolation from Matlab 5).
<figure><img src="image_2.png" /><figcaption>Figure 1: Evolution of the relative error with the number of photons at 0° (left) and 90° (right)</figcaption></figure>

<table><caption>Table 1: Optimal α value, time per iteration and number of iterations required to converge for α̂ at $N^0 = 10^{2.2}$ for two different projection angle θ</caption><thead><tr><th></th><th>θ(°)</th><th>α̂<sub>θ</sub></th><th>Time (sec)</th><th>Iteration number</th></tr></thead><tbody><tr><td rowspan="2">WLS</td><td>0</td><td>10<sup>2</sup></td><td>3.37</td><td>5</td></tr><tr><td>90</td><td>10<sup>1.5</sup></td><td>3.49</td><td>4</td></tr><tr><td rowspan="2">KL</td><td>0</td><td>10<sup>1.4</sup></td><td>2.95</td><td>4</td></tr><tr><td>90</td><td>10</td><td>3.56</td><td>4</td></tr></tbody></table>

## 5. RESULTS AND DISCUSSION

The PMD are first obtained on the DigiMouse phantom. The computation of the relative errors is very relevant to evaluate the quality of the decompositions.

Figure 1 shows the evolution of the relative error ξ on the projections for the optimized α̂ and for the projection angles θ = 0° (left) and θ = 90° (right) with an increasing number of photons $N^0$. At a low number of photons, the WLS is not a good approximation of the KL distance, which takes account of the photonic noise. For both projection angles, KL outperforms WLS at a small number of photons. At a higher signal-to-noise ratio, the quadratic approximation is good, and both WLS and KL lead to the same relative error. Similar curves can be found by considering each material (soft tissues and bones).

Table 1 provides quantitative results from the decompositions at $N^0 = 10^{2.2}$ photons, for θ = 0° and θ = 90°. This table presents the optimal value of α, the time of computation of one iteration and the number of iteration required to converge. KL and WLS are taking almost the same CPU time (3.5 GHz Xeon E5-1650 v3 and 64 GiB of RAM) to make an iteration and also the same number of iterations to converge. This comes from the calculation of the Jacobian and Hessian where they are almost the same.

Previous results show that KL outperforms WLS at a low number of photons. To highlight this trend, and to validate it, decomposition of soft tissues and bones were done for $θ \in [0, 360]$, at $N^0 = 10^{2.2}$ photons. Evolution of relative error ξ, for α̂, for different θ is depicted on Figure 2 (top). It shows that KL outperforms WLS for any projection angle. On the last two rows decomposition of soft tissues (middle row) and bones (bottom row) are shown for $θ = 180°$. The first column corresponds to the ground truth, the second one to WLS and the last to KL. From those images we can see that KL retrieves more faithfully the PMD from the data. A video is available to see the evolution of the PMD along the angle projection.

A tomographic reconstruction was performed from each material. Figure 3 shows the relative error δ for each slice z. This figure presents the slice 274 since it highlights the best performance of KL with respect to WLS. On soft tissues, KL allows having more details inside the volume than WLS. On bones, KL retrieves more faithfully the value of interest $ρ_{bones}$. Yet, bones and soft tissue regions are mixed in some slices. A video is also available to see the evolution of the error with the slice number.
<figure><img src="image_3.png" /><figcaption>Figure 2: Video1: Ground truth of the material density map of soft tissues and bones in the projection domain (left), decomposition with WLS (middle) and with KL (right) for $N^0 = 10^{2.2}$ photons for $\hat{\alpha}$ at $\theta = 180^\circ$. http://dx.doi.org/doi.number.goes.here</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Video2: Ground truth of the material volume of soft tissues and bones in the object domain (left), decomposition with WLS (middle) and with KL (right) for $N^0 = 10^{2.2}$ photons for $\hat{\alpha}$ at the slice number 274. http://dx.doi.org/doi.number.goes.here</figcaption></figure>
## 6. CONCLUSION

In this work, a new algorithm of material decomposition in the projection domain for SPCT is presented. The cost function uses the KL distance instead of the standard WLS data fidelity term. This allows taking into account the Poisson noise statistics. Since WLS is adapted to Gaussian noise and KL to Poisson noise, the latest divergence should produce better decompositions. Numerical experiments on a mouse phantom, with a large range of signal-to-noise ratios (which controls the approximation of the noise as Gaussian), allow comparing the WLS and KL methods. We obtain PMD maps with the two decomposition methods on the projection domain. Results from experiments show that KL outperforms WLS, considering the relative decomposition error, for any angle of projection $\theta$, at a low number of photons. This efficiency comes from the fact that at a low signal-to-noise ratio the KL distance is more adapted to the noise structure. A tomographic reconstruction is carried out to create a 3D volume of materials. Those volumes are reconstructed from low signal-to-noise projection, where KL is a more adapted data term. In the object domain, KL still has a lower relative error for each slice of the object. In the future, we plan to assess the method with other phantoms, a larger number of materials and to develop methods for automatic selection of the regularization parameter.

## Acknowledgement

This work was supported by the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). It was also supported by the ANR grant ANR-11-INBS-0006.

## REFERENCES

[1] Feuerlein, S., Roessl, E., Proksa, R., Martens, G., Klass, O., Jeltsch, M., Rasche, V., Brambs, H.-J., Hoffmann, M. H., and Schlomka, J.-P., "Multienergy photon-counting k-edge imaging: potential for improved luminal depiction in vascular imaging 1," *Radiology* 249(3), 1010-1016 (2008).

[2] Cormode, D. P., Roessl, E., Thran, A., Skajaa, T., Gordon, R. E., Schlomka, J.-P., Fuster, V., Fisher, E. A., Mulder, W. J., Proksa, R., et al., "Atherosclerotic plaque composition: analysis with multicolor ct and targeted gold nanoparticles 1," *Radiology* 256(3), 774-782 (2010).

[3] Maaß, C., Baer, M., and Kachelrieß, M., "Image-based dual energy ct using optimized precorrection functions: A practical new approach of material decomposition in image domain," *Medical physics* 36(8), 3818-3829 (2009).

[4] Long, Y. and Fessler, J. A., "Multi-material decomposition using statistical image reconstruction for spectral ct," *IEEE transactions on medical imaging* 33(8), 1614-1626 (2014).

[5] Schlomka, J., Roessl, E., Dorscheid, R., Dill, S., Martens, G., Istel, T., Bäumer, C., Herrmann, C., Steadman, R., Zeitler, G., et al., "Experimental feasibility of multi-energy photon-counting k-edge imaging in pre-clinical computed tomography," *Physics in medicine and biology* 53(15), 4031 (2008).

[6] Ducros, N., Abascal, J. F. P. J., Sixou, B., Rit, S., and Peyrin, F., "Regularization of nonlinear decomposition of spectral x-ray projection images," working paper or preprint (Nov. 2016).

[7] Roessl, E. and Herrmann, C., "Cramér-rao lower bound of basis image noise in multiple-energy x-ray imaging," *Physics in medicine and biology* 54(5), 1307 (2009).

[8] Schirra, C. O., Brendel, B., Anastasio, M. A., and Roessl, E., "Spectral ct: a technology primer for contrast agent development," *Contrast media & molecular imaging* 9(1), 62-70 (2014).

[9] Xu, Q., Sawatzky, A., Anastasio, M. A., and Schirra, C. O., "Sparsity-regularized image reconstruction of decomposed k-edge data in spectral ct," *Physics in medicine and biology* 59(10), N65 (2014).

[10] Resmerita, E. and Anderssen, R., "Joint additive kullback-leibler residual minimization and regularization for linear inverse problems," *Mathematical Methods in the Applied Sciences* 30(13), 1527-1544 (2007).

[11] Alvarez, R. E. and Macovski, A., "Energy-selective reconstructions in x-ray computerised tomography," *Physics in medicine and biology* 21(5), 733 (1976).

[12] Roessl, E., Brendel, B., Engel, K.-J., Schlomka, J.-P., Thran, A., and Proksa, R., "Sensitivity of photon-counting based-edge imaging in x-ray computed tomography," *IEEE transactions on medical imaging* 30(9), 1678-1690 (2011).

[13] Favati, P., Lotti, G., Menchi, O., and Romani, F., "Performance analysis of maximum likelihood methods for regularization problems with nonnegativity constraints," *Inverse Problems* 26(8), 085013 (2010).

[14] Kullback, S. and Leibler, R., "On information and sufficiency," *The annals of mathematical statistics* 22(1), 79-86 (1951).

[15] Dogdas, B., Stout, D., Chatziioannou, A., and Leahy, R., "Digimouse: a 3d whole body mouse atlas from CT and cryosection data," *Physics in medicine and biology* 52(3), 577 (2007).

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

[17] "Spray matlab toolbox." https://www.creatis.insa-lyon.fr/~ducros/WebPage/spray.html.