# Spectral CT material decomposition in the presence of poisson noise: A Kullback-Leibler approach

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

► To cite this version:

Tom Hohweiller, Nicolas Ducros, Françoise Peyrin, Bruno Sixou. Spectral CT material decomposition in the presence of poisson noise: A Kullback-Leibler approach. Innovation and Research in BioMedical engineering, 2017, 38 (4), pp.214-218. 10.1016/j.irbm.2017.06.002 . hal-01557367v2

HAL Id: hal-01557367

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

Submitted on 15 Nov 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
SPECTRAL CT MATERIAL DECOMPOSITION IN THE PRESENCE
OF POISSON NOISE: A KULLBACK-LEIBLER APPROACH

Tom Hohweiller<sup>a</sup>, Nicolas Ducros<sup>a</sup>, Françoise Peyrin<sup>a,b</sup>, Bruno Sixou<sup>a</sup>

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

<sup>b</sup>The European Synchrotron Radiation Facility, Grenoble, France

**Abstract**

**Context:** Spectral Computed Tomography (SPCT) acquires energy-resolved projections that can be used to evaluate the concentrations of the different materials of an object. This can be achieved by decomposing the projection images in a material basis, in a first step, followed by a standard tomographic reconstruction step. However, material decomposition is a challenging non-linear ill-posed inverse problem, which is very sensitive to the photonic noise corrupting the data.

**Methods:** In this paper, material decomposition is solved in a variational framework. We investigate two fidelity terms: a weighted least squares (WLS) term adapted to Gaussian noise and the Kullback-Leibler (KL) distance adapted to Poisson Noise. The decomposition problem is then solved iteratively by a regularized Gauss-Newton algorithm.

**Results:** Simulations are performed in a numerical mouse phantom and the decompositions show that KL outperforms WLS for low numbers of photons, i.e., for high photonic noise.

**Conclusion:** In this work, a new method for material decomposition in SPCT is presented. It is based on the KL distance to take into account the noise statistics. This new method could be of particular interest when low-dose acquisitions are performed.

**Keywords:** Spectral CT, material decomposition, Kullback-Leibler distance, inverse problem

# 1. INTRODUCTION

Standard Computed Tomography (CT) is currently one of the most useful imaging modalities. However, X-ray detectors integrate the whole energy spectrum and since material attenuation depends on the incident X-ray photons energy, there is a loss of information. By performing energy-resolved acquisition, spectral CT is promising to further improve the interest of CT. In particular, 3D mono-energetic images or 3D density map of materials can be retrieved. SPCT has many potentials for clinical applications, in particular, to image high Z contrast agents [1, 2]. While, dual-energy CT have been proposed to acquire energy-resolved data [3], but with a large irradiation dose. Photon-counting detectors based on a new technology can count photons that are hitting the detector and retrieve their energy [4].

Three different approaches for decomposing

SPCT data into material volumes are available. A first approach consists in the tomographic reconstruction of the energy-resolved measurements followed by a material decomposition step in the object domain [5]. A second approach is to do the tomographic reconstruction and the decomposition at the same time, which is usually referred to as the one-step approach [6]. A third approach is to perform material decomposition in the projection domain followed by a tomographic reconstruction step [7]. The last approach is highly parallelizable (each projection angle can be processed individually) and it embeds all the non-linearities of the forward problem, i.e., tomographic reconstructions can then be performed linearly without any simplifying assumptions. In this work, material decomposition in the projection domain is considered.

Material decomposition is a non-linear ill-posed inverse problem that can be solved by minimizing a cost function composed of a data fidelity term
and a regularization term. A regularization parameter is chosen to balance the two terms of the cost function. SPCT data are corrupted by photonic noise, which is Poisson distributed. Roessl et al. and Schlomka et al. have proposed decomposition methods that take into account the statistics of the noise in the projection domain [8, 7], e.g., using negative log likelihood. However, the latter methods rely on zero-th or first-order minimization algorithms [9, 10] that have slow convergence rates. Recently, a fast method based on a Gauss-Newton algorithm was proposed [11]. Unfortunately, it uses a weighted least squares fidelity term that is not perfectly adapted to Poisson noise. In this work, the Kullback-Leibler distance is proposed to better handle the noise statistics and a Gauss-Newton algorithm is implemented to get fast decomposition.

This paper is organized as follows. In Section II, the forward model that maps the projected mass of the material onto the number of photons detected is presented. In Section III, we detail the regularization functional considered and the minimization algorithms used. Numerical experiments are presented in Section IV before the conclusion.

## 2. THEORY

### 2.1. Object model

Assuming that the object of interest can be decomposed on a basis of M materials [12], the linear attenuation coefficient is written

$$ \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$^{-1}$) at point $\mathbf{x}$ and energy E, $\rho_m(\mathbf{x})$ is the mass density (in g.cm$^{-3}$) of the m-th material at point $\mathbf{x}$, and $\tau_m(E)$ is the mass attenuation (in cm$^2$.g$^{-1}$) of the m-th material.

### 2.2. Forward model

The standard SPCT forward model [11, 13] is written as

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

where

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

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

### 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$, i.e., $s_{i,p} = s_i(\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$, i.e., $a_{m,p} = a_m(\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 projected mass density $\mathbf{a}$ from the data $\mathbf{s}$. Let $\mathbf{s} = \mathcal{F}(\mathbf{a})$ be the forward model defined in (2). The PMD is 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 Tikhonov regularization terms used in [11] is chosen

$$ \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.
### 3. PROPOSED MATERIAL DECOMPOSITION

#### 3.1. Fidelity term

The usual fidelity term is the WLS, defined by

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

where $||.||^2$ is the $l_2$-norm and $\mathbf{W}$ is the weighting matrix defined as $\mathbf{W} = \text{diag}(\frac{1}{\sqrt{s+1}})$ where $\text{diag}(\mathbf{x})$ represents a diagonal matrix with the elements $\mathbf{x}$. However the noise associated with photon detection problem is Poisson noise. The Kullback-Leibler divergence is the distance associated to Poisson noise [14]. The KL distance between two Poisson-distributed data $\mathbf{s}$ and $\mathcal{F}(\mathbf{a})$ is written [15]

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

where $\zeta$ is a small integer that avoids to dividing by zero and having the log of zero.

The corresponding data fidelity terms will be denoted $\mathcal{D}_{\text{KL}}$ and $\mathcal{D}_{\text{WLS}}$.

#### 3.2. Minimization of the cost function

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 (11) $$

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 (12) $$

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. They can be written as follow [16]

$$ \text{WLS}: \begin{cases} \mathbf{Z}_{\text{g}} = 2\mathbf{W}^{\text{T}}\mathbf{W}, \\ \mathbf{Z}_{\text{h}} = 2\mathbf{W}^{\text{T}}\mathbf{W} \end{cases}, \text{KL}: \begin{cases} \mathbf{Z}_{\text{g}} = \mathbf{Y}^{-1}, \\ \mathbf{Z}_{\text{h}} = \mathbf{Y}^{-2}\mathbf{S} \end{cases} $$

where $\mathbf{W}$ is the weighting matrix of the WLS term, $\mathbf{Y} = \text{diag}(\mathcal{F}(\mathbf{a}) + \zeta)$ and $\mathbf{S} = \text{diag}(\mathbf{s} + \zeta)$.

The step length $\lambda^k$ is computed as

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

### 4. NUMERICAL EXPERIMENTS AND IMAGE ANALYSIS

Material decomposition is performed in the numerical phantom DigiMouse [17]. Projections are generated with the radon Matlab function. The goal is to decompose $M = 2$ materials (soft tissues and bones) from projections of $P = 437 \times 992$ pixels.

We define the number of photons $N^0$ send at each pixel detector

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

Since $N^0$ controls the signal-to-noise ratio, a wide range is tested: $N^0 = [10^2, 10^{2.2}, ..., 10^4]$. The X-ray spectrum used is generated with a 90 kV source. Considering a detector with $I = 3$ energy bins (14–39 keV, 40–64 keV and 65–90 keV as energy bins).

The regularization parameter values tested are $\alpha = [10^{-2}, 10^{-1.8}, ..., 10^4]$. Our initial guess $\mathbf{a}^0$ was set as uniform material maps with $\mathbf{a}_{\text{soft}}^0 = 2$ g.cm$^{-2}$ and $\mathbf{a}_{\text{bone}}^0 = 0$ g.cm$^{-2}$. The algorithm is stopped when the relative decrease of the cost function is below 0.1%, if the computed step $\lambda^k$ is below $10^{-2}$ or if the maximum iteration number 50 is reached. The relative error is computed as

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

where $\mathbf{a}_{m}^{\text{dec}}$ is the decomposed map of the $m^{th}$ material and $\mathbf{a}_{m}^{\text{true}}$ is the ground truth for the material $m$. Moreover, $\hat{\alpha}$ denotes the regularization parameter that minimizes the relative error $\xi$

$$ \hat{\alpha} = \underset{\alpha}{\operatorname{argmin}} \xi(\alpha) \quad (16) $$

### 5. RESULTS AND DISCUSSION

In this section, we present results from experiments on DigiMouse.

In Figure 1 (left), the evolution of the fidelity terms $\mathcal{D}_{\text{KL}}$ and $\mathcal{D}_{\text{WLS}}$ is plotted as a function of the number of iterations. There is an apparent decrease of the regularized functionals, and local minima are obtained after four iterations. The algorithm stops since there is a low decrease of the cost function or a small step length $\lambda^k$, which represents a local minimum because $\delta \mathbf{a}^k$ is, by construction, a descent's direction.
<figure><img src="image_2.png" /><figcaption>Figure 1: Evolution of the fidelity terms for $\hat{\alpha}$ and for $N^0 = 10^{2.2}$ (left) and of the decomposition error with the number of photons (right)</figcaption></figure>

<table><caption>Table 1: Mean and standard deviation of time and the minimum and maximum of iterations to converge for $n^0 = 10^{2.2}$ with $\alpha = [10^{-2}, 10^{-1.8}, ..., 10^4]$</caption><thead><tr><th></th><th>Time (sec)</th><th>Iterations</th></tr></thead><tbody><tr><td>WLS</td><td>3.28 ± 1.57</td><td>[3-4]</td></tr><tr><td>KL</td><td>2.45 ± 0.09</td><td>[3-4]</td></tr></tbody></table>

Figure 1 (right) shows the evolution of the decomposition error for $\hat{\alpha}$ with increasing numbers of photons $N^0$ for the fidelity terms $D_{KL}$ and $D_{WLS}$. For real projections, the ground truth is not known and thus from this figure the improvement of the reconstruction results with the KL distance for the lower number of photons is highlighted. For these low numbers of photons the weighted least square is not a good approximation of the KL distance and this statistical divergence is more efficient to take into account the noise statistic. Similar curves are obtained when relative errors for a given material are considered.

Figure 2 shows the ground truth images of the PMD images and decompositions with WLS and KL for $N^0 = 10^{2.2}$. A low number of photons is chosen to highlight the performance of KL. The KL distance clearly retrieves more faithfully soft tissues and bones than WLS at a small number of photons.

In Table 1, the average time required in order to make one iteration (3.5 GHz Xeon E5-1650 v3 and 64 GiB of RAM) as well as the standard deviation for $\alpha = [10^{-2}, 10^{-1.8}, ..., 10^4]$ are shown for $N^0 = 10^{2.2}$. The minimum and the maximum number of iterations required for convergence are also displayed. Table 2 provides the value of the optimal $\alpha$, the computation time and the iterations number in order to converge for $N^0 = 10^{2.2}$. Both KL and WLS converge rapidly for a large range of $\alpha$'s and take almost the same CPU time. KL is found to outperform WLS for a low numbers of photons.

<table><caption>Table 2: Value, time and number of iteration require to converge for $\hat{\alpha}$ at $n^0 = 10^{2.2}$</caption><thead><tr><th></th><th>$\hat{\alpha}$</th><th>Time (sec)</th><th>Iteration number</th></tr></thead><tbody><tr><td>WLS</td><td>10<sup>2</sup></td><td>2.74</td><td>4</td></tr><tr><td>KL</td><td>10<sup>1.4</sup></td><td>2.47</td><td>4</td></tr></tbody></table>

## 6. CONCLUSION

In this work, methods for material decomposition are carried out in the projection domain in SPCT. The non-linearity of the forward problem is taken into account. Two data fidelity terms, one based on weighted least squares and the other on the Kullback-Leibler divergence, are compared. The regularized cost functions are minimized with a Gauss-Newton scheme. The methods are compared for 2D projections of a 3D realistic phantom for different numbers of photons. The Kullback-Leibler distance gives the best results for the lower number of photons. In the future, we plan to assess the method with other phantoms and 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.

Conflict of interest: none

## BIBLIOGRAPHY

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

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

[3] T. R. Johnson, "Dual-energy ct: general principles," *American Journal of Roentgenology*, vol. 199, no. 5_supplement, pp. S3-S8, 2012.

[4] K. Taguchi and J. S. Iwanczyk, "Vision 20/20: Single photon counting x-ray detectors in medical imaging," *Medical physics*, vol. 40, no. 10, 2013.

[5] C. Maaß, M. Baer, and M. Kachelrieb, "Image-based dual energy ct using optimized precorrection functions: A practical new approach of material decomposition
<figure><img src="image_3.png" /><figcaption>Figure 2: Ground truth of the material map of soft tissues and bones (left), decomposition with WLS (middle) and with KL (right) for $N^0 = 10^{2.2}$ photons for $\hat{\alpha}$</figcaption></figure>

in image domain," *Medical physics*, vol. 36, no. 8, pp. 3818-3829, 2009.

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

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

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

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

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

[11] N. Ducros, J. F. P. J. Abascal, B. Sixou, S. Rit, and F. Peyrin, "Regularization of nonlinear decomposition of spectral x-ray projection images," Nov. 2016, working paper or preprint. [Online]. Available: https://hal.archives-ouvertes.fr/hal-01391538

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

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

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

[15] S. Kullback and R. Leibler, "On information and sufficiency," *The annals of mathematical statistics*, vol. 22, no. 1, pp. 79-86, 1951.

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

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