2.1

Photon-counts and attenuations

Assuming a photon counting detector with B bins and according to Equation 2, the photon counts may be corrupted with Poisson noise. The measures become

We only use m_b(θ, p) in the rest of the paper and explicitely indicate wether data are corrupted with noise. The projections are then log-transformed according to

where m_0,b is the number of photons without object.

2.2

The polynomial model

We look for a a simplified model of the inverse mapping (A₁, …, A_m) = Φ⁻¹(s₁, …, s_B). We choose a polynomial model. Each material sinogram A_m is approximated by a polynomial ψ_m,D of degree D in the variables (s₁, …, s_B). Formally,

where k = (k₁, …, k_B) is a multi-index, |k| = k₁ + ⋯ + k_B and . We define N(D, B) the total number of coefficients of a polynomial of degree D in B variables, e.g., N(2,2)=6 and N(3,2)=10. Note also that the coefficients must be determined for each basis material m. If D, B and M are fixed, M × N(D, B) coefficients have to be determined. For example, if B = 2, M = 2 and D = 3, we seek 2 × N(3, 2) = 20 coefficients. Note that the same polynomial is applied to all the pixels of the sinogram A_m, i.e., that the source spectrum and the detector response are uniform over the beam and the detector, respectively.

2.3

The consistency metric

In 2D parallel geometry, the sought A_m are the Radon transform of the material map a_m

where and are perpendicular.

To account for the spectral nature of the decomposition problem, we combine the material sinograms A_m into mono-energetic sinograms C_n, in view of applying the consistency metric to them. We choose N energy levels E_n in the energy range of the source and form the mono-energetic sinograms C_n

The coefficients f_m(E_n) are known (see Equation 1).

A consistency condition of the Radon transform states that the integral of each projection (the order-0 moment) does not depend on the projection angle θ since each of these integrals equals the integral of the attenuation coefficient over the object. We define the moment J_n(θ) of the C_n by

We use the variance of J_n to evaluate if J_n is constant over Θ. The consistency function hence reads:

where denotes the mean of J_n over all θ. Note that the function ℓ_n depends on the polynomial coefficients since all A_m do (see Equation 5). Finally, we define M consistency functions on the material sinograms A_m in a similar way. The final consistency metric accumulates the consistency of all computed mono-energetic sinograms and all material-specific sinograms

The consistency metric would evaluate to zero on perfectly consistent material sinograms.

2.4

Minimization

Due to the hardening of the beam in each bin (see e.g. Ref. 8), the measured attenuations s_b do not satisfy the consistency condition. The loss ℓ is minimized with respect to the coefficients c to achieve mono-energetic and material sinograms which are as consistent as possible.

Since we use only order-0 consistency conditions, a constant sinogram (i.e. A_m(θ, p) = constant for all θ and p) is perfectly consistent. To prevent the minimization to output such undesirable solution, we follow the idea of Ref. 8 and constrain the minimization by some known values. First, if there is no attenuation in all bins (s_b = 0, ∀b ∈ {1, …, B}), we enforce 0 in all material and mono-energetic sinograms by setting for all m. Second, there still is a trivial solution to the minimization of the consistency loss function: the null sinogram. We enforce, for each material m, a particular value in one voxel of each reconstructed material map a_m. To this end, a small sample of each material is placed in the field-of-view and a reconstruction from raw data is computed. The small samples are easily identifiable in the reconstruction. Let be one voxel in each material sample and assume reconstructions are computed with a standard Filtered Backprojection (FBP) algorithm. Since FBP is a linear operation, one has

The values FBP(s^k) are easily computed once, before the minimization. Each material map is constrained by exactly M relations, which take the form

3.1

Simulation of data

Numerical experiments used a 2D phantom made of an outer water disc of diameter 32 mm and five bone inserts (with diameters ranging from 2 to 5 mm), placed inside the water disc (see Figure 2). Two tiny inserts of bone and water (1 mm in diameter each) were placed outside the phantom. One voxel in each insert was chosen for the constraints. The material sinograms A_m were analytically computed with RTK (Ref. 9). The simulated sinograms had 700 pixels with 0.05 mm spacing and 720 projections over a 180° angular range. Two effective spectra were used. The low-energy (LE) and high-energy (HE) spectra had a tube-voltage of 80 keV and 120 keV respectively (Figure 1). Without object, the detector received a total number of photons of 1.3 × 10⁶ and 2.9 × 10⁶ photons for the LE and HE spectra respectively. The photon counts m_b were computed by applying Equation 2, then log-transformed according to Equation 4.

Figure 1.

The low-energy (LE) and high-energy (HE) spectra.

Figure 2.

Top: Material map obtained with the DCC-based decomposition. Grayscale is : 1 ± 0.3. The red line indicates the profile used in Figure 4. Botoom: Profiles along the red line for the DCC-based and the calibrated water (left) and bone (right) maps. The consistency metric at convergence was ℓ(c) = 0.0208. If evaluated on the calibrated sinograms, the consistency function was 0.0018.

The degree of the sought polynomials was fixed to D = 3 and the consistency metric Equation 10 was minimized under the constraints defined above, with the Sequential Least Squares Programming algorithm. The total number of estimated polynomial coefficients was 18. The initial guess was always set to zero for all coefficients.

3.2

Evaluation methods

Our method was compared to the calibration from a set of dedicated measurements with the same LE and HE spectra over a set of water and bone lengths. The set covered all combinations which were present in the phantom. More precisely, 100 equi-spaced lengths of each material were measured. Water lengths ranged from 0 to 32.6 mm and bone lengths ranged from 0 to 10.97 mm. All these combinations were irradiated with the same LE and HI spectra as the phantom. Then the polynomials were fitted to the calibration data and further applied to the phantom data to produce a reference polynomial decomposition.

We compared the poly-energetic reconstructions from raw-data with mono-energetic images computed from our DCC-based material maps and from the reference calibrated material maps.

4.1

Noise-less data

Results of the decomposition are shown in Figure 2. Water and bone are adequately separated. The profiles in Figure 2 (bottom row) indicate residual cross-talk between the two material maps. In the center of the water phantom, the low-constrast feature is visible.

Since the consistency is enforced on the mono-energetic images, Figure 3 compares poly-energetic images, DCC-based and calibration-based mono-energetic images. Poly-energetic images clearly suffer from severe beam-hardening, which is almost completely corrected on both mono-energetic images. The profiles presented in Figure 4 reveal that the DCC-based and calibrated 40 keV images can hardly be distinguished. A slight discrepancy between the 80 keV images still subsists though, especially in the vicinity of the border of the phantom.

Figure 3.

Poly-energetic reconstructions (left) from LE (top) and HE (bottom) data. DCC-based (middle) and calibration-based (right) mono-energetic images at 40 keV (top) and 80 keV (bottom). Grayscale for all images : 0 ± 250HU.

Figure 4.

Profile poly-energetic reconstructions from raw data and from mono-energetic images computed from DCC-based and calibrated material maps.

4.2

Robustness to noise

The photons count measurements m_b were corrupted with Poisson noise according to Equation 3. The reference photon flux is given by the spectra in Figure 1, i.e. 1.3 × 10⁶ and 2.9 × 10⁶ photons for the LE and HE spectra respectively. The noise level was set by reducing the total number of emitted photons by a factor 1, 10 and 100. The influence of noise on the material maps is presented in Figure 5. The quality of the images is significantly degraded. The consistency function value at convergence increases with the level of noise. It is 0.2108 for reference noise level (factor 1), 1.4738 at factor 10 and 12.815 at factor 100. We expect the choice of the reference voxel to play a critical role in the presence of heavy noise.

Figure 5.

Mono-energetic images at 40 keV (top) and 80 keV (bottom), for Poisson noise with photon counts donwscaled by a factor 1, 10 and 100 (from left to right). Grayscale for all images : 0±250HU. The consistency function at convergence was ℓ(c) = 0.2108, 1.4738, 12.815 for factor 1, 10 and 100 respectively.