N Antoni

J.-L Ligier

P Saffre

Cumulative microslip at the interface of mechanical assemblies under cyclic loading

Mechanical contact problems involving friction are widely encountered in many industrial applications. Despite the recent progress in non-smooth mechanics and tribology, a number of these problems are still open as they cannot be properly simulated with current numerical tools. The cumulative microslip phenomenon at the interface of mechanical assemblies subject to thermal or mechanical cyclic loading may represent one of these issues. In practice very few mechanical components are designed to prevent that phenomenon from appearing and it is possible to observe significant relative displacements at the interface of affected assemblies leading to a dysfunction.

Introduction

The cumulative microslip phenomenon occurs in mechanical assemblies subject to thermal or mechanical cyclic loading and it is characterized by the accumulation at the contact interface of tangential microslip of one element relative to another assembly element in a preferred direction. These last features give it a determinist nature which resembles the Ratchet effect in plasticity (see for example [START_REF] Lemaitre | Mécanique des Matériaux Solides[END_REF]). In practice, such cumulated microslip may cause the failure of some part assemblies in reciprocating engines with internal combustion (see [START_REF] Ligier J-L | Analyses avancées et expérimentales[END_REF]). Some of them are particularly subject to this phenomenon ; this is the case of the connecting rod big end associated with its bearing and the connecting rod small end associated with its bushing, for which the cumulative microslip phenomenon has been found fairly often.

For a case of cumulative microslip phenomenon induced by asymmetric friction conditions at the contact interface of mechanical assemblies submitted to uniform thermal loading cycles, we show, in this article, some of the results obtained. First, due to the quasilack of literature on this subject, it was necessary to underscore theoretically the existence of such a phenomenon, well-known in practice. For this, a discrete analytical model of a simplified assembly involving a suitable friction law is proposed. This is special in that it extends Coulomb's law of friction which, in its classical form, does not allow a reproduction of this phenomenon. The analysis of this assembly behaviour under thermal cyclic loading show under which conditions the cumulative microslip phenomenon may occur. Secondly, we analyze an anisotropic friction model that allows the determination of the law typical parameters and shows that the arbitrary assumptions previously used for the coefficient of friction are realistic. As a friction law may not be ordinary, we have bought a highly special attention to demonstrating its thermodynamic admissibility.

Theoretical underscoring

For the sake of simplicity, we shall first concentrate our attention on a discrete model to present the phenomenon and we shall extend it with a continuous one.

One-dimensional discrete model

Let us consider the discrete model of a simplified assembly made of a contact elastic spring (coefficient of linear thermal expansion T α and original length 0 l ) with longitudinal stiffness ( ) k , initially at the reference temperature 0 T , for which both nodes are kept in contact on a fixed rigid support under the action of a constant normal nodal force 0 > N . The degree of freedom in translation of node i is denoted as i u ( )

2 , 1 = i
, see Figure 1. The displacement of each contacting node is obviously controlled by friction conditions at the contact interface. Let i t F be the nodal tangential force applied by node i on the support ( )

2 , 1 = i
. To include frictional effects let us assume that regularized asymmetric Coulomb's law of friction holds pointwise on the contact interface (see Figure 2). It takes two phases into account:

A reversible phase corresponding to the elastic deformation of asperities at the contact interface. Also mentioned is the adherence phase with elastic microslip. Let t k denote the contact tangential stiffness which is the typical parameter of this phase ;

An asymmetric irreversible phase corresponding to an asymmetric macroscopic slip : the value of the friction coefficient is different depending on the direction of slip considered. It is assumed, for the one-dimensional model suggested here, that the value of the friction coefficient in the 0 > x -direction is higher than the one in the 0 < x -direction, i.e. :

-+ > µ µ

(1)

u 1 k N ∆T x u 2 N 1 2
Figure 2 : Regularized asymmetric Coulomb's law of friction.

Let T be the uniform temperature forced by the outside world onto the system. Next, we define the temperature difference 0

T T T - = ∆
. Within the framework of small displacements (implying small slip), analytic developments of this assembly response under quasi-static thermal cyclic loading

[ ] max min , T T T ∆ ∆ ∈ ∆
show that, if condition (1) holds, the cumulative microslip phenomenon occurs as soon as

c T T T ∆ > ∆ - ∆ 2
min max (see Figure 3). In such an inequality, c T ∆ denotes a critical temperature variation corresponding to the contact tangential elasticity limit of the system and is defined as :

        + = ∆ - t T c k k l N T 2 1 0 α µ (2) Figure 3 : Cumulated microslip in the 0 < x -direction.
The cumulative microslip phenomenon is characterized by the accumulation of small displacements for each node in the 0 < x -direction and for each loading cycle, i.e. :

   < ∆ < ∆ 0 0 2 1 u u ( 3
)
where i u ∆ is the incremental nodal microslip of node i for each cycle ( )

2 , 1 = i .
With such assumptions, it can be showed that the incremental nodal microslip is identical for both nodes, constant during (stabilized) cycles and strictly negative, as :

( ) 0 2 max min 0 2 1 < ∆ + ∆ - ∆ = ∆ = ∆ c T T T T l u u α (4) max T ∆ min T ∆ c T ∆ c T T ∆ - ∆ 2 max 0 ∆u 1 ∆u 2 T ∆ i t F N + µ N - -µ i u 0 t k u 1, u 2 c T T ∆ + ∆ 2 min

One-dimensional continuous model

There should be now considered the continuous equivalent of the previous discrete model. A thin ( )

l b h << , elastic ( E is Young's modulus, T
α is the thermal expansion coefficient) strip kept in contact on a fixed rigid support, under the action of a uniformly distributed normal surface traction p (see Figure 4), may be assimilated, for a tension only behaviour, with a serial assembly comprising an infinity of such discrete systems.

Similarly, each particle on the contact interface 0

Ω ∂ C
is assumed to obey regularized asymmetric Coulomb's law of friction previously introduced.

Analytic developments enable again showing the occurrence of the phenomenon under thermal cyclic loading

[ ] max min , T T T ∆ ∆ ∈ ∆ if condition (1) holds and if - ∆ > ∆ - ∆ c T T T 2 min max . In the previous inequality, - ∆ c
T denotes a first critical temperature variation corresponding to the original tangential contact elastic limit of the system and is defined by:

        -                 = ∆ - - 1 ch sh Eh k l Ehk Eh k l p T t t T t c α µ (5)
When this is reached, the end 0 = x starts slipping in the 0 < x -direction whereas the rest of the strip is in the adherence phase with elastic microslip. The cumulative microslip phenomenon is characterized by the accumulation in the 0 < x -direction of small displacements of the strip for each particle on the contact interface 0 Ω ∂ C and for each loading cycle, i.e. :

( )

0 0 . Ω ∂ ∈ ∀ < ∆ C x x e x u r (6)
where u r ∆ is the incremental small displacement for each cycle. Results have been extended to the asymmetric case of two fitted tubes (one being fixed) but the analytical expression is too complicated to be presented in this paper.

At this stage, the cumulative microslip phenomenon may occur when the contact interface behaviour is governed by asymmetric friction conditions. The physical validity of such assumptions is given in the following paragraph.

Sufficient condition for existence

In order to justify the physical admissibility of regularized asymmetric Coulomb's law of friction we have proceeded with two steps. First, a natural approach based on micromechanical observations provide an anisotropic friction model enabling, via averaging, a determination of the law characteristics (in this respect see [START_REF] Mroz | An Anisotropic Friction and Wear Model[END_REF], [START_REF] Mroz | Contact Friction Models and Stability Problems[END_REF]). Then its thermodynamic admissibility is proven.

Anisotropic friction model

Let us consider the idealized surface texture (in practice, there are always dispersions, see Figure 5a) of a fixed rigid body denoted as (2) on which are periodically distributed long wedge-shaped asperities (see Figure 5b). Each of them is characterized by geometrical properties ( ) h , , β α as well as mechanical ( )

β α µ µ ,
. This last couple represents the real friction coefficients for each wedge of one asperity (see Figure 5c). The steady macroscopic sliding of body (1) over body (2) is then considered, assuming first an infinite flexural stiffness : the springs contract and release while sliding along wedge-shaped asperities so that the body (1) motion only occur in the nominal contact plane ( ) b a t t , ˆ. It is then possible to express, averaging over a period, the average friction coefficient in the principle direction a t ˆ of the contact interface. Its expression varies according to the considered path direction and it is independent from n k and from h : [START_REF] Lemaitre | Mécanique des Matériaux Solides[END_REF] where the superscripts "+" and "-" denote the a t ˆ and ( ) a tdirections, respectively, and the function H is defined as :

( ) ( ) ( ) ( )        + - + = + - + = - + β α β α µ α β µ µ β α α β µ β α µ µ α β β α tan tan tan , tan , tan tan tan , tan , H H H H a a
( ) y x y x y x H tan 1 tan , - + = (8)
From equations [START_REF] Lemaitre | Mécanique des Matériaux Solides[END_REF], it is shown that it is possible to reach a macroscopic friction asymmetry in the case of a complete asymmetry of the body (2) surface texture (geometrical and mechanical properties), i.e. :

( ) ( ) - + ≥ ⇒ ≥ ≥ a a µ µ µ µ β α β α and (9)
One of the features of equations ( 7) is that they do not depend on the height h of wedge asperities. So, they remain valid for a body (2) surface texture that is physically more realistic, with asperities of a constant shape (i.e. ( ) β α , constant) but with a Gaussian distribution of their heights (this is the case met most often in practice but other distributions may also be considered, for instance Normal log. distribution). Note that the height h for each wedge asperity must be measured from the same mean line.

An infinite flexural stiffness has been assumed and thus a rigid body motion of the springs during the steady sliding of body (1). If this stiffness is assumed to be large but finite, there is then a bending of the springs prior to their sliding along wedge asperities indeed in relation with a linear elastic behaviour which allows a reversible microslip e t a u of body (1) relative to body (2) in the principle slip direction a t ˆ of the nominal contact plane.

The relation between principal tangential stress a t σ at the nominal contact interface and this reversible displacement is thus linear, i.e. : One generalizes on the one hand equation ( 10) from a two-dimensional point of view, assuming an isotropic bending behaviour of the springs, which then becomes :

t t e t k u σ r r = ( 11
)
The springs start sliding in an elastically deformed state whenever the macro slip condition has been satisfied, which allows on the other hand a justification of the wellknown decomposition (see [2], [START_REF] Curnier | A Theory of Friction[END_REF], [START_REF] Klarbring | Derivation and Analysis of Rate Boundary-Value Problems of Frictional Contact[END_REF], [START_REF] Mroz | Contact Friction Models and Stability Problems[END_REF] and [START_REF] Strömberg | Derivation and Analysis of a Generalized Standard Model for Contact, Friction and Wear[END_REF] for example) of the total relative tangential displacement (observable), into one reversible elastic part denoted as e t u r (however small it is), and one irreversible inelastic part denoted as in t u r and corresponding to the macro slipping :

in t e t t u u u r r r + = (12)
Equation ( 11) is always adopted in the literature but is not basically very realistic if whatever a surface texture is taken into account. A reversible orthotropic phase could be thus considered, characterized by two contact tangential stiffness

Thermodynamic admissibility

Generally speaking, a mechanical constitutive law cannot be any, and it must comply with some principles for it to have a physical sense. It must, in particular, be such that the two fundamental principles of thermodynamics, thus the Clausius-Duhem basic inequality, be satisfied. Starting from experimental observations, a method based on the concept of generalized standard materials is applied in order to rebuild a law valid from a thermodynamic point of view. Such approach is commonly used to derive mechanical constitutive models for material non-linearities (see [START_REF] Lemaitre | Mécanique des Matériaux Solides[END_REF] for example). A friction law may be considered as a mechanical constitutive law for a of continuous medium region that is not of volume any more and only existing in a contact situation. The general and natural framework for writing friction laws is thus that of continuum thermodynamics.

Since the method previously quoted may not be directly particularized for contact interface, it must therefore be adapted. An implementation has been proposed by several authors. There may be quoted, among others, Moreau (see [START_REF] Moreau | Sur les Lois de Frottement, Plasticité et Viscosité[END_REF]) for his work on non-smooth phenomena like plasticity and friction, Felder (see [START_REF] Felder | Formulation Thermodynamique des Interactions Superficielles entre deux Corps[END_REF]) for his work on surface interactions applied to metal forming, or also Strömberg et al (see [START_REF] Strömberg | Derivation and Analysis of a Generalized Standard Model for Contact, Friction and Wear[END_REF]) for studying fretting (wear process arising when contacting surfaces undergo oscillatory small displacements).

The general approach consists first in using the Clausius-Duhem inequality deducted from the expressions in local form of the first and second principles of thermodynamics at the contact interface. In order to prove the thermodynamic validity of regularized asymmetric Coulomb's law of friction, the three following steps are taken successively into consideration :

Choice of a set of independent state variables, internal representation of the state ; Derivation of state laws for reversible interactions from a specific free energy ; Derivation of evolution laws for irreversible interactions from a specific dissipation potential which must be taken as convex (but not necessarily differentiable) in order to ensure satisfaction of the dissipation inequality.

Set of independent state variables :

The following elementary set is selected :

( ) in t t C u u r r , = χ (13)
These two state variables (observable and internal, respectively) are assumed to comply with the decomposition in equation [START_REF] Oden | Nonlocal and Nonlinear Friction Laws and Variational Principles for Contact Problems in Elasticiy[END_REF]. The associated thermodynamic forces are respectively ( )

t t C Y σ σ r r , =
.

Specific free energy and state laws :

The thermodynamic potential to be taken into account is written as:

( ) ( ) 2 2 1 , in t t t in t t C u u k u u r r r r - = ψ (14)
Taking into account the decomposition in equation ( 12), the state law which derives from this thermodynamic potential is given by:

e t t t u k r r = σ ( 15
)
which is simply the reversible phase in equation [START_REF] Mroz | Contact Friction Models and Stability Problems[END_REF].

Specific dual pseudo-potential and evolution laws :

Let us consider the closed convex set 2 R ⊂ K
as follows:

( ) { } 0 , / 2 ≤ ∈ = n t t g K σ σ σ r r R ( 16
)
where g is a convex differentiable function, so-called friction criterion. The proposal dual pseudo-potential, from which derive the complementary laws, is defined as :

( ) ( ) t K t C I σ σ r r = Φ * (17)
where K I denotes the indicator function of the set K .

As this is convex, K I is also convex. This function is also minimum at 0 r r = t σ . The pseudo-potential * Φ C defined above being convex with respect to the dual variable ( ) t σ r and minimum at ( )

0 r r = t σ
, the Clausius-Duhem inequality is therefore satisfied.

Several choices of g are possible. One friction limit criterion, differentiable and convex, directly issued from the previously described anisotropic friction model, is defined by : , the evolution law that derives from the convex (but not differentiable) dual pseudo-potential * Φ C is expressed as : A sliding rule, associated with the friction limit criterion, has been selected here. It is however possible to introduce a family of sliding potentials that are not associated with this criterion (in this respect, see [START_REF] Mroz | Contact Friction Models and Stability Problems[END_REF]).

( )                  ≤ ≥ -                 +         ≤ ≤ -                 +         ≥ ≤ -                 +         ≥ ≥ -                 +         = - + - - + - + + 0 et 0 if 0 et 0 if 0 et 0 if 0 et 0 if , , 2 / 1 2 2 2 / 1 2 2 2 / 1 2 2 2 / 1 2 2 b a b a b a b a b a b a b a b a b a t t

Conclusion

Asymmetric friction conditions at the contact interface which characterize regularized asymmetric Coulomb's law of friction and represent a sufficient condition to obtain the cumulative microslip phenomenon have been achieved from micromechanical considerations on surface textures. This law, for which the thermodynamic validity has been proven, is thus associated with a given surface texture. In practice, we must underline that the available surface texture parameters in the industrial tools do not accurately enable the detection and characterization of surface textures of the type presented in this article. Also, usual FEA do not allow the direct integration of regularized asymmetric Coulomb's law of friction whereas it should be easy for a given situation, to define asymmetry friction levels to be avoided so the phenomenon does not occur.

A comparison of these results with experiments is planned in order to validate this type of approach in spite of the simplifying assumptions used.

Figure 1 :

 1 Figure 1 : One-dimensional discrete model.

Figure 4 :

 4 Figure 4 : One-dimensional continuous model.The calculated response of the strip under thermal cyclic loading clearly exhibits, in a continuous sense, all the expected similarities with the one of the previous discrete model. The expression for the incremental microslip is not given here due to its complexity.

Figure 5 :

 5 Figure 5 : A surface texture model.An elastic body denoted as (1) on which are isotropically and uniformly distributed asperities modelled as identical elastic springs (i.e. same original length, same longitudinal stiffness n k and same flexural stiffness f k ), is kept in contact on the previous one (see Figure6).

Figure 6 :

 6 Figure 6 : A type of contact model.

  principle slip directions a t ˆ and b t ˆ.

Figure 7 :

 7 Figure 7 : Anisotropic friction limit criterion: example of representation.

& & & &(19)

Acknowledgements -This work is supported by Renault S.A.. Pr. J. Pastor, responsible for MMSC research team, is also gratefully acknowledged here.

A write-up that is so general assumes the existence of a surface texture exhibiting such macroscopic friction properties, which is not obvious. The criterion expression has been built from that of the usual orthotropic friction limit criteria, shown as ellipses (see [1], [START_REF] Felder | Thermodynamics of Surface Interactions between two Media with Special Emphasis on Anisotropic Effects[END_REF], [START_REF] Mroz | Contact Friction Models and Stability Problems[END_REF], [START_REF] Panagiotopoulos | Hemivariational Inequalities in Frictional Contact Problems and Applications[END_REF]).

An example of representation ( )

with that of a Coulomb's isotropic friction criterion ( )

ABAQUS THEORY MANUAL, Surface Interactions, Version 6.2, Hibbit, Karlsson & Sornsen, Inc., U.S.A., pp. 5.2.3-1 -5.2.6-3, 2001.

[2] CHENG J-H., KIKUCHI N., An Incremental Constitutive Relation of Unilateral


Cyril Mory

Bruno Sixou

Salim Si-Mohamed

Loïc Boussel

Simon Rit

Salim Si

Comparison of five one-step reconstruction algorithms for spectral CT

 several hundreds, and Barber2016 [3] several thousands. We also sum up other practical aspects, like memory footprint and the need to tune extra parameters.

Introduction

Data acquisition

Dual energy computed tomography (CT) systems are now commercially available, and spectral photon counting CT systems already exist as research prototypes. These new types of scanners provide information on the energy distribution of the X-ray photons that reach the detector: dual energy systems acquire two sets of projections (for low and high energies), while spectral photon counting systems group the incoming photons into so-called "energy bins", which form a partition of the incident spectrum's bandwidth. For each projection's pixel, the "raw" measurements acquired by a spectral photon counting CT scanner are therefore the number of photons in each energy bin, which are commonly referred to as "photon counts".

Attenuation model

Assuming that the object's attenuation can be obtained by linear combination of the attenuations of only a few materials, this energy-resolved information allows to reconstruct multiple volumes, each one representing a different material's concentration map. Alternatively, one can assume that the attenuation stems from only a few physical phenomena, e.g. Compton scattering, photo-electric effect and K-edge [START_REF] Alvarez | Energy-selective reconstructions in X-ray computerized tomography[END_REF][START_REF] Pan | Computed Tomography in Color: NanoK-Enhanced Spectral[END_REF], and reconstruct maps of the characteristic quantities of each of these physical phenomena. For the sake of readability, we present the case of spectral photon-counting CT, and decompose on a basis of materials, but all calculations are also valid for dual-energy CT, and with photo-electric effect and Compton scattering as decomposition basis. It is assumed that the linear attenuation coefficient μje of the imaged object in voxel j and at energy e can be separated into a function of space and a function of energy, i.e. μje = Nm m=1 µ me x jm [START_REF] Alvarez | Energy-selective reconstructions in X-ray computerized tomography[END_REF] where N m is the number of materials, µ me the mass-attenuation coefficient of material m at energy e, and x jm the concentration of material m in voxel j. The attenuation profiles (µ me ) e=1..150keV are known, and the aim is to find the concentration maps x m .

Reconstruction methods

The most widespread methods to reconstruct material-specific volumes from photon counts can be divided into two categories: image-based and projection-based. In two-step image-based methods, each energy bin of the photon counts is log-transformed and reconstructed, just like the projections dataset of a standard CT. This process yields one volume per energy bin and these intermediate volumes are then decomposed into material-specific volumes [START_REF] Maass | Image-based dual energy CT using optimized precorrection functions: A practical new approach of material decomposition in image domain[END_REF][START_REF] Niu | Iterative image-domain decomposition for dualenergy[END_REF]. Unfortunately, the intermediate volumes are corrupted by beam hardening artefacts [START_REF] Taguchi | X-map 2.0" and "iBHC2" for consistent Edema Signal Enhancement for Acute Ischemic Stroke using Non-Contrast-Enhanced Dual-Energy CT Salt Lake City[END_REF], because the energy dependence of the attenuation coefficients is still averaged in each energy bin. Image-based methods often integrate some form of empirical beam-hardening correction [START_REF] Maass | Image-based dual energy CT using optimized precorrection functions: A practical new approach of material decomposition in image domain[END_REF], but perfect beam-hardening correction would require knowledge of the material volumes. In two-step projection-based methods, on the other hand, multichannel projections are first decomposed into material-specific projections, and then reconstructed, independently [START_REF] Alvarez | Energy-selective reconstructions in X-ray computerized tomography[END_REF][START_REF] Schlomka | Experimental feasibility of multi-energy photon-counting K-edge imaging in pre-clinical computed tomography[END_REF] or jointly [START_REF] Sawatzky | Proximal ADMM for Multi-Channel Image Reconstruction in Spectral X-ray CT[END_REF][START_REF] Mory | Generalized least squares for spectral and dual energy CT: A simulation study[END_REF]. It turns out that for typical choices of materials (e.g. water, bone, and a high-Z contrast agent), the normalized attenuation profiles of the materials do not differ enough to yield a robust decomposition: the inevitable statistical noise on photon counts is often sufficient to cause aberrant material line integrals, which result in strong streak artefacts in the reconstructed volumes. Some methods were proposed to regularize the decomposition process [START_REF] Brendel | Penalized Likelihood Decomposition for Dual Layer Spectral CT in[END_REF][START_REF] Ducros | Regularization of Nonlinear Decomposition of Spectral X-ray Projection Images Medical Physics[END_REF], but how to regularize in the projection domain is a difficult problem: there are conditions for projections to be consistent, i.e. to be those of an actual object, called Data-Consistency Conditions (DCC) [START_REF] Clackdoyle | Data Consistency Conditions for Cone-Beam Projections on a Circular Trajectory[END_REF][START_REF] Lesaint | Calibration for Circular Cone-Beam CT Based on Consistency Conditions[END_REF], and regularization with priors used routinely for volumes (e.g. total variation) can generate inconsistent projections. In addition, projection-based methods require the spectral projections to have the same geometry (same source and detector positions), it cannot be applied to fast kV switching [START_REF] Lee | Dual-Energy CT in Cardiovascular Imaging[END_REF] or dual source dual-energy systems [START_REF] Flohr | First performance evaluation of a dual-source CT[END_REF] easily: corresponding projections must be interpolated, which is a potential source of errors. Two-step methods also have a structural drawback: unless the first step provides a one-to-one mapping between inputs and outputs, it implies a loss of information, for which the second step cannot compensate. The aberrant pixels observed after projection-based decomposition illustrate that loss of information: the reconstruction step cannot go back to the photon counts, and therefore must rely on aberrant decomposition results. Recently, several methods have been proposed which reconstruct material-specific volumes directly from the photon counts [START_REF] Cai | A full-spectral Bayesian reconstruction approach based on the material decomposition model applied in dual-energy computed tomography[END_REF][START_REF] Long | Multi-Material Decomposition Using Statistical Image Reconstruction for Spectral CT[END_REF][START_REF] Weidinger | Polychromatic Iterative Statistical Material Image Reconstruction for Photon-Counting Computed[END_REF][START_REF] Barber R F | An algorithm for constrained one-step inversion of spectral CT data[END_REF][START_REF] Mechlem | Joint Statistical Iterative Material Image Reconstruction for Spectral Computed Tomography Using a Semi-Empirical Forward Model[END_REF][START_REF] Chen | Image reconstruction and scan configurations enabled by optimization-based algorithms in multispectral[END_REF]. They are commonly referred to as "one-step inversion", or simply "one-step" methods. All of these methods are iterative: there is currently no analytical inversion formula for the material decomposition problem, let alone for one-step inversion. They consist in combining the forward models of the tomographic reconstruction and the material separation inverse problems, yielding a single (but more complex) forward model, which takes as input a set of material volumes, and yields photon counts as output. In discrete form, the forward model reads

ȳ ib = e s be exp - j a ij m µ me x jm ( 2
)
where ȳib is the expected photon count in pixel i and bin b, s be is the effective spectrum, i.e. the number of photons of actual energy e that are expected to be detected in bin b if there was no attenuating object and a ij is the coefficient of the forward projection matrix A at row i and column j. One-step inversion then means finding the volumes x m for which the forward model yields photon counts y as similar as possible to those measured by the scanner. By construction, the drawbacks of two-step methods are circumvented, and there is no risk of losing information. One-step methods have been reported to yield improved image quality [START_REF] Zhang | Model-Based Iterative Reconstruction for Dual-Energy X-Ray CT Using a Joint Quadratic Likelihood Model[END_REF] and reduced noise-induced bias [START_REF] Mechlem | Joint Statistical Iterative Material Image Reconstruction for Spectral Computed Tomography Using a Semi-Empirical Forward Model[END_REF] with respect to their two-step counterparts, but they are also slower to converge. This is likely to prevent their use in many applications, as it happened for regularized iterative reconstruction techniques on standard CT [START_REF] Pan | Why do commercial CT scanners still employ traditional, filtered back-projection for image reconstruction[END_REF]. Since in the multi-energy CT case the amount of data contained in a single acquisition is multiplied by the number of energy bins, and the one-step problem is more complex than the reconstruction problem alone, the need for fast methods is even more acute than for standard CT. Among the methods studied in this paper [START_REF] Cai | A full-spectral Bayesian reconstruction approach based on the material decomposition model applied in dual-energy computed tomography[END_REF][START_REF] Long | Multi-Material Decomposition Using Statistical Image Reconstruction for Spectral CT[END_REF][START_REF] Weidinger | Polychromatic Iterative Statistical Material Image Reconstruction for Photon-Counting Computed[END_REF][START_REF] Mechlem | Joint Statistical Iterative Material Image Reconstruction for Spectral Computed Tomography Using a Semi-Empirical Forward Model[END_REF][START_REF] Barber R F | An algorithm for constrained one-step inversion of spectral CT data[END_REF], only Cai et al. [START_REF] Cai | A full-spectral Bayesian reconstruction approach based on the material decomposition model applied in dual-energy computed tomography[END_REF] reported a comparison with another one-step method [START_REF] Sukovic | Penalized weighted least-squares image reconstruction for dual energy X-ray transmission tomography[END_REF] and an advanced image-based two-step method [START_REF] Maass | Empirical multiple energy calibration[END_REF]. Barber et al. [START_REF] Barber R F | An algorithm for constrained one-step inversion of spectral CT data[END_REF] showed many results, but only with the proposed method, and the other ones [START_REF] Long | Multi-Material Decomposition Using Statistical Image Reconstruction for Spectral CT[END_REF][START_REF] Weidinger | Polychromatic Iterative Statistical Material Image Reconstruction for Photon-Counting Computed[END_REF][START_REF] Mechlem | Joint Statistical Iterative Material Image Reconstruction for Spectral Computed Tomography Using a Semi-Empirical Forward Model[END_REF] are compared to two-step methods, either projection-based or image-based [START_REF] Mendonca P R S | A Flexible Method for Multi-Material Decomposition of Dual-Energy CT Images[END_REF], with the exception that Mechlem et al. [START_REF] Mechlem | Joint Statistical Iterative Material Image Reconstruction for Spectral Computed Tomography Using a Semi-Empirical Forward Model[END_REF] contains a comparison with the method of Weidinger et al. [START_REF] Weidinger | Polychromatic Iterative Statistical Material Image Reconstruction for Photon-Counting Computed[END_REF] in terms of convergence speed. An independent comparison between many one-step methods is therefore still missing, although it would certainly benefit the spectral CT scanner manufacturers, saving them the hassle of investigating many methods to implement only one. We provide such a comparison in this paper, with experiments conducted on simulated data on a single slice. Since we provide the code we used, it will also help the research community: ranking a new one-step method against the ones studied in this paper will only require to implement the new method in the open source framework we provide.

All methods were first adapted to solve the same one-step problem: spectral photon-counting CT reconstruction from five energy bins using a three materials decomposition basis and spatial regularization. Since the five methods we studied use different data-attachment and regularization terms, as well as different optimization techniques, we defined an image quality objective, independant of the cost functions, and measured the number of iterations required to reach it with each method. These results constitute the main contribution of the present paper. For each method, we also study its memory footprint, state whether it requires to tune additional parameters, and whether or not the derivation of the algorithm uses simplifying assumptions/approximations.

In [START_REF] Barber R F | An algorithm for constrained one-step inversion of spectral CT data[END_REF][START_REF] Fessler | Method for statistically reconstructing images from a plurality of transmission measurements having energy diversity and image reconstructor apparatus utilizing the method[END_REF], it is proposed to perform internal computations not with the actual materials into which one wishes to decompose the object, but with carefully chosen linear combinations of these materials. Since the materials are characterized only by their attenuation coefficient µ as a function of the incident energy, Barber et al. call this step µ-preconditioning. In this paper, we present three different ways to perform µ-preconditioning and evaluate their impact on convergence speed. Since most methods were not explicitly designed for µ-preconditioning, a small adaptation is required (and described in the present paper) for the regularization step.

Material and methods

Compared methods

We have studied five one-step reconstruction methods, adapted them to solve the 3-material, 5-bin, unconstrained and regularized one-step spectral reconstruction problem, and compared their behavior. These methods are:

• Cai2013 [START_REF] Cai | A full-spectral Bayesian reconstruction approach based on the material decomposition model applied in dual-energy computed tomography[END_REF], which uses non-linear conjugate gradient with a heuristical step size,

• Long2014, Weidinger2016 and Mechlem2018 [START_REF] Long | Multi-Material Decomposition Using Statistical Image Reconstruction for Spectral CT[END_REF][START_REF] Weidinger | Polychromatic Iterative Statistical Material Image Reconstruction for Photon-Counting Computed[END_REF][START_REF] Mechlem | Joint Statistical Iterative Material Image Reconstruction for Spectral Computed Tomography Using a Semi-Empirical Forward Model[END_REF], which use Separable Quadratic Surrogates (SQS) to minimize cost functions,

• Barber2016 [START_REF] Barber R F | An algorithm for constrained one-step inversion of spectral CT data[END_REF], which relies on a generalization of the Chambolle-Pock algorithm [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to[END_REF].

All methods use the polychromatic Beer-Lambert law as forward model (the way to estimate the expected photon counts y from material maps x m ) which is reminded in its discrete form in equation 2. They slightly differ by the cost functions they attempt to minimize, which are summarized in Table 1. The data-attachment terms differ because they model the noise differently: Cai2013 optimizes a Gaussian observation model in a Bayesian framework where all the others minimize the Poisson negative loglikelihood. All regularization terms are different but share some similarities, involving the 1-norm (or a twice-differentiable approximation of it) of the spatial gradient of each reconstructed volume.

Since three methods out of five use them, it is useful to remind the definition of Separable Quadratic Surrogates (SQS) and how they are used. Surrogates are defined as follows: let Ψ : R N → R and

x 0 ∈ R N . The function Φ x 0 : R N → R is a surrogate of Ψ at x 0 if and only if ∀x ∈ R N , Φ x 0 (x) ≥ Ψ(x) Φ x 0 (x 0 ) = Ψ(x 0 ), Φ x 0 (x 0 ) = Ψ (x 0 ) (3)
i.e. if Φ x 0 is above Ψ on R N , and tangent to Ψ at x 0 . Note that we have used R N as domain for simplicity, but in Long2014 the domain actually used is R +N . The algorithms based on surrogates consist in finding a surrogate Φ x 0 to the cost function Ψ at the current iterate x 0 , minimizing the surrogate instead of the original cost function, adopting that minimizer as the new iterate x 1 , and starting over. Obviously, the surrogates must be chosen carefully in order to be easier to minimize than the original cost function.

In [START_REF] Long | Multi-Material Decomposition Using Statistical Image Reconstruction for Spectral CT[END_REF][START_REF] Weidinger | Polychromatic Iterative Statistical Material Image Reconstruction for Photon-Counting Computed[END_REF][START_REF] Mechlem | Joint Statistical Iterative Material Image Reconstruction for Spectral Computed Tomography Using a Semi-Empirical Forward Model[END_REF], surrogates are derived one after the other: first a convex one, then a quadratic one, and finally a separable quadratic one, each new one being a surrogate of all the previous ones, and therefore of the original cost function.

Preliminary comments

We have tried to stick as much as possible to the original algorithms, but some modifications were required. In this section, we briefly introduce each method, and describe the modifications we had to perform. Though we have tried to make the present paper as self-contained as possible, it is impossible to explain all the methods in details without ending up duplicating large passages of the original papers.

The detailed maths can be found in our notes (see section 2.7), in which we guide the reader through the equations and re-do most calculations, for all methods but Barber2016. For Barber2016, our notes only provide implementation suggestions, because the original paper is already written in a very didactic way. However, since we had to adapt the five studied methods to deal with a common problem, we provide the pseudo-code of all methods in Appendix. In these pseudo-codes, as much as possible, we used common notations for all methods, in particular for the matrices that describe detector response, incident spectrum and material attenuations. A drawback of this approach is that the pseudo-code notations no longer correspond to those of the original papers. Throughout the paper, when referring to variable names, we use the ones from the pseudo-code. We use a contraction of the first author's name and publication year (e.g. "Cai2013") to refer to each method. When the distinction is relevant, it should be clear from the context whether we mean the method implemented and modified by us, or the one described in the original papers.

2.2.1. Cai2013 Cai2013 [START_REF] Cai | A full-spectral Bayesian reconstruction approach based on the material decomposition model applied in dual-energy computed tomography[END_REF] was originally published for dual energy CT. The cost function is built with a variable y that is the ratios between photon counts with object and photon counts without object (instead of just the photon counts with object, in the other methods). The variance of these ratios y is assumed to be proportional to y, with a factor k d that has to be estimated. We generalized Cai2013 [START_REF] Cai | A full-spectral Bayesian reconstruction approach based on the material decomposition model applied in dual-energy computed tomography[END_REF] to 3 materials, which is straightforward, and k d is determined by a simulation with the same spectrum, detector response and geometry as in the experiments, but without object. It yields k d = 1.59 × 10 -4 . The regularization term is the Huber function of the spatial gradient of the volumes, with the Huber function defined as

φ Huber (x, δ) = x 2 if |x| < δ, 2δ|x| -δ 2 else. ( 4
)
The cost function is minimized by a non-linear conjugate gradient, with a heuristically determined step size. Given the descent direction, the step size is computed as the one that would minimize the cost function's second-order Taylor expansion. Then the cost function is computed at the candidate destination point, and the step size divided by 2 if it causes an increase of the cost instead of a decrease. According to the authors, the case where the step length had to be divided by 2 never occurred in the experiments they conducted (which were limited to dual-energy CT). In ours, it did occur, and we had to implement a reset mechanism: if, after 10 attempts to divide the step length by 2, we still obtain an increase of the cost function, we assume that the descent direction yielded by non-linear CG is erroneous, and try using the opposite gradient instead, as in an ordinary gradient descent. Since it is an ad hoc adaptation, and only for a special case, we have not included it into the pseudo-code in the Appendix.

Long2014

Long2014 [START_REF] Long | Multi-Material Decomposition Using Statistical Image Reconstruction for Spectral CT[END_REF] defines the following hyperbola for approximating the absolute value function in the regularization:

φ Long (x, δ) = δ 2 3 1 + 3 x δ 2 -1 . (5)
The authors sequentially derive the surrogates of the cost function to minimize it: first a convex one, then a quadratic one, and finally a separable quadratic one. Deriving the quadratic surrogate requires to assume that the material maps are non-negative. It was originally published as a method for reconstructing three or more materials from dual energy CT scans. To that end, it adds a lot of a priori information, in the form of equality and inequality constraints on the reconstructed volumes. Performing minimization while complying with all constraints requires an inner loop of quadratic programming using the Generalized Sequential Minimal Optimization (GSMO) algorithm [START_REF] Keerthi | Convergence of a Generalized SMO Algorithm for SVM Classifier Design[END_REF]. Since in our case, there are more bins than the number of materials, we have removed all constraints in our implementation of Long2014, which allowed to replace the quadratic programming inner loop by a single iteration of Newton's method. Long2014 uses Ordered Subsets (OS) to speed up convergence. See section 2.4 for a quick reminder on OS.

Weidinger2016

Weidinger2016 [START_REF] Weidinger | Polychromatic Iterative Statistical Material Image Reconstruction for Photon-Counting Computed[END_REF] uses the Green potential [START_REF] Green | Bayesian reconstructions from emission tomography data using a modified EM algorithm[END_REF] for approximating the absolute value function in the regularization:

φ Green (x) = 27 128 log cosh 16 3 √ 3 x . (6)
In Weidinger2016, as in Long2014, surrogates are derived sequentially. However, contrary to Long2014, Weidinger2016 does not assume that the material maps are non-negative, and therefore cannot find a rigorous quadratic surrogate. The quadratic function is actually not a surrogate of the convex surrogate, but its second-order Taylor approximation. This means that, if its derivative is negative, that quadratic function violates the conditions for being a surrogate. Such a choice simplifies the calculations, but the method is no longer guaranteed to converge. In our experiments, though, it did converge.

Mechlem2018

Mechlem2018 [START_REF] Mechlem | Joint Statistical Iterative Material Image Reconstruction for Spectral Computed Tomography Using a Semi-Empirical Forward Model[END_REF] builds upon Weidinger2016 using the Huber function for the regularization (Equation 4), meaning it also uses the second-order Taylor approximation of the convex surrogate as if it were an actual quadratic surrogate, which implies the same risk on convergence as in Weidinger2016. Still, it adds three features with respect to Weidinger2016: calibration, Ordered Subsets and Nesterov acceleration. The calibration part allows to estimate the product of the incident spectrum and detector response, in cases where it is either not known at all, or known with insufficient accuracy.

In this study, we have assumed perfect knowledge of these system characteristics, and therefore dropped the calibration part. In addition to that, Mechlem2018 [START_REF] Mechlem | Joint Statistical Iterative Material Image Reconstruction for Spectral Computed Tomography Using a Semi-Empirical Forward Model[END_REF] suggests to compute the Hessian (and its inverse) only once, thus updating only the gradient at each iteration. This approach is only valid if the Hessian is near-constant, i.e. if the estimated volumes undergo little changes over the course of the iterations. Therefore, the paper proposes to start the iterations from a previously known approximate result, like a set of two-step reconstructed volumes. Since we start from zero-filled estimates, we do not use this acceleration.

Barber2016

Barber2016 [START_REF] Barber R F | An algorithm for constrained one-step inversion of spectral CT data[END_REF] provides the calculations to minimize either a transmission Poisson likelihood or a least-squares data-attachment term. We chose to focus on the transmission Poisson likelihood term. As suggested by [START_REF] Sidky | Three material decomposition for spectral computed tomography enabled by block-diagonal step-preconditioning[END_REF], we have replaced the exponential function in Equation 2 with the soft exponential defined by

softexp(x) = exp(x) if x ≤ 0, x + 1 else. ( 7
)
Regularization is handled with a Total Variation (TV) term, and the cost function is minimized using the mirrored convex/concave (MOCCA) algorithm [START_REF] Barber | MOCCA: Mirrored convex/concave optimization for nonconvex composite functions[END_REF]. MOCCA is a primal-dual algorithm designed to minimize non-convex functions that can be expressed as the difference of two convex functions. As explained in section 3.1 of [START_REF] Barber R F | An algorithm for constrained one-step inversion of spectral CT data[END_REF], the algorithm should have an outer and an inner loop, but the inner loop is set to perform only a single iteration to speed up convergence. Unfortunately, that modification can cause the algorithm to diverge, so an additional parameter λ has been introduced and tuned to avoid that. Roughly, λ is tuned to the largest value that does not cause divergence. We also adapted the parameter θ of the Chambolle-Pock algorithm (see the last line of equation 7 in [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to[END_REF]), which must lie in [0; 1], and in [START_REF] Barber R F | An algorithm for constrained one-step inversion of spectral CT data[END_REF] is implicitly set to 1. With θ = 1, examining the sequence of iterates, we observed fast and ample oscillations, and artifacts constantly appearing and disappearing. These oscillations and artifacts did not significantly decrease in intensity over the course of the iterations. Setting θ = 0.5 got rid of them, yielding a more stable sequence of iterates. Computing the primal and dual step sizes Σ and T involves the element-wise absolute value of products of matrices, which is not trivial if both matrices have positive and negative values and are too large to be actually stored. The computation of the step sizes proved particularly tedious to implement, so we restricted our study to a µ-preconditioning method (see section 2.3) that does not generate negative attenuation values (i.e. simply normalizing the attenuation coefficient vectors), which makes simplifications possible. This µ-preconditioning is different from the one presented in the original paper [START_REF] Barber R F | An algorithm for constrained one-step inversion of spectral CT data[END_REF] (i.e. orthonormalizing the attenuation coefficient vectors), but our study on Cai2013 in section 3.4 shows a better convergence speed for normalization than for orthonormalization, so it is unlikely that switching to orthonormalization for Barber2016 would greatly improve its convergence speed.

Method

Cost function

D(x) + m R(x m , w m ) Parameters of Data attachment Regularization the cost function D R Cai2013 k d , w m , δ m W (y -ȳ) 2 2 + log(ȳ) with W = diag[k d ȳ] -1/2 w m φ Huber (∇x, δ m ) Long2014 w m , δ m ȳ -y log ȳ w m Nv v=1 ṽ∈N φ Long (z k-1,v -z k-1,ṽ , δ m ) Weidinger2016 w m w m Nv v=1 ṽ∈N φ Green (z k-1,v -z k-1,ṽ ) Mechlem2018 w m , δ m w m Nv v=1 ṽ∈N φ Huber (z k-1,v -z k-1,ṽ , δ m ) Barber2016 w m χ(x, w m ) = 0 if (|∇x|) 1 ≤ w m , +∞ else.

Mu-preconditioning

We form the matrix M as follows:

M =    µ I,1keV µ Gd,1keV µ H2O,1keV . . . . . . . . . µ I,150keV µ Gd,150keV µ H2O,150keV    ( 8
)
where µ me is the mass-attenuation coefficient of material m for incident X-ray photons of energy e, I stands for iodine, Gd for gadolinium and H2O for water. For e ≤ 150 keV, µ H2O,e and µ Gd,e typically differ by 2 orders of magnitude, which causes conditioning issues in the resolution of the one-step inversion problem, and slows down convergence. Section 3.4 of [START_REF] Barber R F | An algorithm for constrained one-step inversion of spectral CT data[END_REF] proposes to modify the material attenuation matrix M to speed up and stabilize the internal computations. They call this step "µ-preconditioning". It can be interpreted as creating synthetic materials, which are linear combinations of the original ones, possibly with negative coefficients. We present three different µ-preconditioning methods. Each one yields a synthetic material attenuation matrix M and a matrix P such that M P = M :

• Normalization: M is obtained by dividing each column of M by its norm. This has the advantage of not introducing negative attenuations. The added complexity associated with the use of negative attenuation coefficients, like the use of absolute values in the calculation of Σ and T in [START_REF] Barber R F | An algorithm for constrained one-step inversion of spectral CT data[END_REF], can be avoided with this synthetic materials basis.

• Orthonormalization: M is obtained either using the Gram-Schmidt algorithm on M , or diagonalizing M T M (as proposed in [START_REF] Barber R F | An algorithm for constrained one-step inversion of spectral CT data[END_REF]). This will generate negative attenuation coefficients (for the dot products between columns to be zero).

• Fessler's method:

K = SM S1
i.e. (K) b,m = e s be µ me e s be [START_REF] Ducros | Regularization of Nonlinear Decomposition of Spectral X-ray Projection Images Medical Physics[END_REF] with s be the effective spectrum at energy e in bin b, and µ me the attenuation of material m at energy e. Each value (K) b,m is a weighted mean of the attenuation coefficients of material m over the whole energy spectrum, the weights being the number of photons seen by the detector in bin b. Roughly, it evaluates how much the material m attenuates the photons that will end up being detected in bin b. Now, the problem of determining which linear combination of materials causes the observed attenuation would be much easier if each material attenuated the photons of only one bin, i.e. if K was diagonal. Therefore, we define M = M K -1 , i.e. P = K -1 . When K is not square, e.g. when there are more bins that materials, K cannot be inverted, but P can be defined as the Moore-Penrose pseudo inverse of K, i.e. P = (K T K) -1 K T . This approach was proposed in [START_REF] Fessler | Method for statistically reconstructing images from a plurality of transmission measurements having energy diversity and image reconstructor apparatus utilizing the method[END_REF] for dual energy CT, but its extension to more materials and more bins is straightforward.

When the algorithm is fed a matrix M = M P of the attenuations of synthetic materials, it reconstructs synthetic material volumes x. No change is required for the data-attachment term, but the regularization must be perfomed on the original materials volumes x = P x, since we want the original material volumes to be regular, not the synthetic ones. The changes in algorithms implied by µ-preconditioning appear in the pseudo code of Cai2013 and Barber2016 (not in the SQS methods, since they are insensitive to µ-preconditioning) in section 5.

Ordered subsets

In some methods, iterations are divided into sub-iterations, and during each sub-iteration, only a subset of the projections through the current volume is computed, compared to the measured ones, and backprojected to obtain an update for the current volume [START_REF] Hudson | Accelerated image reconstruction using ordered subsets of projection data[END_REF]. The subsets form a partition of the full projections dataset. After one iteration, i.e. one cycle through all subsets, each projection has been forward and back projected once and only once. At each iteration, OS-methods therefore perform as many updates of the reconstructed volume as the number of subsets, which results in faster convergence per iteration (or equivalently, per forward/back projection), at the cost of longer iterations, and potential instabilities if the measured data contains inconsistencies, e.g. noise or motion.

Simulation setup

We have designed a simple two-dimensional 3-materials phantom (see Figure 3), consisting of a large square of water at 1 g/ml, a small square of iodine at 10 mg/ml, and a small square of gadolinium at 10 mg/ml. The iodine and gadolinium squares are inside the water square, but do not overlap. Note that the K-edge of iodine is at 33.2 keV, which is too low to reliably use it as a K-edge contrast agent.

On the other hand, the basis water/iodine can be used to represent biological samples [START_REF] Si-Mohamed | Evaluation of spectral photon counting computed tomography K-edge imaging for determination of gold nanoparticle biodistribution: In vivo Nanoscale[END_REF], just like the Compton/photoelectric basis, with the bones appearing in the iodine material map. Therefore, we emphasize that water/iodine/gadolinium is a "two materials for biological tissues, one K-edge contrast agent" situation, and not a dual contrast agent one. The phantom consists of a single slice, which has 256 2 voxels, of size 1 mm × 1 mm. Through this phantom, 725 parallel projections were simulated, with 362 pixels per projection, using the sparse forward projection matrix generated by the AIR toolbox [START_REF] Hansen | AIR Tools -A MATLAB package of algebraic iterative reconstruction methods[END_REF].

Detector pixels too have a size of 1 mm. The line integrals obtained were then converted to photon counts, following the classical polychromatic Beer-Lambert attenuation law. For the incident spectrum, we used the same one as in Ducros et al. 2017 [START_REF] Ducros | Regularization of Nonlinear Decomposition of Spectral X-ray Projection Images Medical Physics[END_REF] which, quoting the paper, "was simulated using the SpekCalc software [START_REF] Poludniowski | SpekCalc: A program to calculate photon spectra from tungsten anode x-ray tubes[END_REF], considering a tube voltage of 120 kV, a 12 • anode angle, and a filter of 1.2 mm of Al". Figure 1 shows the number of incident photons per keV on each pixel. The detector response was simulated according to the model presented in appendix A.2 of [START_REF] Schlomka | Experimental feasibility of multi-energy photon-counting K-edge imaging in pre-clinical computed tomography[END_REF], and is shown in Figure 2. Following the design of the Philips spectral CT scanner prototype installed in Lyon, France, 5 bins were simulated. The thresholds were set to 30, 51, 62, 72 and 83 keV. In the end, the photon counts were corrupted with Poisson noise. This model therefore neglects pile-up, scatter, charge sharing and probably many other complex effects. Since the problem is nonconvex, initialization matters. It is unclear what should be the starting point of the iterations: in the original papers, Cai2013 and Barber2016 use zero-filled material slices as a starting point, while the SQS-based methods are intialized with a pre-existing reconstruction (e.g. one obtained from a two-step method). In this study, we have adopted a zero initialization for all methods.

Iodine Gadolinium Water

Figure 3: Ground truth of the phantom

Comparison method

In any iterative reconstruction method, the convergence speed, image quality at convergence and regularization parameters depend on each other. Therefore, in order to obtain meaningful comparison results on convergence speed, one must first set the regularization parameters to values that yield comparable image quality at convergence for all methods. Only then can one measure how fast a given convergence criterion is reached. To set the regularization parameters, we ran many simulations, exploring the space of parameters with multiplicative steps of 3 or 10. We then selected parameters yielding similar standard deviations in the iodine, gadolinium and water ROIs, as shown in Table 3. The regularization parameters are summarized in Table 2. Finally, the criterion we used for the convergence speeds comparison is the number of iterations each method required for all materials to reach their target concentration in their ROI, within either 20% or 10% tolerance. Besides that, we also provide:

• the concentrations of iodine and gadolinium in their respective ROIs over the course of the iterations,

• the time per iteration with our implementation. Measurements were performed on a machine equipped with an Intel Xeon E5-2620 v4 CPU running at 2.10GHz and an Nvidia GeForce GTX Titan X, and running OpenSuse Leap 42.2 and Matlab R2017b,

• the mean and standard deviation (std) in each ROI at the last iteration, which were used to adjust regularization parameters,

• the theoretical memory footprint, since the actual memory footprint depends too much on the implementation.

Since storing all iterates is impractical, for all methods except Mechlem2018 we stored one iterate every ten, which is why Table 4 contains mostly multiples of ten.

Software

For this study, we re-implemented all of the five methods compared here in a lightweight Matlab framework.

That code is open source, available at https://github.com/SimonRit/ OneStepSpectralCT, and we appreciate any relevant contribution to it. In particular, we encourage the authors of the five methods we have studied to check that code against their own implementation. 1). For Barber2016, the "regularization" parameters are actually the TV limits for iodine, gadolinium and water. The δ parameter is used in the Huber function or its equivalent, in the regularization terms of the methods.

Method

Together with the code, we provide detailed notes on how to re-do most of the calculations of the original papers.

Results

Figure 4 shows the last iterate computed for each of the methods presented here. For Cai2013, we used Fessler's method for µ-preconditioning. A video of the iterates, on which it is easy to evaluate the relative convergence rates of the various methods, is available as supplementary material. It has the same grayscale and is arranged the same way as Figure 4.

Iodine Gadolinium Water

Algorithm parameters

In addition to the quantitative criteria on convergence speed and memory footprint, a few other important aspects of the algorithms deserve to be mentioned here. All methods require parameters for regularization and Cai2013 needs a parameter k d characterizing the detector for its data attachment term (see Table 1). The rest of the parameters are free parameters for optimizing the cost function.

Our implementation of Cai2013 requires a threshold on how many times one tries to halve the step size before resetting the descent direction (see section 2.2.1). Barber2016 has a parameter λ = 10 -4 , and a parameter θ = 0.5 (see section 2.2.5). Long2014 splits the projections into 20 subsets, and Mechlem2018 into 4 subsets, in order to speed up convergence. Mechlem2018 also uses Nesterov's momentum technique, which speeds up convergence even more, but makes it a bit unstable: when splitting the projections into too many subsets (typically more than 6), Mechlem2018 sometimes diverges: the iterates accumulate strong streak artefacts, then they reach a point where the Hessian of the SQS is not invertible, which stops the algorithm.

Quantitative criteria

We show in Figure 5 the evolution of the concentrations of iodine and gadolinium over the course of iterations, for all five methods, with a logarithmic scale only on the x-axis. The concentrations are measured in the square where iodine or gadolinium is expected, eroded by 2 pixels in order to avoid partial volume effects. In all methods, the pixels at the border of the region of interest (ROI) have significantly lower concentrations. The full line is the ground truth concentration value in the ROI.

Since tracking the mean concentration in a ROI can only capture information on how well the low frequencies are reconstructed, Figure 6 shows, for each method, the normalized 2 difference with the last iterate, computed as follows:

Normalized 2(x k ) = Nm m Nv v=1 (x k,v,m -x N iterations ,v,m ) 2 N m Nv v=1 GroundTruth 2 v,m (10)
where k, v and m are respectively the indices for iteration, voxel and material, N v is the number of voxels of one reconstructed material map and N m the number of materials. Table 3 contains the mean and standard deviation of iodine, gadolinium and water concentrations in the ROI where they are expected eroded by 2 pixels, on the final iterate. Standard deviations of each material are roughly the same between all algorithms, yet there are some outliers: the stds of Cai2013 in gadolinium and of Barber2016 in iodine. This is a consequence of the fact that those two methods are the slowest ones, and consequently we were not able to sample the parameter space for them with a fine enough grid. While it is unfortunate, it highlights the fact that slowly converging methods are hard to tune, and that additional parameters (λ in Barber2016) imply a longer tuning time.

Table 4 shows the total number of iterations performed, the time per iteration, and the number of iterations it took for all materials to reach target concentration in their respective ROIs, within 20% or 10% tolerance. On our test case, Mechlem2018 converges 1 to 3 orders of magnitude faster than its competitors.

Memory footprint of each method

We provide here a simple analysis of the memory footprint of each method, by listing the variables that have to be stored and recalling their size. The results are summarized in Table 5, and the detailed calculations for each method are presented in this section. The measured photon counts, which constitute the main input of these reconstruction algorithms, are assumed to be available on disk, so they never have to be fully loaded in memory and are not mentioned in this section. Obviously, re-reading them from disk at each iteration is highly inefficient, and one should load them in memory if possible. We used the notations of the pseudo codes in the Appendix for the names of the variables, and the following notations for various quantities: N v is the number of voxels of one reconstructed material map, N m the number of materials, N p the number of pixels of all projections for one bin, N b the number of energy bins in photon counts measurements and N e the number of energies used in internal calculations. Despite the fact that we performed 2D simulations, we assume for these calculations that the reconstructed volumes are 3D (which matters for the computation of the spatial gradient). Note that the code we provide keeps many intermediate variables in memory, for convenience and speed considerations, and therefore uses much more memory than the minimum required.

Method

Minimal memory footprint (in number of floating point values)

Cai2013 4N v N b with Fessler's preconditioning 4N v N m otherwise Long2014 4 + Nm+1 2 N v N m Weidinger2016 Mechlem2018 6 + Nm+1 2 N v N m Barber2016 17N v N m + 3N p N b
Table 5: Summary of the minimal memory footprints of the five studied methods.

Cai2013

The following variables must be stored and updated at each iteration:

x k , d k , g k and g k-1
, which each contain a set of material volumes. Since all computations can be performed pixel-wise, the expected photon counts ȳ do not have to be fully stored in memory, and neither does any of the intermediate variables in the calculation of the gradient and Hessian of the data-attachement term, so we do not count them. Lastly, the spatial gradient of x k , which is 3 times as large as x k , is used in the calculation of g and of d T Hd. More precisely, it is the divergence of the first and second derivative of the Huber function of the gradient that is used. These calculations can be computed piecewise, with overlapping pieces since the gradient and divergence require neighboring voxels. In total, at least 4N v N m floating point values must be kept in memory, and storing 3N v N m additional values for the gradient of x k avoids redundant computations, and the problem of having to divide the volumes into overlapping pieces.

Long2014, Weidinger2016 and Mechlem2018

In all three SQS-based methods we have implemented, the following variables must be stored and updated at each iteration: the current iterate (x k as in [START_REF] Long | Multi-Material Decomposition Using Statistical Image Reconstruction for Spectral CT[END_REF]), of size N v N m , and the gradient and Hessian of the surrogate at x k , of size N v N m and N v N 2 m , respectively. Since the Hessian is known to be symmetric, it can be stored in N v N m (N m + 1)/2 floating point values. The gradient and the Hessian of the surrogate of the data-attachment term both result from one forward and one back projection. By performing these one angle at a time, they can be computed with negligible amounts of memory. The gradient and Hessian of the surrogate of the regularization term both require a memory size of N v N m (the Hessian has N v N 2 m elements, but is known to be diagonal), and internal computations can be done in place. Since typically they are calculated separately, and then added to those of the regularization terms, their requirements in memory add up. In total, the memory footprint is at least 4 + Nm+1

Barber2016

This method makes use of more intermediate variables than the other four. Table 6 lists the variables used in the pseudo code given in Appendix, whether they must be stored or not, why, and how many floating point values are therefore required. In total, it uses at least 17N v N m + 3N p N b floating point values.

Variable

Stored ? Size in memory (in floating point values)

x 0 yes N v N m x (k-1) x (k) x (k+1) only x (k-1) N v N m Σ sino yes N p N b Σ grad w only one of both 3N v N m
T yes: used only once, but computed as the sum of two very different calculations, unlikely to occur voxelwise in parallel

N v N m z no: used only once, pixelwise 0 u (k-1) sino u (k) sino u (k+1) sino not u (k+1)
sino , because it replaces u

x (k) x (k+1) yes 2N v N m
Table 6: Memory footprint of Barber2016, for each variable used in appendix C.6 of [START_REF] Barber R F | An algorithm for constrained one-step inversion of spectral CT data[END_REF]. N v , N m , N p and N b are defined in section 3.3 of the present paper. The words "pixelwise" and "voxelwise" mean that the variable can be computed one pixel (or voxel) at a time, thus its computation requires negligible amounts of memory.

Impact of mu-preconditioning

None of the SQS-based methods is impacted by µ-preconditioning: once mapped back to the original materials, the iterates are almost identical to the ones obtained without µ-preconditioning, with difference maps two orders of magnitude less intense than the reconstructed slices. The only impact is with Fessler's method [START_REF] Fessler | Method for statistically reconstructing images from a plurality of transmission measurements having energy diversity and image reconstructor apparatus utilizing the method[END_REF]: since in our case (three materials, five energy bins) it artificially increases the number of materials to that of the number of bins, the Hessians of the SQS are no longer invertible, so this µ-preconditioning cannot be used at all with SQS-based methods. On Cai2013, however, µpreconditioning has a strong impact. 7 shows the evolution of the cost function over the course of the iterations minus its minimum, using Cai2013 without µ-preconditioning, with normalization, orthonormalization and Fessler's method. Since all four of these experiments aim to minimize the same cost function, the minimum we subtracted is the minimum over all iterations and all four experiments (reached, as can be seen on the graph, with µ-preconditioning by normalization). It clearly shows that Fessler's method has the fastest initial convergence, which is confirmed by looking at the iterates (not shown here). After more than 1000 iterations, the experiment using normalization reached a lower value of the cost function than the one using Fessler's method.

Discussion

We have presented a comparison of five one-step inversion methods for spectral CT, focused on two aspects important for implementation into a real scanner: convergence speed and memory footprint. This choice of criteria is in part due to the initial aim of our study, which was to determine which algorithm should be implemented to process large real datasets in reasonable time. But more importantly, it was dictated by the behavior of the algorithms themselves: while it seemed impossible to make a convincing case on which method performed best in terms of image quality, the differences in convergence speed were striking. We provide the complete Matlab code we have used for this study (see section 2.7), so that readers interested in conducting comparisons on other aspects have the means to do so.

All algorithms were adapted to find the three material maps that minimized a cost function composed of a data attachment term for five energy bins and a spatial regularization minimizing the 1-norm (or a twice-differentiable approximation of it) of the spatial gradient of each reconstructed volume. The original algorithms were generally designed for other problems and our conclusions hold only for our specific problem. For examples, the following modifications of the problem would require another comparison: reconstructing two maps from dual-energy measurements as in [START_REF] Cai | A full-spectral Bayesian reconstruction approach based on the material decomposition model applied in dual-energy computed tomography[END_REF] or constraining the solutions with positivity as in [START_REF] Long | Multi-Material Decomposition Using Statistical Image Reconstruction for Spectral CT[END_REF].

There are other one-step inversion methods in the literature than the ones we selected, e.g. [START_REF] Chen | Image reconstruction and scan configurations enabled by optimization-based algorithms in multispectral[END_REF][START_REF] Sukovic | Penalized weighted least-squares image reconstruction for dual energy X-ray transmission tomography[END_REF], and it is likely that more will be proposed in the future. We hope that others will be interested in implementing their method in the Matlab framework we provide, to extend the comparison.

A limitation of the present study is the small size of the data used in the experiments, which was dictated by our implementation using a sparse system matrix stored in computer memory. It is therefore possible that on larger datasets, e.g. on a full real helical acquisition, the relative convergence speeds of the studied methods would be different. However, the other observations made in the paper would remain valid.

Another limitation is that since one-step inversion is a non-convex problem, the results may strongly depend on the initialization. Mechlem et al. recommend to initialize the algorithms with reconstruction results from a two-step method. We have not tried that approach, and have instead initialized all reconstructions with zero-filled slices. Different initialization strategies may lead to different conclusions on relative convergence speeds.

In several ways, the simulation we have studied is simplistic: in a real acquisition, the spectrum is usually filtered to remove most low-energy photons, made pixel-dependent by a bow tie filter, and the pixels and voxels are much finer than the ones we used (around 0.2 mm instead of 1 mm). Since the same sparse system matrix was used both for simulation and reconstruction, this work adopts an inverse-crime setup [START_REF] Wirgin | The inverse crime[END_REF], with statistical noise. It does not take into account any advanced physical effect like scatter, pile-up or charge sharing. However, these observations hold for all methods, so the comparison remains valid.

Sidky et al. have recently proposed a more efficient µ-preconditioning method [START_REF] Sidky | Three material decomposition for spectral computed tomography enabled by block-diagonal step-preconditioning[END_REF], better-suited to three-materials reconstructions than the one described in Barber2016. In the experiment they present, this new µ-preconditioning method speeds up convergence by a factor 10. The slow convergence of Barber2016 in our experiments must therefore be mitigated by the fact that we did not implement that method. However, a factor 10 would only make its convergence speed similar to Cai2013, Long2014 and Weidinger2016, and our other observations on Barber2016 (mostly its high memory footprint) remain valid.

Our implementation of Mechlem2018 differs from that of Weidinger2016 only by the use of Nesterov's momentum technique and a small number of ordered subsets, and yet greatly outperforms Weidinger2016. Therefore, we confirm the observation made in the discussion of Mechlem2018, and in several other papers [START_REF] Kim | Combining Ordered Subsets and Momentum for Accelerated X-Ray CT Image Reconstruction[END_REF][START_REF] Hansen | Fast 4D cone-beam CT from 60 s acquisitions[END_REF]: Nesterov's momentum technique combined with ordered subsets highly accelerates the convergence speed per iteration. This acceleration might cause numerical instabilities which would make Mechlem2018 more difficult to use with real data.

Conclusion

We have adapted five one-step methods for spectral CT (Cai2013 [START_REF] Cai | A full-spectral Bayesian reconstruction approach based on the material decomposition model applied in dual-energy computed tomography[END_REF], Long2014 [START_REF] Long | Multi-Material Decomposition Using Statistical Image Reconstruction for Spectral CT[END_REF], Weidinger2016 [START_REF] Weidinger | Polychromatic Iterative Statistical Material Image Reconstruction for Photon-Counting Computed[END_REF], Barber2016 [START_REF] Barber R F | An algorithm for constrained one-step inversion of spectral CT data[END_REF], Mechlem2018 [START_REF] Mechlem | Joint Statistical Iterative Material Image Reconstruction for Spectral Computed Tomography Using a Semi-Empirical Forward Model[END_REF]) to reconstruct three material maps from five energy bins with spatial regularization in each material map. We have compared them in terms of practical aspects: convergence speed to reach a fixed image quality threshold, and memory footprint. They are obviously not the only criteria for the choice of a method, and many other parameters should probably be considered such as image quality at convergence or numerical stability. Yet, Mechlem2018 solved this one-step problem in 1 to 3 orders of magnitude faster than its competitors, and such a vast difference should not be overlooked. Alternative optimizations or modifications of the problem would require reinvestigating this conclusion, and we provide the source code of the software used in this study to facilitate such investigations.

5.1), already contains the regularization weights (one per material) and, if applicable, the thresholds (one per material). For SQS-based methods, the regularization is computed by examining, for each voxel, the surrounding voxels that belong to a neighborhood N . For Long2014 and Mechlem2018, which use subsets, we define A s = sub(A, s) which selects rows in the projection matrix A corresponding to the projections in subset s. Which projection goes into which subset is determined by the following procedure: we compute a random permutation of [1, 2, . . . , N projections ], and divide the result into N subsets equal parts.

Pseudo-code of Cai2013

Contrary to the other four methods, Cai2013 defines S as the normalized spectrum instead of the absolute one, and y as transmission ratios instead of photon counts. Let us define w = (w 1 , . . . , w Nm ) and δ the vectors of regularization weights, and of Huber function thresholds, each holding one value per material. In Cai2013, the regularization term is the sum over all voxels of the Huber function of the spatial gradient of the material volumes, i.e.

Φ :R 3NvNm × R Nm × R Nm → R (t, w, δ) → Nm m=1 w m 3Nv i=1 φ(t i , δ m ) (12)
where in practice, t is always the spatial gradient of the material volumes. The gradient and Hessian of Φ are denoted Φ and Φ , respectively, and are of size

[3N v N m , 1] and [3N v N m , 3N v N m ].

Pseudo-code of Long2014

It uses the optimal curvature function, which is defined on R + as follows:

optimal curvature(x) = 2 1 -e -x -xe -x x 2 if x > 0, 1 if x = 0 (13)
When the input to optimal curvature is a vector, it is applied element-wise. All other four methods use a weighted regularization term, with weights controlling the strength of the regularization: the larger the weight, the stronger the regularization. Barber2016 works differently: it enforces a constraint, stating that the TV on each material must not exceed a threshold γ m (one per material). Therefore, the larger the γ m s, the looser the TV constraint. x(k+1) ← x (k+1) + θ(x (k+1) -x (k) )

27 end 28 return P x (k) Algorithm 5: Barber2016 pseudo code

Figure 1 :

 1 Figure 1: Incident spectrum on each pixel.

Figure 2 :

 2 Figure 2: Spectral response of each detector bin.

Figure 4 :

 4 Figure 4: Reconstructed slices with the five methods, with grayscale ranges equal to [0; 15] mg/ml for iodine and gadolinium, and [0; 1.5] g/ml for water.

Figure 5 :

 5 Figure 5: Concentration of iodine (left), gadolinium (middle) and water (right) over the course of the iterations for the five methods, in the ROI where it is expected eroded by 2 pixels.

Figure 6 :

 6 Figure 6: Normalized 2 difference with the last iterate.

2 N

 2 v N m floating point values for Long2014 and Weidinger2016. For Mechlem2018, two additional Nesterov variables must be stored (z k and u k ), taking the memory footprint to a minimum of 6 + Nm+1 2 N v N m floating point values.

Figure 7 :

 7 Figure 7: Cost function of Cai2013 over the course of the iterations minus its mininum, depending on the µ-preconditioning method used.

1 P 11 z ← y 2 -ȳ 2 k d ȳ2 - 1 ȳ 12 x 13 g 2 15 18 d real ← P d k 19 d 27 Q 1 : 8 for s ← 1 to N subsets do 9 A 11 Q 12 ȳ ← SQ 13 g 18 H 19 for v ← 1 to N v do 20 g regul,v ← 2 N 21 H 22 g v ← g data,v + g regul 23 Hv g v 25 end 26 end 27 end 28 return x k Algorithm 2 : 7 Q 8 ȳ ← SQ 9 g 10 H 11 for v ← 1 to N v do 12 g 13 Hv g v 17 end 18 end 19 return x k Algorithm 3 : 7 for s ← 1 to N subsets do 8 A 9 Q 10 ȳ ← SQ 11 g 12 H 13 for v ← 1 to N v do 14 g 15 Hv g v 19 end 20 t k ← 1 2 1 + 1 + 4t 2 k- 1 21α k ← z k -d k 22 u k ← u k -t k-1 d k 23 z k ← α k + t k k l=0 t l (u k -α k ); 24 end 25 end 26 return z k Algorithm 4 :

 11121121321819271891112131819202212223252627278910111213378910111213141519201121222324254 = P ⊗ I Nv 2 A = I Nm ⊗ A 3 M = M ⊗ I Np 4 S = S ⊗ I Np 5 x 0 ← P -1 X 6 g 0 ← 0 7 d 0 ← 0 8 Q ← exp (-M Ax 0 ) 9 ȳ ← SQ 10 for k ← 1 to N iterations do real ← P x k-1 k ← A T M T QS T z + P T ∇ T Φ ∇x real , W , δ 14 β ← <g k -g k-1 ,g k-1 > g k-1 2 if β < 0 or k = 1 then β ← 0 16 d k ← -g k + βd k 17 v ← diag[ 2y 2 k d ȳ3 -1 ȳ2 ] T k Hd k ← d T k A T M T Q(S T vSQ -diag[S T z])M Ad k 20 +d T real P T ∇ T Φ ∇x real , W , δ ∇P d real 21 α ← -<g k ,d k > d T k Hd k 22 x k ← x k-1 + αd k 23 while cost(x k ) >cost(x k-1 ) do 24 α ← α/2 25 x k ← x k-1 + αd k ; 26 end ← exp (-M Ax k ) 28 ȳ ← SQ 29 end 30 return P x k Algorithm Cai2013 pseudo code 1 π = A * (A1 Nv,1 ) 2 A = I Nm ⊗ A 3 L = (M ⊗ 1 1,Nm ) * (1 1,Nm ⊗ M ) 4 M = M ⊗ I Np 5 S = S ⊗ I Np 6 x 0 ← X 7 for k ← 1 to N iterations do s = sub(A, s) 10 π s = sub(π, s) ← exp (-M A s x k-1 ) data ← A T s ( y ȳ -1) * SQM 14 c = optimal curvature(M Rx k-1 )15 for p ← 1 to Npixels in subset s do16 C p = N b b=1 (S p (c p * L)) 17 end data ← reshape(π T s C, [N v , N m , N m ]) subsets ṽ∈N φ (x k-1,v -x k-1,ṽ ) regul,v ← 4 diag ṽ∈N φ (x k-1,v -x k-1,ṽ ) v ← H data,v + H regul 24 x k,v ← x k-1,v -H -1 Long2014 pseudo code 5.3. Pseudo-code of Weidinger2016 1 A = I Nm ⊗ A 2 L = (M ⊗ 1 1,Nm ) * (1 1,Nm ⊗ M ) 3 M = M ⊗ I Np 4 S = S ⊗ I Np 5 x 0 ← X 6 for k ← 1 to N iterations do ← exp (-M Ax k-1 ) data ← A T ( y ȳ -1) * (SQM ) data ← A T ((A1 Nv,1 ) * (SQL)) regul,v ← 2 ṽ∈N φ (x k-1,v -x k-1,ṽ ) regul,v ← 4 diag ṽ∈N φ (x k-1,v -x k-1,ṽ ) 14 g v ← g data,v + g regul 15 H v ← H data,v + H regul 16 x k,v ← x k-1,v -H -1 Weidinger2016 pseudo code5.4. Pseudo-code of Mechlem2018It differs of that of Weidinger2016 only by the use of subsets and Nesterov's momentum technique.1 A = I Nm ⊗ A 2 L = (M ⊗ 1 1,Nm ) * (1 1,Nm ⊗ M ) 3 M = M ⊗ I Np 4 S = S ⊗ I Np 5 z 0 ← X 6 for k ← 1 to N iterations do s = sub(A, s) ← exp (-M A s z k-1 ) data ← A T s ( y ȳ -1) * (SQM ) data ← A T s ((A s 1 Nv,1 ) * (SQL)) regul,v ← 2 N subsets ṽ∈N φ (z k-1,v -z k-1,ṽ ) regul,v ← 4 diag ṽ∈N φ (z k-1,v -z k-1,ṽ ) 16 g v ← g data,v + g regul 17 H v ← H data,v + H regul 18 d k,v ← z k-1,v -H -1 Mechlem2018 pseudo code 5.5. Pseudo-code of Barber2016

1 Z 8 x 0 ← x(k) 9 Bpp 10 with I pp = 1 if p = p and 0 otherwise 11 K ← B(x 0 )Z 12 ȳ 13 E 15 Σ 16 Σ grad ← 1 λ 17 T 18 ∀m, w m ← Σ 1 25 x

 18910111213151611718125 pe,mv ← M me A pv 2 A = I Nm ⊗ A 3 M = M ⊗ I Np 4 S = S ⊗ I Np 5 P = P ⊗ I Nv 6 x (0) ← P -1 X 7 for k ← 1 to N iterations do bp,p e (x) ← s bpe exp(-(Zx)pe) e s bpe exp(-(Zx) pe ) I ← S exp (-M Ax 0 ) ← max(ȳ -y, 0) 14 b ← (ȳ -E)Kx 0 + ȳ -y sino ← 1 λ diag [K1] -1 diag [|∇|P 1] -1 ← λ diag K T 1 (|∇|P 1) T -1 sino ← (ȳ + Σ sino ) -1 ȳ(u (k) sino + Σ sino K x(k) ) -Σ sino (b + Ez) 21 ∀m, u + grad,m ← u (n) grad,m + Σ grad,m ∇(P x(k) ) m 22 ∀m, g + m ← |u + grad,m /w m | 23 ∀m, ĝ+ m ← sign(u + grad,m /w m ) 24 ∀m, u(k+1)grad,m ← u + grad,m -w m ĝ+ m Proj(g + m ; g, g/w m ≤ γ m ) (k+1) ← x(k) -

Table 1 :

 1 Summary of the cost functions minimized by the studied methods. φ

Huber

, φ Long and φ Green are twicedifferentiable approximations of the absolute value function, defined in Equation 4, Equation 5 and Equation

6

, respectively. ∇ is the 3-components spatial gradient, and the neighborhood N is arbitrary, but always set to all 26 immediate neighbors of the voxel v. The regularization of Barber2016 uses the convex indicator function χ to enforce an upper bound on the total variation (|∇x|) 1 of x. The numerical values of the parameters of the cost function used for the simulated experiment are listed in Table

2

.

Table 2 :

 2 δ iodine δ gadolinium δ water w iodine w gadolinium w water Regularization parameters of the five methods used to obtain the reconstructions (Table

	Cai2013	0.001	0.001	0.1	100000	100000	30
	Long2014	0.001	0.001	0.1	100000	100000	10
	Weidinger2016 N/A	N/A	N/A 30000	30000	3
	Barber2016	N/A	N/A	N/A	100	100	5000
	Mechlem2018 0.001	0.001	0.1	30000	30000	3

Table 3 :

 3 Mean and standard deviation of iodine, gadolinium and water concentrations in the ROI where they are expected eroded by 2 pixels on the final iterate.

	Method	Iterations performed	Time per iter (s)	Within 20% tolerance	Within 10% tolerance
	Cai2013	5000	5.2	270	430
	Long2014	5000	1.33	140	280
	Weidinger2016	5000	0.48	190	390
	Barber2016	20000	1.02	3850	8300
	Mechlem2018	200	0.62	5	10

Table 4 :

 4 Number of iterations performed and required to reach target concentrations within 20% and 10% tolerance.

Acknowledgements

This work was performed within the framework of the EU's H2020 research and innovation program under the grant agreement No. 668142, the SIRIC LYriCAN Grant INCa-INSERM-DGOS-12563, and the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the ANR. We gratefully acknowledge the support of NVIDIA Corporation with the donation of the Titan X GPU used for this research.

Appendix

Since we performed some modifications on almost all methods, we provide in this section the complete pseudo code of our implementations. This pseudo code makes extensive use of matrix notation, so we define some matrices and their size, in format [lines, columns] (see section 3.3 for the definitions of the N values): M is the (preconditioned) material mass-attenuation matrix, and P the µ-preconditioning matrix, defined in section 2. In some algorithms, those matrices are redefined using Kronecker products (noted ⊗). They still perform the same operation, only extended to several materials, bins, pixels or voxels. When the original matrix is needed, it is indexed by the pixel or voxel. Some calculations cannot be easily translated into Kronecker products, therefore we sometimes use the symbol * , which * stands for element-wise product with implicit extension, e.g.

Of course, it would have been highly inefficient to actually implement the presented methods using Kronecker products, so the pseudo code does not reflect the way things are actually implemented in the Matlab code we provide. In vectors storing volumes, the order of the dimensions is "x, y, z, material". In vectors storing photon counts, it is "x, y, projection, energy bin". If µ-preconditioning is used, the initial guess, which we denote X, must be multiplied by P -1 to obtain x 0 . In all methods except Barber2016, regularization involves the use of a twice-differentiable approximation of the 1 norm. Cai2013 and Mechlem2018 use the Huber function, Weidinger2016 uses the Green prior, and Long2014 a custom hyperbola. We denote all of these functions with a unique symbol φ which, as for Cai2013 (see section