Odran Pivot
email: odran.pivot@creatis.insa-lyon.fr

Clarisse Fournier

Joachim Tabary

Jean Michel Létang

Simon Rit

J M L Étang

Scatter correction for spectral CT using a primary modulator mask

Keywords: Photon counting detectors, Primary modulator mask, Scatter correction, Spectral CT

come

Scatter correction for spectral CT using a primary modulator mask

I. INTRODUCTION

Photon counting detectors (PCD) classify photons according to their energy. Their emergence has triggered the development of energy-resolved computed tomography (or spectral CT). This new modality gives the possibility to retrieve information about the chemical composition of the inspected object [START_REF] Taguchi | Vision 20/20: Single photon counting x-ray detectors in medical imaging[END_REF]. When considering flat-panel detectors, the scattered radiation, which induces a bias on the detected images and leads to a loss of contrast and artifacts in CT, cannot be neglected [START_REF] Siewerdsen | Cone-beam computed tomography with a flat-panel imager: Magnitude and effects of x-ray scatter[END_REF]. In order to fully benefit from the spectral information, scattercorrected images are required.

Since the early developments of x-ray imaging with integration-mode detectors, the problem of scattered radiation correction gave rise to various methods [START_REF] Rührnschopf | A general framework and review of scatter correction methods in x-ray cone-beam computerized tomography. Part 1: Scatter compensation approaches[END_REF]. These methods can be classified into two groups: hardware-based and softwarebased methods.

The group of hardware-based methods includes approaches based on geometry modifications or the use of additional devices to reduce the part of scattered photons reaching the detector. The most common approach consists in collimating the pixels of the detector by using anti-scatter grids, the disadvantage of which is to increase the dose while not rejecting totally the scatter [START_REF] Siewerdsen | The influence of antiscatter grids on soft-tissue detectability in cone-beam computed tomography with flat-panel detectors[END_REF]. The amount of scattered radiation reaching the detector can also be lowered by increasing the air-gap [START_REF] Sorenson | Scatter rejection by air gaps: An empirical model[END_REF] [START_REF] Persliden | Scatter rejection by air gaps in diagnostic radiology. Calculations using a Monte Carlo collision density method and consideration of molecular interference in coherent scattering[END_REF] or by limiting the field-of-view with additional collimation [START_REF] Chen | Feasibility of volume-of-interest (VOI) scanning technique in cone beam breast CT-a preliminary study[END_REF] [START_REF] Lai | Reduction in x-ray scatter and radiation dose for volume-of-interest (VOI) cone beam breast CT-a phantom study[END_REF].

The second group, software-based methods, covers techniques aiming at estimating the scattered radiation map and using it to correct the measured images [START_REF] Rührnschopf | A general framework and review of scatter correction methods in cone beam CT. Part 2: Scatter estimation approaches[END_REF]. This group can again be divided into two sub-families: model-based and measurement-based methods. Model-based methods estimate the scatter map from the measured data, by modelling the physics of the acquisition. The models used in the literature are based on scatter-kernels [START_REF] Rinkel | A new method for xray scatter correction: first assessment on a cone-beam CT experimental setup[END_REF] [START_REF] Maltz | Algorithm for X-ray Scatter, Beam-Hardening, and Beam Profile Correction in Diagnostic (Kilovoltage) and Treatment (Megavoltage) Cone Beam CT[END_REF] [START_REF] Sun | Improved scatter correction using adaptive scatter kernel superposition[END_REF] or more complex Monte-Carlo approaches [START_REF] Jarry | Scatter correction for kilovoltage cone-beam computed tomography (CBCT) images using Monte Carlo simulations[END_REF] [14] [START_REF] Bootsma | Efficient scatter distribution estimation and correction in CBCT using concurrent Monte Carlo fitting[END_REF]. Measurement-based methods seek to estimate the scatter distribution with the aid of a physical beam modifier. Methods using static or moving fully-attenuating beam-blockers directly estimate scatter in beam-blocked regions of the detector [START_REF] Ning | X-ray scatter correction algorithm for cone beam CT imaging[END_REF] [17] [START_REF] Siewerdsen | A simple, direct method for x-ray scatter estimation and correction in digital radiography and cone-beam CT[END_REF] [START_REF] Ouyang | A moving blocker system for cone-beam computed tomography scatter correction[END_REF]. Several years ago, an approach consisting in using a semitransparent modulator mask with a perfectly uniform highfrequency pattern between the source and the object has been introduced [START_REF] Maltz | Cone beam X-ray scatter removal via image frequency modulation and filtering[END_REF] [START_REF] Zhu | Scatter Correction Method for X-Ray CT Using Primary Modulation: Theory and Preliminary Results[END_REF]. Assuming that the primary map is modulated in the high-frequencies and the scatter mostly has low-frequencies, the corrected primary map is recovered after low-pass filtering in the Fourier domain and demodulation. This method is hard to implement in practice, due to the difficulty to manufacture such a primary modulator. More recently, a promising alternative has been proposed which allows the use of a mask with a non uniform pattern. Instead of removing the scatter components by filtration in the Fourier domain, the scatter map is estimated by solving a patch-based local optimization problem in the projection domain [START_REF] Yan | Iterative scatter correction for x-ray cone-beam CT with semi-transparent beam stop array[END_REF] [23] [START_REF] Bier | Scatter correction using a primary modulator on a clinical angiography C-arm CT system[END_REF], under the assumptions that the scatter map is smooth and the primary map is locally smooth [START_REF] Chen | Optimization-based scatter estimation using primary modulation for computed tomography[END_REF]. Alternatively, a global formulation of the optimization problem has been proposed, as well as an edge-preserving weighting which increases the robustness when the primary map contains some discontinuities [START_REF] Chen | Optimization-based scatter estimation using primary modulation for computed tomography[END_REF].

The previously described scatter correction methods were designed for integration mode detectors and do not exploit the spectral information provided by energy-resolved PCD. A modified version of the beam-stop method, adapted to multi-energy imaging, has been proposed [START_REF] Sossin | A novel scatter separation method for multi-energy x-ray imaging[END_REF] and shows a very good accuracy in practice [START_REF] Sossin | Experimental validation of a multi-energy x-ray adapted scatter separation method[END_REF]. However, it needs an extra acquisition which increases the dose and the exposure.

This study aims at presenting a method inspired by [START_REF] Ritschl | Robust primary modulation-based scatter estimation for cone-beam CT[END_REF] and [START_REF] Chen | Optimization-based scatter estimation using primary modulation for computed tomography[END_REF] and adapted to spectral CT. Our main contributions are a correction matrix which compensates for the effect of the modulator mask taking into account the spectral information, a scatter model based on B-spline basis functions and a novel cost function.

II. MATERIALS AND METHODS

The first section (II-A.1) presents the basic theory of scatter estimation using a primary modulator mask as introduced in [START_REF] Ritschl | Robust primary modulation-based scatter estimation for cone-beam CT[END_REF], but adapted to spectral CT imaging. We then detail our contributions: a model-based correction matrix which uses the spectral information to compensate for the mask attenuation (II-A.2), a B-spline model for the representation of the scatter map (II-A.3), specific weighting matrices aiming at improving the method accuracy (II-A.4), the complete cost function (II-A.5) and the optimization algorithm (II-A.6).

A. Scatter correction method 1) Principle: Throughout this paper, the energy dimension is discretized and we note

a u,v,θ = l∈L u,v,θ µ(λ θ + lζ u,v,θ ) dl (1)
the vector (length: N E ) of integrals of the energy-dependent linear attenuation coefficient (LAC) map µ of the object, along line L u,v,θ going through the 3D source position λ θ in direction ζ u,v,θ of pixel with coordinates (u, v). Each element of the vector a u,v,θ corresponds to a sample of the energy dimension. The total x-ray spectrum N t u,v,θ (number of photons, length: N E ) reaching the detector at pixel (u, v) and projection angle θ can be decomposed as the sum of the scatter spectrum N s u,v,θ and the primary one N p u,v,θ . The latter follows the Beer-Lambert law:

N p u,v,θ = N 0 u,v exp(-a u,v,θ ), (2)
N 0 u,v being the x-ray spectrum incident on the detector at pixel position (u, v) without object.

We consider a spectral CT scanner with a flat panel photon counting detector of N B energy bins (N B is typically smaller than N E ). The total transmission spectrum t u,v,θ (dimensionless, length: N B ) detected in pixel (u, v) and at projection angle θ is computed as the raw data normalized by the ones without object (both in detected number of photons, length: N B ):

t u,v,θ = (D u,v N t u,v,θ ) (D u,v N 0 u,v ) (3)
= diag D u,v N 0 u,v ) -1 D u,v N t u,v,θ , (4)
where D u,v ∈ R N B ×N E denotes the pixel-dependent detector response matrix (DRM) which columns give the probability distribution for a photon at a given energy to be detected in an energy-bin and the Hadamard division operator. It can be further decomposed as the sum of the primary and the scatter components (p u,v,θ and s u,v,θ respectively, each of length N B ):

t u,v,θ = p u,v,θ + s u,v,θ , (5)
where the primary transmission spectrum p u,v,θ is given by the following expression:

p u,v,θ = diag D u,v N 0 u,v ) -1 D u,v N p u,v,θ . (6)
When a primary modulator mask made with an homogeneous material m is placed between the source and the object (as shown in Figure 1), the primary map is modulated and becomes:

pu,v,θ = diag D u,v N 0 u,v ) -1 D u,v N p u,v,θ , (7)
with

N p u,v,θ = N 0 u,v exp(-a u,v,θ -µ m l m u,v ) (8)
= diag exp(-µ m l m u,v ) N p u,v,θ (9)
where l m u,v is the length of modulator mask crossed by the lines L u,v,θ and µ m the LAC of the material constituting the primary modulator. Note that some rays do not intersect the mask since it has holes and l m u,v = 0 for those rays. An example of such primary modulator mask is shown in Figure 2. Furthermore, the mask modifies the x-ray beam and we note s u,v,θ the scatter transmission spectrum obtained with a beam modulated by the mask. The measured transmission map tu,v,θ is then: The objective of the proposed method is to retrieve the demodulated primary, and it is necessary to provide a correction scheme allowing to recover it given the measurement tu,v,θ and the scatter s u,v,θ . If the detector response matrices D u,v were invertible, one could recover p u,v,θ from Equations 9 and 10 with the knowledge of the scatter spectrum:

tu,v = pu,v + s u,v . (10
p u,v,θ = C u,v ( tu,v,θ -s u,v,θ ), (11)
where C u,v ∈ R N B ×N B is the correction matrix associated with the pixel (u, v), which expression is:

C u,v = diag D u,v N 0 u,v ) -1 D u,v diag(exp(µ m l m u,v )) D -1 u,v diag D u,v N 0 u,v ). (12
) This matrix therefore transforms a measured spectrum with mask to the spectrum which would have been measured without mask. One can verify that the correction matrix is the identity matrix when l m u,v = 0. In practice, D u,v is not square and therefore not invertible, and we discuss its practical computation in section II-A.2.

The measured, primary and scatter spectral sinograms can be flattened into vectors denoted t, p and s, respectively. For example, t is defined as follows: t = [ tT u0,v0,θ0 , tT u1,v0,θ0 , . . . , tT

u N U ,v N V ,θ N Θ ] T , (13)
with T the transpose operator and N U , N V and N Θ the numbers of pixels in the two directions of the detector and projection angles acquired, respectively. Similarly, the correc-

tion matrix C ∈ R N B N U N V NΘ×N B N U N V NΘ is constructed by concatenating the matrices C u,v such that p = C( t -s). (14)
As in [START_REF] Ritschl | Robust primary modulation-based scatter estimation for cone-beam CT[END_REF], we now consider an estimation of the scatter map s with an error ε given by ε = ss. By developing Equation 14, we obtain the estimated primary map p:

p = C( t -s) = C( t -s + ε) = p + Cε. ( 15
)
One can observe that an error ε on the estimated scatter map induces an error Cε on the estimated primary map. Assuming that both the scatter map and its estimate are smooth (this point is explained in section II-A.3), the error ε is also smooth and therefore the error Cε on the estimated primary map contains the same discontinuities as the correction matrix C.

In other words, the modulation pattern of the mask remains visible in the estimated primary map as long as ε is not null. Considering a primary modulator with sharp edges, one can say that the presence of the modulation pattern increases the spatial gradient of the estimated primary. A direct adaptation of [START_REF] Ritschl | Robust primary modulation-based scatter estimation for cone-beam CT[END_REF] to spectral CT is to estimate the smooth scatter map by minimizing the weighted gradient of the scatter-corrected primary:

s ∈ arg min s∈R N B N U N V N Θ s smooth W ∇C( t -s) 1 , (16)
where

W ∈ R 2N B N U N V NΘ×2N B N U N V NΘ is a diagonal weighting matrix and ∇ ∈ R 2N B N U N V NΘ×N B N U N V NΘ
is the numerical spatial gradient of the image along the two dimensions u and v of the detector.

After detailing the correction matrix in section II-A.2, several additions to this simple cost function are proposed in the following sections. First, in section II-A.3, a B-spline model of the scatter map is proposed to reduce the dimensionality of the scatter estimation problem. Then, in section II-A.4, specific weighting matrices W are introduced to improve the validity of Equation 16. The complete cost function (Equation 25) with a regularization term and specific weights is presented in section II-A.5.

2) Correction matrix: The correction matrix C u,v associated with the pixel (u, v) cannot be computed from Equation 12because the DRM is not invertible and generally not known. We propose instead a model-based estimation using an empirical model for the rows of the correction matrix C u,v defined as:

c b,u,v (b ) = h(b -b) c 0 b,u,v δ(b -b)-c 1 b,u,v exp(-c 2 b,u,v (b -b)) , ( 17
) with c b,u,v ∈ R N B the row of C u,v corresponding to energy bin b, [c 0 b,u,v , c 1 b,u,v , c 2 b,u,v
] a triplet of scalar parameters and h and δ the Heaviside and Dirac functions.

The parameters of the proposed empirical model for each pixel and energy bin are then estimated through an offline calibration procedure using a scatter-free dataset of N cal pairs of spectral transmission images with various attenuation levels, both with and without primary modulator mask ( pk and p k respectively):

[c 0 b,u,v , c1 b,u,v , c2 b,u,v ] ∈ arg min c 0 b,u,v ∈R * + c 1 b,u,v ∈R * + c 2 b,u,v ∈R * + N cal k=1 (p k,b,u,v -c T b,u,v pk,u,v ) 2 .
(18) In this work, the minimization is performed independently for each pixel with Nelder-Mead's downhill simplex algorithm. An example of such a correction matrix with N B = 8, for a pixel corresponding to an x-ray attenuated by the mask, is presented in Figure 3. The estimation of the correction matrix is object independent and therefore performed only once in an offline calibration step, before the acquisition of the object images.

3) Scatter model: It is well known that high frequencies of x-ray scatter maps are negligible in the two spatial directions of the flat-panel detector and in the angular direction of the projections [START_REF] Bootsma | Spatial frequency spectrum of the x-ray scatter distribution in CBCT projections[END_REF]. Aiming at an accurate representation of the scatter map with a minimal number of parameters, we chose to represent it using B-spline basis functions, which are commonly used for modelling smooth maps [START_REF] Unser | B-spline signal processing. I. Theory[END_REF]. In this paper, we used three dimensional quadratic B-spline basis functions (two dimensions along the flat detector and one along the projection angles), and we recover the scatter map knowing the N K parameters x e,k for each energy bin with the following expression:

s b,u,v,θ = N K k=1 β k (u, v, θ)x b,k , (19)
with β k (u, v, θ) the k-th B-spline basis function evaluated at pixel (u, v) and projection θ. The definition of the B-spline basis function is

β k (u, v, θ) = β 2 u -u k ∆ u β 2 v -v k ∆ v β 2 θ -θ k ∆ θ , (20)
where (u k , v k , θ k ) and (∆ u , ∆ v , ∆ θ ) are the B-splines knots positions and spacings respectively. β 2 (z) is the quadratic Bspline basis function which expression is given by:

β 2 (z) =          1 2 (z 2 + 3z + 9 4 ), ∀z ∈ [-3 2 , -1 2 [, 3 4 -z 2 , ∀z ∈ [-1 2 , 1 2 [, 1 2 (z 2 -3z + 9 4 ), ∀z ∈ [ 1 2 , 3 2 [, 0, ∀z ∈] -∞, -3 2 [∪[ 3 2 , +∞[. (21)
The number of parameters N K needed to represent the scatter map in one energy bin corresponds to the number of B-spline knots chosen, and depends on the spacing between two knots in each of the three dimensions. The choice of the optimal steps is discussed in section IV.

We denote B ∈ R N B N U N V NΘ×N B N K the matrix containing the B-spline basis functions and x the vector of the model parameters x b,k . The scatter map is then obtained from x with

s = Bx. ( 22
)
4) Weights: As explained in II-A.1, the scatter can be estimated by minimizing the primary gradient. However, it is not desirable to estimate a scatter which would compensate for the spatial gradient of the primary map. The use of novel specific weighting matrices in the cost function tackles this issue.

The primary gradient evaluated in pixels out of the modulator mask edges (either in the holes or in regions with constant thickness) contains only information about the object structures, and nothing about the persistence of the mask footprint in the primary map. As a consequence, we only need to evaluate the spatial gradient of the primary map on the mask edges. This can be done by weighting the spatial gradient of the primary map by the spatial gradient of the mask image, which can be achieved using the following weighting matrix:

W 1 = diag(∇ p0 ). (23)
However, object structures may also overlap with the modulator mask edges and the cost function should account less for the spatial gradient at these locations. We reused the weighting introduced in [START_REF] Chen | Optimization-based scatter estimation using primary modulation for computed tomography[END_REF] to solve the same issue:

W 2 = diag (exp(-η∇ p)) , (24)
with η a scalar parameter and ∇ p ∈ R 2N B N U N V NΘ an estimation of the primary gradient map from the measured sinogram m. In our work, ∇ p was directly estimated from ∇C t in pixels where there is no mask structure (i.e. where ∇ p0 = 0), thus neglecting the spatial gradient of the scatter map, and it was linearly interpolated between those pixels. 5) Cost function: Following the B-spline model detailed in section II-A.3, the scatter estimation comes down to the estimation of the B-spline coefficients at each energy bin. We stabilize the estimation by adding a Tikhonov regularization, on the parameters directly as spatial smoothness is inherently constrained by the B-spline model. We obtain the following estimation of the scatter map parameters:

x ∈ arg min x∈R N B N K W 1 W 2 ∇C( t -Bx) 1,σ + T x 2 2 , (25)
where

• W 1 and W 2 are the weighting matrices described above and • 1,σ is a smooth approximation of the 1 norm for this data-fidelity term [START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF] defined as

y 1,σ = n y 2 n + σ 2 -σ, (26)
with σ a strictly positive smoothing scalar parameter; • T is the diagonal matrix which energy-dependent elements τ b control the strength of Tikhonov's regularization in each energy bin. From the estimate of the scatter map parameters x, the primary map is deduced by incorporating Equation 22into Equation 14:

p = C( t -B x). (27)
6) Optimization: We propose to solve Equation 25 using Newton's method, starting with an initial guess x (0) = 0 and building new estimates with the following update rule:

x (n+1) = x (n) + δ (n) x , (28)
where δ

(n)

x ∈ R N B N K denotes Newton's step at iteration (n). Newton's step is obtained solving

H[F](x (n) ) + T T T δ (n) x = -∇[F](x (n) ) -T T T x (n) , (29)
with ∇[F](x (n) ) and H[F](x (n) ) the gradient and the Hessian of the data-fidelity term at iteration (n), respectively. Their expressions are given by

∇[F](x (n) ) = A T M (n) (d -Ax (n) ) H[F](x (n) ) = A T (M (n) -N (n) )A ( 30
)
where

           d = W 1 W 2 ∇C t A = W 1 W 2 ∇CB M (n) = diag (d -Ax (n) ) 2 + σ 2 -1 2 N (n) = diag (d -Ax (n) ) 2 (d -Ax (n) ) 2 + σ 2 -3 2 . ( 31
)
In practice, one iteration of Equation 28 was sufficient to recover the scatter map in our experiments which indicates that the cost function is nearly quadratic.

B. Reconstruction

For CT reconstruction, the transmission sinograms (t, p and p) are converted to attenuation sinograms (g t , g p and g p respectively) by computing their negative natural logarithm. For example, the attenuation sinogram of the estimated primary map is given by

g p = -ln( p). (32)
The attenuation sinograms at each energy bin are then reconstructed with the Reconstruction Toolkit (RTK) [START_REF] Rit | The Reconstruction Toolkit (RTK), an open-source cone-beam CT reconstruction toolkit based on the Insight Toolkit (ITK)[END_REF] using Feldkamp's filtered backprojection algorithm [START_REF] Feldkamp | Practical cone-beam algorithm[END_REF]. For the parallel fan-beam geometry (section II-C.2), each slice is reconstructed separately using the same implementation.

The reconstructed volumes of total, reference and estimated primary are denoted f t , f p and f p, respectively.

C. Test cases

Both simulations and acquisitions were realized to evaluate the proposed method. The corresponding setups are described in II-C.1 and II-C.2, respectively.

1) Simulations: The medical CT scanner simulated, which geometry is described in Figure 1 and Table I, had a tungsten x-ray tube with an anode angle of 30 • operated at 120 kV and 0.2 mAs. A 1024 × 1024 pixels flat-panel energy-resolved photon counting detector with a pitch of 400 µm and a 3 mm thick CdTe sensor has been simulated using an internal software, imitating an existing energy-resolved detector [START_REF] Brambilla | CdTe linear pixel X-ray detector with enhanced spectrometric performance for high flux X-ray imaging[END_REF]. The same simulated DRM was used for all pixels, which has been binned in the energy domain in 8 energy bins as a compromise between the number of bins and the noise level within each bin. The DRM is shown in Figure 4. The pixels were also binned such that one projection pixel regroups 16 detector pixels (4 × 4). Pile-up effects were not acquired offline without object and large doses can therefore be used. We simulated a CT scan of head of the ICRP realistic anthropomorphic phantom [START_REF] Menzel | ICRP Publication 110. Realistic reference phantoms: an ICRP/ICRU joint effort. A report of adult reference computational phantoms[END_REF] with 360 deterministic primary projections and 36 scatter projections (to limit the computation burden) using fixed-forced detection [START_REF] Poludniowski | An efficient Monte Carlo-based algorithm for scatter correction in keV cone-beam CT[END_REF] with one million particles per projection, both with Gate [START_REF] Jan | GATE V6: a major enhancement of the GATE simulation platform enabling modelling of CT and radiotherapy[END_REF]. The scatter maps were then interpolated with quadratic B-splines in order to obtain the 360 projections corresponding to the simulated primary maps. The primary modulator mask consists of a 5 mm thick graphite slab, drilled by 3 mm wide square holes every 5 mm.

The various attenuation levels taken for the calibration dataset (used to estimate the correction matrix) correspond to the combination of nine 2 cm thick polymethyl methacrylate slabs and a 1 mm aluminum slab.

Second order B-splines, with 128 pixels (≈ 20 cm) between two knots in the two directions of the detector and 15 • in the angular direction were set. The choice of the spacings between B-spline knots and the scalar parameters is discussed later in section IV. The scalar parameters η and σ were also tuned manually, as well as the regularization parameters τ b . Their values are given in table II. 2) Physical experiments: The physical experiments were performed using a parallel fan-beam set-up. The system is made of a tungsten x-ray source, powered with 120 kV. The detector is a line-detector without anti-scatter grid composed of 5 commercial detectors (Detection Technology ME100) for a total of 640 pixels with a pitch of 0.8 mm. 16 fan-beam acquisitions were stacked such that each projection had 640 × 16 isotropic pixels. The ME100 has 64 energy-bins of about 2 keV width, which are binned into 8 energy-bins described in Table III. The primary modulator mask is a 5 mm thick graphite slab drilled with 3 mm diameter holes each 5 mm (see Figure 2). The design of the mask is discussed later in section IV. The phantom is placed on a vertical translation table (allowing us to acquire various slices). The support of the mask is attached to the translation table. The system contains three collimators: the detector case and two additional collimators placed at the level of the source and just before the object. Unlike the detector case, the additional collimators are removable. The two imaged phantoms are a 33 × 27 cm 2 water-equivalent plastic ellipsoid (mimicking a pelvis) with tissue-equivalent plastic inserts of 3 cm diameter (CIRS phantom, model 062MQA, stoichiometry available in [START_REF] Hudobivnik | Comparison of proton therapy treatment planning for head tumors with a pencil beam algorithm on dual and single energy CT images[END_REF]) and a realistic human head phantom (RANDO phantom). Pictures of the two phantoms are shown in Figure 5. The calibration dataset was acquired with the combination of three 5 cm thick polypropylene slabs and two 5 mm thick aluminum slabs. The image without object, the mask image p0 and the calibration images were acquired with 20 mAs while the 180 projections of the tomographic acquisitions of the phantoms were acquired with 2 mAs. For each phantom, three CT scans were acquired: two without collimation (with and without the modulator mask) and one without mask and with fan-beam source collimation to irradiate only the one-line detector. The latter will be considered as a scatter-free reference acquisition for the rest of this paper following, e.g., [START_REF] Ritschl | Robust primary modulation-based scatter estimation for cone-beam CT[END_REF], [START_REF] Chen | Optimization-based scatter estimation using primary modulation for computed tomography[END_REF]. The scalar parameters used in the method are listed in Table II. The spacings between two knots on the detector were set to 64 pixels for the CIRS case (≈ 5 cm) and 128 pixels in the RANDO case (≈ 10 cm). In both cases, the spacing between two knots in the angular dimension of the projections was set to 10 • and second order B-splines were used.

1 × 10 1 1 × 10 1 σ 3 × 10 -3 1 × 10 -4 1 × 10 -4 τ 1 5 × 10 0 1 × 10 -6 1 × 10 -6 τ 2 2 × 10 -1 1 × 10 -5 6 × 10 -7 τ 3 1 × 10 -1 4 × 10 -5 4

D. Evaluation

The quantitative evaluation in the projection domain was performed on simulations with the mean relative error (M RE) on the transmission sinograms. For example, the mean relative error associated with the estimated primary at energy bin b is computed as follows:

M RE p b = 100 × 1 N U N V N Θ N U u=1 N V v=1 NΘ θ=1 |p b,u,v,θ -pb,u,v,θ | p b,u,v,θ , (33)
where p is the noise-free primary transmission sinogram. The errors M RE p b and M RE t b associated with the reference primary and the total sinograms are computed by replacing p by p and t in Equation 33, respectively. The transmission profiles of total, simulated primary, estimated primary, simulated scatter and estimated scatter at energy bins 2 and 6 are plotted on Figure 7. The profiles location is shown on the top-left sub-figure of Figure 6. Table IV shows the mean relative error on simulated data for each one of the 8 energy bins. The first column shows the error on the simulated total map, the second one on the simulated primary map and the third one on the estimated primary map. The results in the object domain, for the simulation case, are presented in Figure 8. The top-left sub-figure shows a CT slice and is divided into two parts: the simulated total map on the left and the estimated primary map on the right. The top-right sub-figure shows the profiles at energy bin 2, which location is shown on the top-left sub-figure. The bottom sub-figures show the mean spectra in voxels of two tissues of interest (TOI) corresponding to human tissues: the brain (TOI1) and the cranium spongiosa (TOI2). The blue, green and red curves represent the simulated total map, simulated primary map and estimated primary map, respectively.

III. RESULTS

A. Simulations

B. Physical experiments

The results on the CIRS phantom are shown in Figure 9, where the first, second and third rows correspond to the total map, reference primary map and estimated primary map, respectively. The first column shows attenuation sinograms and the second the corresponding CT slices. Each column is divided in two, where the left part is the second energy bin and the right part the sixth one.

Figure 10 shows the mean spectra of three regions of interest (ROI) of the CIRS phantom which locations are shown Mean spectra of TOI2 f t f p f p Fig. 8. Top left: CT slice at bin 2 of the simulated total map on the left and the estimated primary map on the right; top right: profiles at bin 2; bottom left: mean spectra of TOI 1 (brain); bottom right: mean spectral of TOI 2 (cranium spongiosa); the color code is the same for all graphs, i.e., simulated total map in blue, reference primary map in orange, estimated primary map in green.

on the top-left figure. The three ROIs are tissue-like inserts representing lung (ROI1), soft tissue (ROI2) and bone (ROI3). Blue, orange and green curves are the total, reference primary map and estimated primary map respectively.

Finally, the results on the RANDO anthropomorphic phantom are presented in Figure 11, where the top-right sub-figure shows a CT slice of the total (left) and estimated primary (right) maps. The top-right sub-figure shows the mean spectra of a cranium region of interest which location is shown on the top-left sub-figure. The two bottom sub-figures show the profiles (see location on top-left sub-figure) at bin 2 (bottomleft) and bin 6 (bottom right). The blue curves represent the total map, while the orange and green ones represent respectively the reference and estimated primary maps.

IV. DISCUSSION

This paper presents a scatter correction method adapted to spectral CT which uses a primary modulator mask. It has been tested on simulated and experimental data, and shows a good qualitative accuracy in terms of contrast enhancement and scatter artifacts removal. In addition, the analysis of the mean relative error in the simulated case reveals a good quantitative accuracy since the error after correction is close to the error induced by the statistical noise (third vs second column of Table IV). Furthermore, we have observed a good correspondence between the reference and estimated primary in the experimental cases.

The following section aims at discussing the design of the primary modulator mask, noise and ring artifacts increase induced by the method and the choice of the parameters for the B-spline scatter model and the cost function. The choice of the mask material for such scatter correction methods has been studied in [START_REF] Gao | Modulator design for x-ray scatter correction using primary modulation: Material selection[END_REF] considering energy-integrating detectors. It states that to compensate the beam-hardening effect, the material must be chosen such that its K-edge is near the mean energy of the source spectrum. We have chosen a mask material which LAC is slowly-variant along the energy range of the source spectrum because a K-edge would require having a different model for the correction matrix for the rows affected by the K-edge. Low-Z materials such as carbon or aluminum are the best choice for the mask material because their LAC is quite smooth over the considered energy range (30-120 keV).

Having selected the material of the primary modulator mask, the choice of its thickness follows a trade-off between the requirements to be thin enough to limit the parallax effect, which tends to lower the sharpness of the mask edges seen in the images, and sufficiently thick to induce a visible footprint on the estimated primary map when the estimation of the scatter map is wrong.

Finally, the mask must have as many edges as possible in order to have enough pixels carrying information about the presence of the mask footprint. As a consequence, the pattern should be as small as possible given the mechanical constraints of the mask material for the selected thickness. In our physical experiments, the pattern of the primary modulator mask was the smallest possible for the 5 mm thick graphite The color code is the same for all graphs: total map (blue), reference primary map (orange) and estimated primary map (green).

slab according to its manufacturer.

In the presented method, the correction matrix is calibrated for given source, mask and detector positions. In practice, they may be each subject to different small movements during the rotation due to, e.g., vibrations and gravity, which would compromise the calibration of the correction matrix. This will strongly depend on the design of the cone-beam CT scanner which was out of the scope of this study.

b) Noise and ring artifacts: The method does not induce ring artifacts but it increases them as well as the noise level on the images. Due to the smooth behavior of the proposed scatter model, by subtracting it to scatter-corrupted images, the absolute noise level in the transmission sinograms is unchanged, but its relative level increases, particularly after taking the logarithm. As a result, the noise and the vertical stripes are more visible in the estimated primary attenuation sinograms than in the total ones. This is also the case in CT images where the vertical stripes are replaced by ring artifacts.

The observed ring artifacts were more severe without collimation than with collimation, even before scatter correction and without modulation mask (Figure 9, top vs middle). Scatter correction further increases the contrast of those rings. These rings are mostly visible at the edges of the detector crystals (four per Detection Technology ME100 detector) but we do not have an explanation for this physical effect probably due to having a larger part of the detector irradiated.

The scatter-corrected spectral projection may be used in a material decomposition algorithm [START_REF] Alvarez | Energy-selective reconstructions in X-ray computerized tomography[END_REF]. Material decomposition is known to be an ill-conditioned problem on which scatter has a strong impact [START_REF] Sossin | Influence of scattering on material quantification using multi-energy xray imaging[END_REF]. Residual noise and rings after scatter correction may also be a problem, in which case the acquisition parameters can be adjusted and denoising algorithms used prior to the material decomposition. High energy profiles (bin 6)

Fig. 11. Top left: CT slice of RANDO at bin 2 of total map on the left against the estimated primary map on the right. Top right: mean spectra of ROI (location on CT slice). Bottom left and bottom right: profiles at bin 2 and bin 6 respectively. The color code is the same for all graphs: total map (blue), reference primary map (orange) and estimated primary map (green).

c) Choice of the parameters: The spacing between two Bspline knots has to be the largest one allowing recovery of a scatter map which fits well with the real (or simulated) one to minimize the number of parameters. In this work, we performed a grid-search with the objective of minimizing the squared difference between the reference and estimated primary images. The variations of the scatter map (in the two directions of the detector) is highly dependent on the acquisition geometry as it lowers when the air-gap between object and the detector increases. Highly-attenuating objects also increase the high frequencies of the scatter map, and the spacing between two knots on the detector should be smaller when the attenuation of the object increases. This justifies the smaller spacing between B-spline knots of the CIRS pelvis phantom which is much more attenuating than the RANDO head phantom.

Aiming at a smooth representation of the scatter map, the B-splines should be at least quadratic and no improvement has been observed using cubic B-splines. Note that the zeroth order case induces a piece-wise constant representation of the scatter map, as proposed in [START_REF] Ritschl | Robust primary modulation-based scatter estimation for cone-beam CT[END_REF] and [START_REF] Chen | Optimization-based scatter estimation using primary modulation for computed tomography[END_REF]. Therefore, the representation of scatter with B-spline basis functions can be seen as a generalization of previous works.

The cost function presented in this paper in Equation 25 depends on scalar parameters which have to be tuned. The strengths of the structure weighting and the 1 norm smoothing are controlled by η and σ respectively, and their optimal values may be found with a grid-search. It was done in this work using the simulated primary map and the reference primary map (for simulated and experimental cases respectively). Then, the parameters of the Tikhonov regularization (one for each energy bin) were tuned manually. As the problem is poorly conditioned, the optimal solution might have too high magni-tude and lead to negative transmission values on the estimated primary. The selected regularization penalizes large scatter estimates to avoid such negative values.

In this work, the spatial gradient of the primary map was visually preserved (e.g. in Figure 7) using a dedicated weighting in the cost function (Equation 24). Since it is exponential, small variations of the parameter η induce great variations of dynamics in the weighting matrix W 2 and its value has therefore to be tuned quite precisely. On the contrary, the influence of the smoothing parameter σ is not very sensitive and small variations of its value does not visually affect the results provided by the method. The sensitivity of the method to the Tikhonov regularization parameters depends on the condtionning of the method which is related to the statistical noise, the mask attenuation, and the ratio between the number of pixels at mask edges and the number of knots in the Bspline scatter model. Indeed, as explained in section II-A, an error ε = ss induces an error Cε on the estimated primary. The scatter can therefore be detected if ∇Cε is larger than the gradient of the statistical noise and the weight associated with a given B-spline knot can be estimated precisely if the mask attenuation and the number of edge pixels in the area of this knot are large enough. Since the impact of the scatter bias on the attenuation images depends on the SPR, one can conclude that the method is sensitive to the Tikhonov parameters in high SPR regions and where the statistical noise is large. In our experiments, these parameters had to be tuned precisely in the lower energy-bins, where the SPR and the statistical noise are larger.

V. CONCLUSION

In this paper, we presented a scatter reduction technique for spectral CT using a primary modulator mask. The accuracy of the technique has been demonstrated on simulated and experimental spectral CT datasets, both in the projection and in the object domains. The initial mean relative error of around 20% in the lower energy-bins has been reduced to around 4%, which is close to the error induced by the statistical noise.

Fig. 1 .

 1 Fig. 1. Top: schematic of the CT scanner geometry ; bottom: experimental set-up.

Fig. 2 .

 2 Fig. 2. Left: picture of the primary modulator mask; right: details of the mask geometry.

Fig. 3 .

 3 Fig. 3. Left: example of a correction matrix for a pixel attenuated by the primary modulator mask; right: second row's profile ([c 0 , c1 , c2 ] = [1.35, 0.0952, 1.71]).

Fig. 4 .

 4 Fig. 4. Left: simulated detector response matrix (DRM); right: source spectrum.

Fig. 5 .

 5 Fig. 5. Left: Picture of the CIRS phantom; right: picture of the RANDO phantom.

Figure 6

 6 Figure 6 shows the results in the projection domain. The topleft sub-figure is divided into two parts: the left one represents the simulated total attenuation while the right part shows the estimated primary. In the same way, the left part of the topright sub-figure shows a projection of the simulated scatter map and the right part the estimated scatter map. The two bottom sub-figures show the spectra of the two pixels of interest (POI) shown on the top-left sub-figure. Blue curves represent the total attenuation, while the orange and green ones represent the simulated and estimated primary, respectively.

Fig. 6 .

 6 Fig.6. Top: attenuation image (left) and scatter image (right) for one projection angle; both images represent the second energy bin. Each image is subdivided in two: the simulation (resp. total or scatter) on the left against the estimation (resp. primary or scatter) on the right. Bottom graphs share the same color code: total (blue), reference (orange) and estimated (green) spectra of the two pixels of interest (POI).

Fig. 7 .

 7 Fig.7. Profiles of total (blue), reference primary (orange), estimated primary (green), scatter (red) and estimated scatter (magenta) with the same color code in each graph. Left: bin 2, right: bin 6 (profile location on Figure6).

Fig. 9 .

 9 Fig. 9. Left: sinograms of CIRS at bin 2 (left) against bin 6 (right). Right: CT slices of CIRS at bin 2 (left) againt bin 6 (right). Top: total; center: reference primary; bottom: estimated primary.

Fig. 10 .

 10 Fig. 10. Top left: location of regions of interest (ROI); Top right, bottom left, bottom right: mean spectra of ROI 1, ROI 2 and ROI 3 respectively.The color code is the same for all graphs: total map (blue), reference primary map (orange) and estimated primary map (green).

TABLE IV MEAN

 IV RELATIVE ERRORS ON SIMULATED DATA (%).

	Bin # M RE t b	M RE p b	M RE p b
	1	18.55	2.64	4.08
	2	21.27	2.06	3.45
	3	12.95	1.37	2.70
	4	9.26	1.24	2.36
	5	7.42	1.43	2.46
	6	6.21	1.63	2.56
	7	5.30	2.05	3.23
	8	4.81	3.21	4.02

ACKNOWLEDGEMENTS

The authors would like to thank Michel Garcin for his help provided for the realization of the physical experiments.

was supported by the LABEX PRIMES (ANR-11-LABX-0063) of Universit é de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX