N Ducros

J F P J Abascal

B Sixou

S Rit

F Peyrin

Regularization of Nonlinear Decomposition of Spectral X-ray Projection Images

Purpose: Exploiting the X-ray measurements obtained in different energy bins, spectral computed tomography (CT) has the ability to recover the 3D description of a patient in a material basis. This may be achieved solving two subproblems, namely the material decomposition and the tomographic reconstruction problems. In this work, we address the material decomposition of spectral x-ray projection images, which is a non-linear ill-posed problem. Methods: Our main contribution is to introduce a material-dependent spatial regularization in the projection domain. The decomposition problem is solved iteratively using a Gauss-Newton algorithm that can benefit from fast linear solvers. A Matlab implementation is available online. The proposed regularized weighted least squares Gauss-Newton algorithm (RWLS-GN) is validated on numerical simulations of a thorax phantom made of up to five materials (soft tissue, bone, lung, adipose tissue and gadolinium), which is scanned with a 120 kV source and imaged by a 4-bin photon counting detector. To evaluate the method performance of our algorithm, different scenarios are created by varying the number of incident photons, the concentration of the marker and the configuration of the phantom. The RWLS-GN method is compared to the reference maximum likelihood Nelder-Mead algorithm (ML-NM). The convergence of the proposed method and its dependence on the regularization parameter are also studied. Results: We show that material decomposition is feasible with the proposed method and that it converges in few iterations. Material decomposition with ML-NM was very sensitive to noise, leading to decomposed images highly affected by noise and artifacts even for the best case scenario. The proposed method was less sensitive to noise and improved contrast-to-noise ratio of the gadolinium image. Results were superior to those provided by ML-NM in terms of image quality and decomposition was 70 times faster. For the assessed experiments, material decomposition was possible with the proposed method when the number of incident photons was equal or larger than 10 5 and when the marker concentration was equal or larger than 0.03 g.cm -3 . Conclusions: The proposed method efficiently solves the nonlinear decomposition problem for spectral CT, which opens up new possibilities such as material-specific regularization in the projection domain and a parallelization framework, in which projections are solved in parallel.

I. INTRODUCTION

Spectral computed tomography (CT) is gaining increasing attention. Exploiting X-ray measurements acquired at multiple energies, this new imaging modality has the ability to recover the concentration maps of the constituents of the tissues in a quantitative manner. While the principle of dual-energy CT has been known for more than 30 years, recent developments in energy selective detectors have boosted the research in this area [START_REF] Polad | Energy-resolved computed tomography: first experimental results[END_REF][START_REF] Overdick | Status of direct conversion detectors for medical imaging with x-rays[END_REF][START_REF] Brunner | First k-edge imaging with a micro-CT based on the xpad3 hybrid pixel detector[END_REF][START_REF] Ballabriga | Medipix3: A 64k pixel detector readout chip working in single photon counting mode with improved spectrometric performance[END_REF][START_REF] Gimenez | Characterization of medipix3 with synchrotron radiation[END_REF][START_REF] Cheng Xu | Evaluation of a second-generation ultra-fast energy-resolved asic for photon-counting spectral CT[END_REF][START_REF] Wang | Mi-croCT with energy-resolved photon-counting detectors[END_REF][START_REF] Wang | Material separation in x-ray CT with energy resolved photon-counting detectors[END_REF][START_REF] Taguchi | Vision 20/20: Single photon counting x-ray detectors in medical imaging[END_REF]. In particular, it has been shown that photon counting detectors can be used to image high Z contrast agents by exploiting the K-edge energy discontinuity of their linear attenuation coefficient (LAC) [START_REF] Roessl | K-edge imaging in x-ray computed tomography using multi-bin photon counting detectors[END_REF][START_REF] Schlomka | Experimental feasibility of multi-energy photon-counting k-edge imaging in pre-clinical computed tomography[END_REF][START_REF] Schirra | Statistical reconstruction of material decomposed data in spectral CT[END_REF]. Therefore, coupling spectral CT acquisition to the injection of contrast agent is promising not only to increase the sensitivity of standard CT but also to open new clinical applications for X-ray imaging as a functional imaging tool, e.g., for the characterization of the atherosclerotic plaque [START_REF] Feuerlein | Multienergy photon-counting k-edge imaging: Potential for improved luminal depiction in vascular imaging1[END_REF][START_REF] Cormode | Atherosclerotic plaque composition: Analysis with multicolor CT and targeted gold nanoparticles1[END_REF][START_REF] Pan | Computed tomography in color: Nanok-enhanced spectral ct molecular imaging[END_REF].

As in standard CT, spectral CT reconstruction requires to solve an inverse problem. However, spectral CT reconstruction can exploit measurements in different * nicolas.ducros@creastis.insa-lyon.fr energy channels, which allows to recover the concentration of different materials [START_REF] Alvarez | Energy-selective reconstructions in x-ray computerised tomography[END_REF]. A wide variety of approaches for material reconstruction have been investigated. Spectral CT material reconstruction can be divided into image-based, direct inversion, and projection-based approaches. The image-based approach is a pragmatic approach that consists in the tomographic reconstruction of each of the energy sinogram followed by a decomposition step. Although it provides valuable results for multi-material decomposition [START_REF] Huy | Least squares parameter estimation methods for material decomposition with energy discriminating detectors[END_REF][START_REF] Clark | In vivo characterization of tumor vasculature using iodine and gold nanoparticles and dual energy micro-CT[END_REF][START_REF] Faby | Performance of today's dual energy CT and future multi energy CT in virtual non-contrast imaging and in iodine quantification: A simulation study[END_REF], it makes substantial assumptions concerning the forward model that may results in image artifacts and quantification inaccuracy. Development of sophisticated correction schemes to alleviate this problem is a topic of ongoing research [START_REF] Persson | Energy-resolved CT imaging with a photon-counting silicon-strip detector[END_REF][START_REF] Mendonca | A flexible method for multi-material decomposition of dual-energy CT images[END_REF][START_REF] Touch | A neural network-based method for spectral distortion correction in photon counting x-ray CT[END_REF]. Direct inversion and projection-based approaches both rely on a mathematical description of the physics of image formation. Direct inversion methods recover the material density directly from the measured energyresolved sinograms while the projection-based methods first decompose the energy-resolved sinograms into material-specific sinograms and then perform a tomographic reconstruction of each of the decomposed sinograms. Direct inversion is probably the most elegant approach [START_REF] Long | Multi-material decomposition using statistical image reconstruction for spectral CT[END_REF][START_REF] Zhang | Model-based iterative recon-struction for dual-energy x-ray CT using a joint quadratic likelihood model[END_REF][START_REF] Zhao | An extended algebraic reconstruction technique (E-ART) for dual spectral CT[END_REF][START_REF] Barber | An algorithm for constrained onestep inversion of spectral CT data[END_REF]. However, the resulting inverse problem is computationally demanding since data from all views have to be processed at the same time. On the contrary, the projection-based approach allows independent decomposition of each projection view, which can be parallelized and is computationally more tractable.

The decomposition step has led to various approaches, some rely on parametric models while the others are based on a physical model. On the one hand, a parametric model avoids measuring the source spectrum and the energy responses of the detector, which may be difficult [START_REF] Sidky | A robust method of x-ray source spectrum estimation from transmission measurements: Demonstrated on computer simulated, scatter-free transmission data[END_REF]. Despite its excellent decomposition quality [START_REF] Alvarez | Estimator for photon counting energy selective x-ray imaging with multibin pulse height analysis[END_REF][START_REF] Maaß | Empirical multiple energy calibration (emec) for materialselective CT[END_REF][START_REF] Kevin | Experimental comparison of empirical material decomposition methods for spectral CT[END_REF][START_REF] Alvarez | Efficient, non-iterative estimator for imaging contrast agents with spectral x-ray detectors[END_REF][START_REF] Wu | A weighted polynomial based material decomposition method for spectral x-ray CT imaging[END_REF], this approach is more difficult to implement when three materials or more are to be recovered. This can be a serious drawback for K-edge imaging applications. On the other hand, approaches relying on a physical model of the imaging setup can be more easily applied to arbitrary numbers of materials [START_REF] Schlomka | Experimental feasibility of multi-energy photon-counting k-edge imaging in pre-clinical computed tomography[END_REF].

In this paper, we present a new framework for the decomposition of a single projection, assuming a physical image formation model is known. The main contribution of this work is to introduce material-specific regularization in the projection domain.

Regularization is not a new idea in spectral CT which is an ill-posed inverse problem. Different regularization strategies have been investigated in order to improve the signal-to-noise ratio of the reconstructions, e.g., statistical penalties and spatial priors. Direct inversion approaches have the ability to incorporate spatial prior in the image domain. Previous works include non-linear edge-preserving statistical reconstruction [START_REF] Long | Multi-material decomposition using statistical image reconstruction for spectral CT[END_REF], Gaussian Markov random field regularization with positivity constrains [START_REF] Zhang | Model-based iterative recon-struction for dual-energy x-ray CT using a joint quadratic likelihood model[END_REF], extension of the ART algorithm [START_REF] Zhao | An extended algebraic reconstruction technique (E-ART) for dual spectral CT[END_REF], and TV-regularized primal-dual algorithm [START_REF] Barber | An algorithm for constrained onestep inversion of spectral CT data[END_REF]. However, it may be noted that these approaches may be computationally expensive in terms of memory requirements and computation time, and have mainly been limited to the reconstruction of 2D slices exploiting a low number of energy bins. Moreover, simplified forward models are often used to moderate the computational cost of the reconstruction. Similarly, projection-based approaches focus on penalties in the image domain which are imposed regularizing the tomographic problem. Edge preserving regularization based on the Huber potential function has been considered in [START_REF] Schirra | Statistical reconstruction of material decomposed data in spectral CT[END_REF][START_REF] Sawatzky | Proximal admm for multi-channel image reconstruction in spectral x-ray CT[END_REF] while total variation/ 1 sparsity constrains were used in [START_REF] Xu | Sparsity-regularized image reconstruction of decomposed k-edge data in spectral CT[END_REF]. However, to the best of our knowledge, regularization of the decomposition step, which is nonlinear and perhaps more difficult to solve than the tomographic step, was only recently addressed [START_REF] Brendel | Penalized likelihood decomposition for dual layer spectral CT[END_REF][START_REF] Ducros | Non-linear regularized decomposition of spectral x-ray projection images[END_REF].

As proposed by Schlomka et al. [START_REF] Schlomka | Experimental feasibility of multi-energy photon-counting k-edge imaging in pre-clinical computed tomography[END_REF], our forward model is nonlinear and takes into account the spectral response of the detector. In this paper, we also regularize each material projection image. Instead of considering the processing of a single projection line, we consider the decomposition of a 2D projection which is of interest for multi-line and 2D energy resolved detectors. In order to take into account the noise statistics, we have used a weighted least square data term. The inverse

R(R, Ω) ρ Material mass density R(Ω) τ Material basis function R(R) Ω Object volume R 3 u Detector pixel location R 2 x Object voxel location R 3 a
Projected mass density image

R P M s Photon counts image R IP H Hessian of R R M P ×M P J Jacobian of F R IP ×M P E Output energy R R Regularizer R(R P M ) S Detection surface R 2 F Forward model R IP (R P M )
problem is solved using a Gauss-Newton algorithm that exploits second-order derivatives, which converges faster than the most commonly used first-or zeroth-order algorithms. Gauss-Newton minimization was successfully applied to several medical imaging inverse problems [START_REF] Vauhkonen | Tikhonov regularization and prior information in electrical impedance tomography[END_REF][START_REF] Schweiger | Gauss-Newton method for image reconstruction in diffuse optical tomography[END_REF][START_REF] Abascal | Fluorescence diffuse optical tomography using the split bregman method[END_REF][START_REF] Ducros | Reconstruction of an optical inhomogeneity map improves fluorescence diffuse optical tomography[END_REF]. In this paper, it is implemented using explicit computation of the Hessian, which is a sparse and structured matrix. Hence, our algorithm can benefit from efficient linear equation solvers. The formulation is presented in a general framework, providing analytical formulas for general regularization functionals. Our algorithm, which is referred to as regularized weighted least square Gauss-Newton (RWLS-GN), is validated on synthetic data and the projection of a voxelized thorax phantom based on a clinical CT volume. We provide results for different scenarios created by varying the number of incident photons and the amount of contrast agent. To promote reproducibility of this research, the computer code (MATLAB scripts) related to this study is available online [41].

II. THEORY

As depicted in figure 1, we consider a 3-dimensional (3-D) object in a bounded domain Ω that is imaged by a 2-D detector with a surface S. A summary of the notation used throughout the paper is provided in table I. Here, a conical geometry is chosen as an illustration. Ω is the object volume, S is the detector surface, ρm is the mass density of the mth material, and am is projected mass density of the mth material.

A. X-ray/matter interaction

Let n(E, u) denotes the number of photons of energy E that reaches the detector at the pixel position u ∈ S. Neglecting scattering within the object, n follows the Beer-Lambert law. Mathematically,

n(E, u) = n 0 (E, u) exp - L(u) µ(E, x) d ( 1
)
where n 0 (E, u) is the source spectrum, L(u) is the acquisition line that depends on the source emission geometry (parallel, cone-beam, etc), and µ(E, x) is the local LAC of the object at energy E and point x ∈ Ω. For simplicity, the source spectrum is assumed to be the same at each pixel, i.e.,

n 0 (E, u) = n 0 (E), ∀u ∈ S (2)

B. Detection model

The mean signal recorded by a photon counting detector may be modelled by

s(E, u) = R d(E, E) n(E, u) dE (3)
where d(E, E) accounts for the detector response function, E being the incident photon energy and E the detected energy. The photons detected within the ith energy bin [E i , E i+1 ] are accumulated electronically thanks to a counting circuit. The mean number of photons detected within the i-th energy bin is given by si

(u) = R d i (E) n(E, u) dE (4)
where

d i (E) = Ei+1 Ei d(E, E) dE (5)
is the response function of the i-th bin of the detector.

Note that the detection model assumes that the value measured at a given pixel is not correlated to the value measured at neighbor pixels, i.e. charge sharing can be neglected or is corrected for [START_REF] Koenig | How spectroscopic x-ray imaging benefits from inter-pixel communication[END_REF].

C. Object decomposition

It is assumed that the LAC can be decomposed as the sum of M basis functions that are separable in energy and space. We have

µ(E, x) = M m=1 ρ m (x)τ m (E), ∀x ∈ Ω ( 6
)
where the τ m are some well-chosen basis functions and the ρ m are the corresponding weights in the decomposition. Following the work of Alvarez and Macovski [START_REF] Alvarez | Energy-selective reconstructions in x-ray computerised tomography[END_REF], two approaches have emerged for the choice of the basis functions: i) physics-based where τ m models the physical effects, e.g. photoelectric, Compton scattering, and, ii) material-based where τ m is the mass attenuation of the constituents of the objects (in cm 2 .g -1 ). In the latter approach, ρ m is the mass density (in g.cm -3 ) of material m. For the sake of simplicity but with no loss of generality, ρ m is referred to as mass density in the following. K-edge materials may be added with the two types of basis [START_REF] Schlomka | Experimental feasibility of multi-energy photon-counting k-edge imaging in pre-clinical computed tomography[END_REF].

D. Forward model

Substituting ( 6) into (1) and the obtained result into (4), the mean measured signal may be written

si (u) = R n 0 (E)d i (E) exp - M m=1 a m (u)τ m (E) dE (7)
where

a m (u) = L(u) ρ m (x) d (8)
is the projection of the mass density ρ m along the line L(u) (see figure 1) that is expressed in g.cm -2 when a material-based decomposition is performed. In the following, a m is referred to as the projected mass density (PMD).

E. Discretization

Let us assume that the detector consists of an array of P pixels located at u p , p ∈ {1, . . . , P }. Each of the pixels counts photons in I different energy bins. Let s ∈ R P I and s ∈ R P I be the mean and measured photon counts vector, which are defined by s = s1,1 . . . sI,1 . . . . . . s1,P . . . sI,P

s = s 1,1 . . . s I,1 . . . . . . s 1,P . . . s I,P .

In a similar way, let a ∈ R P M be the (unknown) vector containing the PMD of each material in each pixel, which is defined by

a = a 1,1 . . . a M,1 . . . . . . a 1,P . . . a M,P . (11)
The discrete forward model can be formalized as

s = F (a) (12)
where F (a) denotes the non-linear mapping induced by [START_REF] Wang | Mi-croCT with energy-resolved photon-counting detectors[END_REF].

The material decomposition problem is to recover the PMD vector a from the measured data s given the known forward mapping F . This is an ill-posed problem that requires prior knowledge about the solution to stabilize the inversion in the presence of noise.

F. Noise model

The measurements are assumed to be corrupted by independent Poisson noise, i.e., s i,p ∼ P (λ = si,p ) ,

where P (λ) denotes the Poisson distribution of mean λ.

In order to investigate the influence of noise which is related to the value of si,p , we introduce the total number of photons emitted by the source during the acquisition. It is defined by

N 0 = n 0 (E) dE (14)
The larger N 0 is, the better the signal-to-noise ratio of the measurements.

III. PROPOSED DECOMPOSITION METHOD

A. Variational approach

In this manuscript, we propose to recover the projected mass vector within a variational framework, i.e., minimizing the cost function q for bq and as |bq| -for bq where W is a weighting matrix, R is the regularization functional, and α is a global regularization parameter. The weighting matrix is chosen as

C(a) = 1 2 s -F (a) 2 W + αR(a) (15)
W = diag 1 s (16)
such that s -F (a) 2 W is a weighted least squares fidelity term adapted to data corrupted by Poisson noise.

One main difficulty of the inverse problem is the non linearity of the forward mapping, which is addressed by means of the Gauss-Newton algorithm described in Section III C.

B. Regularization functional

Priors must be suited to the images of materials to be recovered. For instance, a soft tissue image is expected to be smoother than a bone image. We adopt the following regularizing functional, which allows the prior of the different materials to be tuned independently:

R(a) = M m=1 β m R m (b m ), with b m = L m a m , ( 17
)
where β m , R m , a m = [a m,1 . . . a m,P ] ∈ R P , and L m ∈ R Qm×P are the regularization parameter, the regularization functional, the PMD, and a linear transform of the m-th basis function (e.g., the derivative or the laplacian), respectively.

Many different regularizing functionals, together with many different transforms have been considered in the past. In this manuscript, we limit our study to the following cases.

a. Material regularization functional We consider the functionals that can be modeled by

R m (b m ) = Qm q=1 ψ(b q ) (18)
where ψ : R → R + is a twice differentiable potential function chosen to promote some image prior and b q = (b m ) q is the q-th component of b m . In particular, the following two potential functions, illustrated in figure 2, are considered:

• Quadratic. The quadratic potential function

ψ quad (b q ) = b 2 q (19)
can be chosen to promote smooth solutions. It leads to the standard 2 penalization also known as Tikhonov regularization:

R quad m (b m ) = b m 2 2 (20)
• Huber. The Huber penalty functional can be chosen to preserve edges in the solution. In particular, its smooth approximation is defined by

ψ huber m (b q ) = b 2 q + 2 -. (21)
where is a parameter that can be tuned knowing the order of magnitude of the b q 's. This choice for ψ results in a 1 -like penalization term:

R huber m (b m ) ≈ b m 1 (22)
b. Linear transform The choice for the linear operator L m is restricted to discrete differential operators up to second order, i.e., L m ∈ {D (0) , D (1) , D (2) } [START_REF] Long | Multi-material decomposition using statistical image reconstruction for spectral CT[END_REF] where D (0) ∈ R P ×P is the identity matrix, D (1) ∈ R 2P ×P is the 2-D discrete gradient, and D (2) ∈ R P ×P is the 2-D discrete Laplacian. Any other linear operators could be used as well. However, these choices lead to sparse L m matrices, which constitute a computational advantage as detailed in Section VI.

C. Optimisation algorithm

We propose to minimize (15) using Gauss-Newton (GN)'s method, which is a classical iterative method for non-linear minimization. It starts with an initial guess a (0) and builds new estimates with the following update rule:

a (k+1) = a (k) + λ (k) δa (k) . ( 24
)
where δa (k) ∈ R P M is the so-called GN step and λ (k) ∈ R is the step length.

The GN step is obtained solving

(J (k) WJ (k) + αH (k) )δa (k) = -g (k) , (25)
where J (k) is the Jacobian matrix of F at the iterate a (k) , H (k) is the Hessian matrix of R at a (k) , and g (k) is the gradient of C at a (k) that is given by:

g (k) = -J (k) W 1 2 s -F (a (k) ) + α∇R(a (k) ). ( 26
)
The step length is obtained minimizing the cost function along the direction provided by the GN direction. Mathematically, it is defined as

λ (k) = arg min C(a (k) + λ δa (k) ) (27)

D. Computation of the Jacobian

The Jacobian matrix of the forward model introduced in ( 25) is the IP × M P real-valued matrix that maps variations in the unknown vector a onto variations in the forward model s = F (a). It is defined by

J = ∂s ∂a (28)
Differentiating the forward model given in [START_REF] Wang | Mi-croCT with energy-resolved photon-counting detectors[END_REF], it can be shown that ∂s i,p ∂a m,p = 0, for p = p j i,m,p , for p = p [START_REF] Maaß | Empirical multiple energy calibration (emec) for materialselective CT[END_REF] with

j i,m,p = - R τ m (E)n(E, u p )d i (E) dE. ( 30
)
As a result, the Jacobian can be expressed as the block diagonal matrix

J = diag J 1 . . . J P (31)
where

J p ∈ R I×M is such that (J p ) i,m = j i,m,p .
It is noteworthy that J is a sparse matrix having only IM P non-zero entries out of IM P 2 . Physically, the many zeros of J indicate that a change of the PMD in a particular pixel does not affect the value measured at another pixel.

E. Computation of the gradient of the regularizer

Computation of the gradient of the regularizer R is required by [START_REF] Barber | An algorithm for constrained onestep inversion of spectral CT data[END_REF]. By linearity of the gradient, we obtain from the definition of the regularizer [START_REF] Huy | Least squares parameter estimation methods for material decomposition with energy discriminating detectors[END_REF] that

∇R(a) = M m=1 β m ∇R m (L m a m ), ( 32
)
where ∇ is the gradient with respect to a, which is defined by

∇ = ∂ ∂a1,1 . . .
and to remark that the full gradient ∇R m ∈ R M P is simply obtained by inserting zeros into the partial gradients

∇ m R m ∈ R P . Precisely, ∇R m = ∇ m R m ⊗ e m ,
where e m ∈ R M is the natural basis vector given by e m = [δ 1,m . . . δ M,m ] where δ m ,m is the Kronecker delta. Therefore, we have

∇R(a) = M m=1 β m ∇ m R m (L m a m ) ⊗ e m . ( 35
)
The latter equation is of particular interest when closed-form formulas for ∇ m R m are available. Applying the appendix formula (A2) to the potential functions considered in ( 19) and ( 21), the following simple identities are obtained:

∇ m R quad m (b m ) = 2L m b m (36)
and

∇ m R huber m (b m ) = L m b1 √ b 2 1 + 2 . . . b Q √ b 2 Q + 2 (37)
where b q = (L m a m ) q .

F. Computation of the Hessian of the regularizer

The Hessian of the regularizer is required to compute the GN update that solves [START_REF] Zhao | An extended algebraic reconstruction technique (E-ART) for dual spectral CT[END_REF]. It is defined by:

H = ∂ 2 R ∂a∂a . ( 38
)
With the same reasoning as in the previous section, we have

H(a) = m β m H m (b m ) ⊗ e m ⊗ e m . ( 39
)
where H m ∈ R P ×P is a shorthand for the Hessian of R m w.r.t. b m . At this point it is important to notice that H is built from the M matrices H m , each of size P × P . Moreover, closed-form equations for H m are obtained applying the appendix equation (B2) to the potential functions given in [START_REF] Faby | Performance of today's dual energy CT and future multi energy CT in virtual non-contrast imaging and in iodine quantification: A simulation study[END_REF] and [START_REF] Mendonca | A flexible method for multi-material decomposition of dual-energy CT images[END_REF]. We have :

H quad m (b m ) = 2L m L m (40)
(b m ) = 2 L m diag (b 2 1 + 2 ) -3 2 . . . (b 2 Q + 2 ) -3 2 L m (41)
where b q = (L m a m ) q .

IV. NUMERICAL SIMULATIONS AND IMAGE ANALYSIS A. Acquisition parameters

The source spectrum n 0 (E) was simulated using the spekCalc software [START_REF] Gavin | Calculation of x-ray spectra emerging from an x-ray tube. part ii. x-ray production and filtration in x-ray targets[END_REF][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. It is depicted in the top row of figure 3. Measurements were performed in I = 4 energy bins. The response function d i (E) of each of the 4 bins was computed according to [START_REF] Gimenez | Characterization of medipix3 with synchrotron radiation[END_REF], where d(E, E) was simulated using the parametric model provided by Schlomka et al. [11, (see appendix A.2)] and choosing the following energy thesholds: E 1 = 15 keV , E 2 = 36 keV, E 3 = 60 keV, and E 4 = 91 keV (see figure 3). In order to investigate the robustness of our decomposition method to noise, different total numbers of photons were considered in the range N 0 ∈ [10 4 , 10 7 ] photons.

The projections were simulated with RTK [START_REF] Rit | The reconstruction toolkit (RTK), an open-source cone-beam CT reconstruction toolkit ba-sed on the insight toolkit (ITK)[END_REF] using a point x-ray source at 1.5 m from a 41 × 41 cm 2 flat panel with 256 × 256 pixels. The ray casting used bi-linear interpolation in each slice of the phantom following and 16 × 16 rays were calculated and averaged per pixel.

B. Phantom generation

We used the thoracic-abdominal CT scan depicted in figure 4a), which is part of the 3D-IRCADb database [START_REF]3d image reconstruction for comparison of algorithm database[END_REF]. It is a 480 × 370 × 167 volume with a 0.961 × 0.961 × 1.8 mm 3 voxel size and an 8-bit grayscale. We also used the anatomical atlas represented in figure 4d). It was obtained and provided by Kéchichian et al. [START_REF] Kéchichian | Shortest-path constraints for 3d multiobject semiautomatic segmentation via clustering and graph cut[END_REF] who segmented the previously mentioned CT scan via vicinity prior. Given both the CT and atlas volumes, the two material mass density volumes depicted in figure 4b) and c) were created as follows. First, segments were grouped. For the three-material phantom, all the segments corresponding to soft tissue [START_REF] Namely | lung (R), lung (L), liver, spleen, gallbladder, pancreas, kidney (R), kidney (L), adrenal gland (R), adrenal gland (L), hyperplasia, aorta, vena cava, portal vein, skin, fat, visceral fat, muscle, marrow, renal pyramids (R), renal pyramids (L), digestive tract tissue, digestive tract cavity, pulmonary arteries (R), pulmonary arteries (L), bronchi(R), bronchi (L), intervertebral discs, hepatic arteries, hepatic vena cava, hepatic portal veins[END_REF] were grouped into the same segment Ω soft . For the five material phantom, Ω soft did not include adipose tissue [START_REF]Namely, fat and visceral fat[END_REF] and lung tissue [START_REF][END_REF] that were assigned to the segments Ω fat and Ω lung , respectively. For both phantoms, the bone segment Ω bone remained unchanged. Second, the mass density in each voxel was estimated from the CT images. A linear mapping between the 8-bit grayscale of the CT scan and the mass density was assumed, i.e., ρ m = λ µ CT . The constant λ (in g.cm -3 ) was calibrated such that mean value of the mass density of visceral fat is set to 0.95 g.cm -3 . Third, the portal vein segment Ω portal was marked a uniform concentration c gd of gadolinium, which is a promising spectral CT contrast agent [START_REF] Feuerlein | Multienergy photon-counting k-edge imaging: Potential for improved luminal depiction in vascular imaging1[END_REF][START_REF] Sawatzky | Proximal admm for multi-channel image reconstruction in spectral x-ray CT[END_REF]. The previous three steps result in the computation of three volumes that were each projected onto the detector according to [START_REF] Wang | Material separation in x-ray CT with energy resolved photon-counting detectors[END_REF] and the acquisition geometry L(u) described in Section IV A. The resulting projected material mass density image is displayed in figure 4d).

The mass attenuations of soft tissue, bone, and gadolinium, as required in [START_REF] Wang | Mi-croCT with energy-resolved photon-counting detectors[END_REF], were taken from the ICRU report 44 [START_REF]International Commission on Radiation Units and Measurements (ICRU), Report 44 of the Tissue Substitutes in Radiation Dosimetry and Measurement[END_REF]. They are depicted in figure 3 (bottom row).

C. Image quality metrics

The quality of the decomposed images has been assessed by means of two different objective performance metrics. First, the normalized euclidean distance has been considered. Given the ground-truth image a true m and the decomposed image a dec m , the normalized Euclidean error for the m-th material is defined as

error m = a dec m -a true m 2 a true m 2 . ( 42
)
The lower the normalized Euclidean error, the better. We also define the total normalized Euclidean error as

error tot = 1 M M m=1 error m . (43)
Second, the contrast-to-noise ratio (CNR) has been considered. Given a region of interest S roi ∈ S and a background region S back = S \ S roi , we define the CNR of the decomposed image a dec m as

cnr m = āroi -āback (w roi σ 2 roi + w back σ 2 back ) 1/2 , (44)
where āroi (ā back ) and σ roi (σ back ) are the average and standard deviation of a dec m in the region-of-interest S roi (background S back ) and w roi (w back ) is the partial surface of S roi (S back ) in the image. The larger the CNR, the better.

D. Algorithm parameters

The quality of the reconstructed images depends dramatically on the choice of the regularization parameters α introduced in [START_REF] Pan | Computed tomography in color: Nanok-enhanced spectral ct molecular imaging[END_REF]. In order to fairly compare decompositions obtained from different acquisition scenarios, each decomposition was performed for a set of regularization parameters in the range log(α) ∈ [-2, 2], with a 0.5 step. The best regularization parameter α has been chosen to obtain the lowest reconstruction error, for the phantom with three materials, and then the same regularization parameter has been used to assess the method for different phantom configurations. For the selection of the material specific regularization parameters β soft , β bone , β gd , we evaluated the algorithm for values between 0.1 and 10 and found that the simple choice β soft = β bone = β gd = 1 provided results close to the best ones (see supplemental data in table V at the end of the manuscript). The hyperparameter was set to a small value = 0.01.

Our RWLS-GN algorithm has two stopping criteria. It is interrupted either if i)the length step is small enough, i.e., λ (k) < λ stop or ii) the cost function decrease is small enough, i.e., 1 -C (k) /C (k-1) < δ stop . We chose λ stop = 5 × 10 -3 and δ stop = 10 -3 . The initial guess a 0 was chosen as constant images with values in the range of the expected values for each material, i.e., a 0 soft = 10 g.cm -2 , a 0 bone = 1 g.cm -2 , and a 0 gd = 0 g.cm -2 .

V. RESULTS

A. Comparison between ML-NM and RWLS-GN methods

Our algorithm was compared to the reference method described in [START_REF] Schlomka | Experimental feasibility of multi-energy photon-counting k-edge imaging in pre-clinical computed tomography[END_REF], which relies on the maximization of Ground truth images of the PMD of soft tissue, bone and gadolinium are shown in the top row of figure 5. For the best case scenario (number of incident photons of 10 7 and concentration of the marker 1 g.cm -3 ), decomposed images provided by both the RWLS-GN and the ML-NM methods are given in the middle row of figure 5. Difference between decomposition and ground truth images is shown in the bottom row of figure 5. While materials were well separated when using RWLS-GN, the images decomposed using ML-NM have poor quality. Large artifacts are observed in the soft and bone tissue images and the gadolinium image is very noisy. For ML-NM, the largest errors are related to the recovery of gadolinium, which results in a dramatic cross-talk with the soft and bone tissue images. For RWLS-GN, it can be seen that the gadolinium image presents the lowest error and that the largest error occurred at borders of the spine. Quantitatively, RWLS-GN led to error soft = 0.014, error bone = 0.271, error gd = 0.071, and cnr gd = 3.42 for α = 10 -0.5 , while ML-MN led to error soft = 1430, error bone = 1558, error gd = 15, and cnr gd = 0.3627.

For the considered 256 × 256 images, a RWLS-GN iteration requires about 2s on a standard laptop (2.6 GHz i7 CPU and 16 GiB of RAM) which leads a decomposition time of about 20s assuming 10 iterations before convergence as observed in all the assessed experiments. The ML-NM is much slower with a decomposition time of about 1500s for the same images.

B. Convergence of the proposed RWLS-GN algorithm

Convergence of the algorithm is shown in figure 6. The method converged in around ten iterations. The step length was above one in almost all iterations and the decrease of the cost function was close to one in all except for the last iterations. Number of iterations for convergence was between 10 and 15 for all scenarios considered in this work.

C. Influence of the regularization parameter

The influence of the regularization parameter α is illustrated in figure 7. Figure 7a shows the solution error for the different materials as a function of α. The lowest error for all materials was attained for similar values of α in the range [0.1, 1]. Figure 7b and 7c shows underregularized and overregularized images, respectively. When regularization is not sufficient, images are noisy and the method fails to decompose materials. On the contrary, excessive regularization leads to oversmoothened images, failing to properly decompose soft and bone tissue. Overall, the gadolinium image presented the best image quality and was more robust to the selection of the regularization parameter, and separation of bone and soft tissue were more sensitive to the selection of α.

The optimal regularization parameter for all scenarios is shown in table II. For a sufficiently large number of photons (N 0 ≥ 10 5 ) the optimal value log(α * ) was found to be almost independent of the number of photons and concentration of the marker (log(α * ) = 0 or -0.5, see table II). For a low number of photons (N 0 ≤ 10 5 ) the method converged to unsatisfying solutions. In addition, it was found that values of α in the range 0.1 and 1 led to similar results in terms of image decomposition and recovery of the concentration of the marker.

For the rest of the analysis in this work, the regulari-

(α ) X X X X X X X X X N0 c gd (g.cm -3 ) 1 0.3 0.1 0.03 0.01 10

D. Influence of noise

The influence of noise has been assessed by varying the number of incident photons from N 0 = 10 7 to N 0 = 10 4 for a concentration of the marker of c gd = 1 g.cm -3 . For N 0 = 10 7 materials were well separated and quantification of the marker concentration was correct (see figure 5). Lowering photons from N 0 = 10 6 to N 0 = 10 5 decreased image quality for bone and soft tissue while the gadolinium image was less affected (see top and middle rows in figure 8). For a number of photons N 0 = 10 4 , material decomposition was not feasible.

These results are confirmed by the Euclidean distance error (first column in table III) and contrast-to-noise ratio of the marker (first column in table IV). Lowering the number of photons from N 0 = 10 7 to N 0 = 10 5 increased the error between two and three times but has limited influence on the contrast-to-noise ratio of the marker.

E. Influence of the contrast agent marker

The influence of the contrast agent marker was studied by varying the marker concentration c gd from 1 to 0.01. Decomposed gadolinium images for different marker concentration and number of photons are shown in figure 9. Marker concentrations of 1 g.cm -3 were well recovered even for a low number of photons (N 0 = 10 5 ). Lowering the marker concentration to 0.1 had a large effect on image quality even for a large number of photons, which is indicated by an increased decomposition error. For concentrations below 0.01 g.cm -3 , recovery of the marker was not feasible. Contrast-to-noise of the decomposed marker shows that the recovery of the marker was more affected by marker concentration than by noise (see table IV). Contrastto-noise was similar to the best case scenario for c gd = 1 when N 0 ≥ 10 5 , c gd = 0.3 when N 0 ≥ 10 6 and for c gd = 0.1 when N 0 ≥ 10 7 .

F. Influence of the phantom

Robustness of the algorithm with respect to the phantom is shown in figure 10. Results correspond to different phantom compositions (soft tissue, bone tissue, and gadolinium only or together with lung and adipose tissue) and projection angles. In all cases, three materials (soft tissue, bone tissue and gadolinium) are considered for decomposition. Decomposition of the sagittal view is found to be feasible considering the algorithm parameters used to decompose the coronal view. For a phantom with five materials (soft tissue, bone, lung tissue, adipose tissue and gadolinium) decomposed images are visually similar to those from a phantom with only three materials (compare figure 5b and 10a and figure 10b and10c). While the bone error is found to be higher in a five-material phantom, it is interesting to note that the gadolinium error is of the same order in both phantom (about 0.07 for the coronal view and 0.2 for the sagittal view).

VI. DISCUSSION AND CONCLUSION

This work presented a new and efficient method to solve the nonlinear decomposition problem in spectral X-ray imaging. The problem was solved using a Gauss-Newton step that exploited sparsity of Jacobian and Hessian operators.

It has been validated on numerical experiments of a thorax phantom and has been compared to the ML-NM algorithm. We showed that the proposed RWLS-GN method solves the nonlinear decomposition problem in an efficient manner. The algorithm converged in around ten iterations for all scenarios. It was also shown to improve the accuracy of the projected material maps. Decomposed material images presented high image quality, low error and accurate concentration of the marker. On the contrary, ML-NM decomposition was highly affected by noise, leading to artifacts in bone and soft tissue images and to a very noisy gadolinium image. For the assessed experiments, RWLS-GN decomposition was feasible for a number of incident photons larger or equal than N 0 = 10 5 and for a marker concentration above 0.03 g.cm -3 . In the context of spectral CT, our method could exploit data redundancy in the projection domain and be parallelized across projections. Most previous studies in spectral CT included regularization in the tomographic problem [START_REF] Schirra | Statistical reconstruction of material decomposed data in spectral CT[END_REF][START_REF] Sawatzky | Proximal admm for multi-channel image reconstruction in spectral x-ray CT[END_REF][START_REF] Xu | Sparsity-regularized image reconstruction of decomposed k-edge data in spectral CT[END_REF]. These methods decompose the data in a previous step and then perform the tomographic problem by regularizing in the image domain, which is a simpler but computationally more expensive problem as it requires to project and backproject the data at each iteration. In contrast, in this work we tackle regularizing the nonlinear decomposition problem, which has not been widely addressed. This problem is nonlinear, harder to solve and it presents some differences but it also offers some advantages that can be exploited. The two main advantages are that the problem is smaller (i.e. one projection at a time) and that the forward operator (and also gradient and Hessian operators) is very sparse. For instance the number of nonzero elements in the Hessian is of the order of the number of pixels instead of the square of the number of pixels. The proposed algorithm exploits these features by using sparse matrices to explicitly represent gradient and Hessian operator and a Gauss-Newton step, which has a higher convergence rate than the most widely used linearized methods. At each iteration, the linear system (25) can be efficiently solved using a Cholesky decomposition. Fast solvers for the factorization of sparse matrix are available. Their complexity depends on the sparsity (percentage of non-zero entries) of the matrix and also the sparsity pattern (location of the non zero entries). These solvers are perfectly adapted to the resolution of (25) since the matrix H = 2J J + αH is both sparse and structured. Indeed, it follows from (31) that J J = diag J 1 J 1 . . . J P J P , i.e., J J is block diagonal with a sparsity of M 2 P/M 2 P 2 = 1/P . Considering P = 256 × 256 pixels, the sparsity (portion of non-zero values) of J J is 0.0015 %. Moreover, the matrix H is build from the products of sparse banded matrices H m , hence it is also sparse banded.

This work has focused on the decomposition of spectral X-ray images acquired under a particular projection angle and may find direct applications in interventional radiology. In this context, it could avoid the need of image subtraction that is subject to motion artifacts. In the context of spectral CT, its main advantage is that it can be parallelized, processing each projection in parallel. This has the potential to provide a fast method for sinogram decomposition. Although it has been shown that RWLS-GN decomposition is feasible for different projection views and that it improves standard ML-NM decomposition, its extension to spectral CT requires further investigation.

One of the main differences with previous studies in the literature is that most of the approaches regularize in the spatial domain and very few works have addressed regularization in the projection domain. The proposed algorithm included regularization on the projection domain and allowed the use of specific regularization functional for each material. The problem was solved using a Gauss-Newton step with efficient computation of operators. The selection of the regularization parameter was studied and the improvement brought by regularization was highlighted for different experimental scenarios. The gain of including regularization in the material decomposition step was shown by comparing the proposed algorithm to ML-NM, which led to unsatisfactory decomposed images highly affected by noise and artifacts.

The proposed method allows the use of different regularization functionals for each material and it is formulated in a general framework for any choice of the operator. In this study, we selected Total Variation for the contrast material and first-order Tikhonov and secondorder Tikhonov for bone and soft-tissue, respectively. These choices are based on assumptions that gadolinium projection images should be piecewise constant and projection image for soft tissue is smooth. Nevertheless, other choices of regularization functionals could be more suitable and provide better results, and should be investigated in future work.

FIG. 1 .

 1 FIG.1. Principle of acquisition and notations. Here, a conical geometry is chosen as an illustration. Ω is the object volume, S is the detector surface, ρm is the mass density of the mth material, and am is projected mass density of the mth material.

FIG. 2 .

 2 FIG.2. Potential functions used in this manuscript. The quadratic potential function leads to the 2, while the Huber function approximates the 1-norm . The Huber function is depicted for = {0.5, 0.2, 0.05, 0.02}. Note that ψ huber (bq) behaves as1 2

∂

  ∂a M,1 . . . . . . ∂ ∂a 1,P . . . ∂ ∂a M,P . (33) Since regularization is performed on each material independently, i.e. R m only depends on a m , we have ∂Rm ∂a m ,p = 0, ∀m = m. It is therefore convenient to introduce the partial gradient w.r.t. the PMD of the m-th material ∇ m = ∂ ∂am,1 . . . ∂ ∂a m,P .

FIG. 3 .

 3 FIG.3. Top row: Spectrum of the x-ray source. Middle row: detector response for each energy bin. Bottom row: mass attenuation of the three materials considered in the numerical experiments (gadolinium in green, bone in grey, lung tissue in red, soft tissue in yellow, adipose tissue in cyan; lung and soft tissues are superimposed).

FIG. 4 .

 4 FIG. 4. Numerical phantom. (a) Slice of the 3D-IRCADb [46] thoracic-abdominal CT scan (in Hounsfield units). (b) Slice of the 3-material phantom showing the mass density ρm (g.cm -3 ) of soft tissue in yellow, bone in gray, and gadolinium in green. (c) Slice of the 5-material phantom showing the mass density ρm (g.cm -3 ) of soft tissue in yellow, bone in gray, gadolinium in green, lung tissue in red, and adipose tissue in blue. (d) Anatomical atlas obtained segmenting the 3D-IRCADb volume via vicinity prior [47], different colors indicate different labels. (e) Projection of the 3-material phantom in the coronal plane. The projected material mass density am (g.cm -2 ) is displayed with the color code used in (b). (f) Projection of the 5-material phantom in the coronal plane. Same color code as in (c).

FIG. 5 .

 5 FIG. 5. Typical material decompositions. The projected mass densities (g.cm -2 ) are obtained considering the 3-material phantom shown in figure 4b/4e. The gadolinium concentration c gd was set to 1 g.cm -3 the number of incident photons N0 to 10 7 photons. (a) Ground truth images. (b) Images recovered by our GN-RWLS method for α = 10 -0.5 (best regularization parameter). (c) Images recovered by the maximum likelihood Nelder-Mead (ML-NM) method described in [11]. (d) GN-RWLS error. (e) ML-NM error.

  parameter α has been set to 10 -0.5 .

FIG. 6 .

 6 FIG. 6. Convergence of the algorithm for α = 10 -0.5 , c gd = 1 g.cm -3 , and N0 = 10 7 photons. (Top left) cost function C as a function of the iteration number k. (Top right) total normalized Euclidean error errortot as a function of the iteration number k. (Bottom left) step length λ (k) vs iteration number k. (Bottom right) cost function decrease C (k-1) -C (k) C (k-1) vs iteration number. The red dashed lines indicate the values of the chosen stopping criteria.

FIG. 7 .

 7 FIG. 7. Influence of regularization. Measurements were simulated for N0 = 10 7 photons and c gd = 1 g.cm -3 . (a) Normalized euclidean distance as a function of the regularization parameter α. (b) Underregularized images recovered for α = 10 -4 . (c) Overregularized images recovered for α = 10 2 .

FIG. 8 .

 8 FIG. 8. Influence of noise. Decomposed images for decreasing number of incident photons: N0 = 10 6 photons (top), N0 = 10 5 photons (middle), and N0 = 10 4 photons (bottom). The marker concentration was fixed to 1 g.cm -3 .

4 FIG. 9 .

 49 FIG. 9. Decomposed gadolinium images obtained for decreasing marker concentrations c gd = {1, 0.1, 0.01} g.cm -3 (from left to right) and decreasing number of incident photons N0 ∈ {10 7 , 10 6 , 10 5 } from top to bottom. In all cases α = α .

FIG. 10 .

 10 FIG. 10. Robustness of the decomposition setting α = 10 -0.5 . From left to right, recovered soft tissue, bone, and gadolinium images. (a) Decomposition of the 5-material phantom (coronal view) leading to error bone = 0.697, error gd = 0.066, and α = 10 -0.5 . (b) Decomposition of the 3-material phantom (sagittal view) leading to error soft = 0.019, error bone = 0.330, error gd = 0.212, and α = 10 -1 . (c) Decomposition of the 5material phantom (sagittal view) leading to error bone = 0.701, error gd = 0.201, and α = 10 -1 . In all cases, measurements were simulated for N0 = 10 7 photons and c gd = 1 g.cm -3 .

TABLE I .

 I Main mathematical notations

	Variable Meaning	Dimension
	a	Projected mass density (our unknown) R(S)
	d	Detector response function	R(R, R)
	i	Output energy bin index	N
	j	input energy bin index	N
	k	iteration index	N
	m	Basis functions index	N
	n	Number of input photons	R(R, S)
	s	Number of output photons	R(R, S)
	E	Input energy	R
	I	Number of output energy bins	1
	J	Number of input energy bins	1
	M	Number of basis functions	1
	P	Number of pixels	1
	µ	Linear attenuation

TABLE II .

 II Best regularization parameter. The table displays log

TABLE III .

 III Total normalized Euclidean distance errortot for different scenarios

	X N0 X X	X c gd (g.cm -3 ) X X X X X 10 7	1 0.12	0.3 0.13	0.1 0.21	0.03 0.49	0.01 1.26
		10 6	0.19	0.21	0.33	0.75	1.97
		10 5	0.33	0.39	0.61	1.45	4.00

TABLE IV .

 IV Contrast-to-noise ratio of the decomposed marker image.

	X N0 X X	X c gd (g.cm -3 ) X X X X X 10 7	1 3.42	0.3 3.41	0.1 3.05	0.03 1.86	0.01 0.82
		10 6	3.48	3.34	2.73	1.31	0.56
		10 5	3.66	2.91	1.82	0.65	0.24

TABLE V .

 V (Supplemental data). Quality of the decomposition for different hyperparameters {β soft , β bone , β gd } and the best global regularization parameter α . In the simulations, we set N0 = 10 7 photons and c gd = 1 g.cm -3 .

	β soft	β bone	β gd	log(α )	error	cnr gd
	1	1	1	-0.5	0.12	3.42
	1	0.1	0.1	0	0.12	3.42
	1	0.1	10	-1	0.17	3.64
	1	10	0.1	-1	0.14	3.33
	1	10	10	-1.5	0.13	3.40
	0.1	1	0.1	0	0.14	3.35
	0.1	1	10	-1	0.15	3.68
	10	1	0.1	-1	0.12	3.41
	10	1	10	-1	0.12	3.52
	0.1	0.1	1	0	0.13	3.57
	0.1	10	1	-0.5	0.16	3.42
	10	0.1	1	-0.5	0.13	3.45
	10	10	1	-1.5	0.12	3.39

ACKNOWLEDGMENTS

This work was supported by the grant ANR-11-INBS-0006 of the French National Research Agency (ANR) and was performed within the framework of the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program ANR-11-IDEX-0007. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement number 701915.

The authors have no relevant conflicts of interest to disclose.

Appendix A: Gradient of the regularizer

We seek an expression for the gradient of the regularizer R(b), b = La, with respect to a. Applying the chain rule leads to:

where ∇ a is the gradient w.r.t. the variable a and ∇ b is the gradient w.r.t. b. Differentiating R(b) = q ψ(b q ) with respect to b, we obtain

Appendix B: Hessian of the regularizer

Applying the chain rule twice, it may be shown that the Hessian of R(b) with respect to the variable a, for a linear change of variable b = La, is given by (see equation 2.1 of Lebl [START_REF] Lebl | Tasty Bits of Several Complex Variables[END_REF])

where


Agnes Jocher

Heinz Pitsch

Thomas Gomez

Guillaume Legros
email: guillaume.legros@upmc.fr

Modification of sooting tendency by magnetic effects

Keywords: sooting tendency, laminar non-premixed flame, magnetic field

= 88 Fig. 1 = 407 Fig. 2 = 99 Fig. 3 = 495 Fig. 4 = 280 Fig. 5 = 145 1

Introduction

Major air pollution disasters in industrialised cities, e.g. London smog disaster in 1952, have highlighted the harmful effects of pollutant emissions by combustion devices on human health. Increased risks of cardiac arrest, pulmonary diseases, pregnancy complications and asthma are now associated with the presence of particulate matter, especially fine-particulate soot, in urban air and higher mortality rates [1][2][3].

Although soot formation and oxidation processes in non-premixed flames have been extensively investigated, they are not fully understood due to the numerous chemical and physical phenomena that drive them. Still, there is now a general understanding that soot production along the combustion process is the result of the following major competing processes: pyrolysis of fuel leading to the formation of soot molecular precursors, nucleation, growth of particle mass via surface reactions and coagulation, and oxidation.

Oxidation mainly proceeds through the removal of carbon atoms from soot particles, by reaction with O 2 and OH radicals [4,5]. While O 2 is considered the most important species in the soot oxidation process, the role of OH radicals is more controversial. Conducting numerical simulations, Frenklach [4] stated that oxidation by OH radicals is relatively weak, while Kim et al. [5] showed that these radicals are highly soot-oxidising species.

Standard strategies of soot production control involve post-treatment methods, such as filters [6] and/or catalysts [7]. However, using in-situ strategies looks like an increasingly promising method. Examples are the design of new fuels [8], addition of exhaust gases to the combustion process [9],

or combustion control by electrical fields [10]. Interestingly, none of these strategies utilises magnetic fields.

Yamada et al. [11] observed experimentally and numerically that the O 2 and OH radicals concentration fields are highly sensitive to the gradient of the square of the magnetic flux density (∇(B 2 )). Among the species contributing to the combustion process, O 2 and OH indeed exhibit relatively high paramagnetic susceptibilities. The magnetic force on the species is directly proportional to the susceptibilities, but indirectly proportional to the temperature. As OH radicals mainly appear in the hot regions of the flame, the magnetic force on OH is negligible small compared to the one acting on O2 in the coflow. Shinoda et al. [12] studied a premixed methaneair flame and found that the magnetic force acts mainly on O 2 , due to its higher mass fraction. The high sensitivity of OH radicals on ∇(B 2 ), shown by Yamada et al. [11], is mainly a result of the modified convection velocity induced by the magnetic force on O 2 .

While such a magnetic force can also play a significant role on the stability of non-premixed flames [13,14], its influence on the sooting behaviour of nonpremixed flames has never been assessed. Consequently, the present study is a contribution to an original in-situ strategy of soot production control by magnetic fields. A steady, laminar, non-premixed ethylene flame is established over the axisymmetric burner, prescribed by Santoro et al. [15]. This configuration is selected due to the available literature about soot production, which delivers both numerical and experimental databases. A non-intrusive optical technique, based on line-of-sight attenuation measurements, provides the maps of soot volume fraction in the flames.

As a result, the modification of the inner flame structure by magnetic effect, is reported in terms of soot volume fraction.

Experimental setup

Figure 1 exhibits a schematic of the experimental setup.

Burner configuration

The non-premixed flames are established over an axisymmetric coflow burner, identical to the one described by Santoro et al. [15]. The burner's axis of symmetry is vertical. This burner was also used in previous investigations [14,16]. The configuration has delivered a well documented literature on the sooting behaviours of a wide range of fuels [17,18].

The ethylene stream flows through a vertical axial brass duct, which has an effective diameter of injection d F , of 11 mm. The coflowing oxidizer mixture is introduced into a concentric brass cylinder of 102 mm inner diameter.

The oxidizer consists of a O 2 and N 2 mixture. Two Bronkhorst EL-FLOW mass flow controllers enable the variation of the oxygen molar content, X O 2 , in the oxidizer stream. Over the range of investigated flow rates, both fuel and oxidizer flows are steady and laminar at the height of the burner tip.

Further details about the burner can be found in Ref. 14.

In the following, the burner's axis of symmetry is (Oz) and its origin is located at the burner tip, defining the height above the burner (HAB). The cross-stream coordinate is r, which is the distance from the axis of symmetry.

In order to expose the flames to fairly uniform upward ∇(B 2 ), the burner tip is located 0.13 m below the coils' axis of symmetry.

None of the rig components is magnetic to avoid mechanical interferences.

The rig enables the adjustment of the burner position inside the electromagnet.

Magnetic field generation

The electromagnet used is dedicated to combustion studies, however slightly different in size from the one described in Ref. 14. It consists of two identical coils of inner radius R = 0.12 m, that face each other, to create the magnetic field. Both coils have a common horizontal axis and are separated by a distance similar to R. An equal, constant, electrical current, I, flows through the turns that constitute each coil. Such a configuration increases the volume of the maximum magnetic field magnitude that a single coil can create and promotes the field uniformity within the system core. An iron magnetic circuit then confines the magnetic field lines.

A two-dimensional numerical simulation has been carried out to compute the magneto static field. A set of measurements by a Hall effect sensor enable the calibration of the computed fields. Details of this procedure are given in Ref. 14. As the isolevel lines in the air gap between both coils are found to be fairly horizontal, the ∇(B 2 ) is directed towards the air core, i.e. upwards at the location of the burner tip. Furthermore, for any current in the electromagnet, the absolute gradient peak appears as a plateau, located at heights between 0.11 m and 0.14 m below the coils' axis. Therefore, within a cylinder of 50 mm radius, whose axis is the burner's axis of symmetry (i.e. containing the coflow at the burner exit), included between the aforementioned heights, the gradient magnitudes are fairly uniform.

This academic configuration will provide clean conditions to be incor-porated into future numerical simulations. These will then allow further optimisation of a practical magnet design.

The level of magnetic gradient is adjusted by a power supply that delivers the current flowing into the coils up to I max = 60 Amperes. The z-component magnitudes of the gradient together with the computed spatial fluctuations within the cylindrical domain and the current required to generate the different gradients are reported in Tab.1.

The horizontal component magnitudes of the gradients were checked to be negligible. Therefore, the gradients' direction and the burner's axis of symmetry are identical. The absence of magnetic field is called MagF 0.

Diagnostics

The Laser Extinction Measurement (LEM) technique, which provides a two-dimensional soot volume fraction field, has been described in details in Ref. 16. It has been shown to be capable, of providing both fine temporal and spatial resolutions. The arrangement of the optical diagnostics to perform LEM, specifically through a flame located in the electromagnet, is shown in As a line-of-sight technique, LEM needs to be combined with a subsequent deconvolution method, to compute the local soot volume fraction f v (r, z).

An Onion-Peeling method is performed here, followed by a Tikhonov regularisation, which stabilised the deconvolution process [16]. This procedure allows the local extinction coefficient κ ext λ (r, z) to be inferred from the extinction measurements provided by the camera and to be integrated over the line-of-sight. The Mie theory then allows the distribution of κ ext λ (r, z) to be transformed into the soot volume fraction, assuming that particles are in the Rayleigh limit:

f v (r, z) = λ κ ext λ (r, z) 6 π E(m) (1)
where E(m) is a function of the complex refractive index m of soot. To overcome the controversial issue of the evaluation of this parameter, together with the relative magnitude of scattering in the extinction phenomenon, we proceeded as Arana et al. [19] by adjusting E(m), to reproduce the peak soot volume fraction, measured by Santoro et al. [15] at HAB = 50 mm. The adjusted E(m) (0.42) is close to the one found by Arana et al. [15] (0.38) who performed a 2D-LEM at a slightly lower wavelength (632.8 nm). The procedure then induces a discrepancy lower than 8 % in magnitude and 5 % in location over the remaining peaks.

Results and discussion

Constants and parameters

The following investigations are restricted to steady-state flames established at atmospheric pressure and constant room temperature of 295K.

Ethylene is chosen as fuel for its extensively documented sooting be-

Assessment of the magnetic effect

The laser beam intensity captured over the CMOS sensor without any flame was first checked not to be influenced by any level of the magnetic field investigated.

The influence of the three different levels of the ∇(B 2 ) is first assessed on the raw integrated extinction coefficients measured inside the ethylene non-premixed flame burning into air. In this way, no modification due to the magnetic field can be attributed to the post-processing of the raw LEM data. Figure 2 shows the distributions of the integrated extinction coefficient as a function of the magnetic gradient that the flame experiences. For MagF 0, the peak integrated extinction coefficient (0.3) is located on the flame's edge at about HAB = 40 mm. From MagF 1 to MagF 3, the peak stays on the edge and monotonically increases in magnitude, reaching 0.32, 0.36, and 0.39, respectively. A similar increase is globally observed in the whole flame.

The oxygen molecules are attracted towards the axis of symmetry of the two iron coils, whereas diamagnetic species like N 2 or C 2 H 4 are rejected by the increasing magnetic field [START_REF] Wakayama | [END_REF]. Nevertheless, the absolute value of the oxygen magnetic susceptibility per unit weight can be assumed several times larger than those of diamagnetic gases [START_REF] Yamaguchi | Magneto-Science[END_REF]. Consequently, the main influence of the magnetic field on the non-premixed flame is expected to be caused by the oxygen concentration field. The flame is here submitted to an upward gradient of the square of the magnetic flux density. As explained by Yamaguchi and Tanimoto [START_REF] Yamaguchi | Magneto-Science[END_REF], the resulting thermo-magneto convection is therefore opposed to buoyancy. For relatively long flames, buoyancy plays an important role. However with increasing ∇(B 2 ), the thermo-magneto convection is enhanced. Consequently the global residence time inside the flame is increased with increasing ∇(B 2 ). Furthermore, the thermo-magneto convection will also reduce the velocity of the fresh oxidizer flow. The shear stress between the oxidizer and fuel flow will therefore be weakened. Both effects then promote soot formation. Indeed, higher integrated extinction coefficients are reported both on the edge and on the axis of the flame when the magnitude of the upward magnetic gradient is increased.

At constant fuel and oxidizer flow rates of 3.83 cm 3 /s and 717 cm 3 /s, respectively, the oxygen content in the oxidizer stream is varied from 21 % O2 (air) to 50 %. Figure 3 exhibits the distributions of the soot volume fraction for different values of the magnetic gradient and for different oxygen contents. To provide comparable maps, every distribution of the soot volume fraction is shown in a non-dimensional form, the reference value being the peak f v,max obtained among the distributions at a given oxygen content.

These reference values are reported on the right of every row.

At any magnetic gradient, when increasing the O 2 content, the flame length is reduced due to enhanced transport of oxygen towards the fuel side of the flame. In Fig. 3, this effect is somewhat less obvious, due to the variable scale of the z-axis. The peak soot volume fraction is then increased with increasing oxygen content since in hotter and smaller flames, the soot formation rate is hastened on the rich side of the flame.

At a constant O 2 content, the influence of the magnetic gradient on the non-premixed flame is significant but ambiguous. For X O 2 =21 %, the peak soot volume fraction is increased and shifted towards the flame axis when increasing the magnetic gradient. However, for higher X O 2 , the peak soot volume fraction is shifted away from the flame axis and now decreases in magnitude with increased magnetic gradient. This is in agreement with the longer characteristic residence time mentioned above. Yet, this decrease is balanced as the soot layer also spreads slightly over when increasing the magnetic gradient.

To further elucidate these observations, net soot production is now characterised by the integrated soot volume fraction F v (z) defined as follows:

F v (z) = 1 π (d F /2) 2 d F /2 0 2 π r f v (r, z) dr (2)
F v (z) appears as an evaluation of the actual soot load at a given HAB in the flame. Its maximum value F v,max can therefore be considered a relevant characteristics of the competition between soot formation and oxidation along the globally upward trajectories that particles experience. For this reason, the height z max at which F v,max occurs is also a characteristic length of the soot formation and oxidation processes. Figure 4 then exhibits the profiles of non-dimensional F v /F v,max as a function of non-dimensional height z/z max .

For every oxygen content, a specific F v,max is defined as the maximum value of F v for MagF 0 and the specific z max is the height where F v,max is reached.

For X O 2 =21 %, the upper parts of the profiles are missing because the soot volume fraction fields are not fully encompassed into the laser beam section.

However, the central parts around F v,max are well characterised.

At a given magnetic gradient, the profiles for the different oxygen contents investigated look very similar, further supporting z max and F v,max as characteristic quantities of soot related processes into these non-premixed flames.

Per se, this finding is not marginal. Among others, Roper [22] showed that non-premixed flame sheets are similar. However, soot particles contribute to phenomena that are essentially volumetric due to the relatively long characteristic times of formation and oxidation. The flame sheet assumption can then be considered all the more questionable when the flame is sooty.

Nonetheless, the quantity F v (z) appears here to decently follow a similarity law.

At any oxygen content, the effect of the magnetic gradient on the net soot production does not reverse as the aforementioned ambiguous trends may have led to conclude. Indeed, increasing the magnetic gradient monotonically enhances soot production for every oxygen content. Similar trends have been reported in non-buoyant non-premixed flames [START_REF] Legros | Proc. Combust. Inst[END_REF], which further supports the above reasoning about the balance between buoyancy and thermomagneto convection. Interestingly, when increasing the oxygen content, the soot production enhancement by the magnetic gradient is emphasised and can be considered significant. Therefore, the modification of soot production by magnetic gradients could constitute a relevant strategy of reacting flows control, especially in the field of oxyfuel combustion where oxygen content is relatively high.

Magnetic scaling of soot production

As an illustration of such a strategy, we investigated the potential to scale the soot production by the magnetic gradient. To address a practical concern, we sought for reproducing the profile of non-dimensional integrated soot volume fraction experimentally measured for Mag 0 (I = 0) and X O 2 = 50 % with a lower oxygen content, i.e. X O 2 = 30 %, and with the help of the magnetic gradient. This strategy would potentially lead to lower the oxygen flow rate in an oxyfuel combustion device.

To this end, the database composed of the experimental z/z max (0.5 ≤ z/z max ≤ 1.2), F v /F v,max , X O 2 , and I/I max was approximated by an interpo-lation function G interp , as schematically represented by Eq.( 3):

F v /F v,max = G interp (z/z max , X O 2 , I/I max ) (3)
G interp was finally a second-order polynomial which was found to adequately model the experimental F v /F v,max with an error range lower than 3 %. The ten coefficients of this polynomial were determined using a trustregion-reflective algorithm [START_REF] Coleman | [END_REF] that minimises the least-squares between the results of Eq.( 3) and the set of experimental F v /F v,max .

We then selected a specific z/z max , i. where soot formation dominates soot oxidation. Downstream z/z max = 1, soot oxidation starts prevailing, potentially due to higher levels of OH radicals [5]. Yet, these radicals also exhibit relatively high magnetic susceptibility.

Therefore, a sophisticated coupling between the transport of these radicals and the thermo-magneto convection conveying oxygen may reveals, damaging the similarity of the F v /F v,max profiles.

Conclusion

The present study may be considered an original and significant step towards the identification of a possible strategy that would enable the magnetic control of soot production in a non-premixed flame. The original experimental setup included the coflow burner prescribed by Santoro, together with an optical arrangement that allows the measurements of two-dimensional soot volume fraction fields in the flame. Thus, the influence of the magnitude of an upward magnetic gradient on the soot production inside the ethylene non-premixed flame was investigated.

Increasing the magnitude of the magnetic gradient allowed the modification of soot production in the flame. This modification is enhanced with increasing oxygen content, as oxygen exhibits a relatively high paramagnetic susceptibility. Furthermore, this modification was shown to enable the shift among similar profiles of soot concentration in the flame, as the variation of oxygen content does.

This latest finding may drive further optimisation of the strategy. CO2 addition to the coflowing air will induce a decrease of the adiabatic flame temperature due to the higher specific heat of CO 2 compared to that N 2 . Since the Curie law predicts a magnetic susceptibility that is inversely proportional to the temperature, this should globally enhance the magnetic effect. For this reason, the scope of the strategy could be extended to the field of oxyfuel combustion.

While reversing the magnetic field itself has no effect on ∇(B 2 ) and hence on the soot volume fraction maps, a reduction of soot production might be possible when shifting the burner to the upper core of the electromagnet to meet downwards ∇(B 2 ). For these further investigations, the whole setup will be carefully modified.

Fig. 1 .

 1 Fig.1. The system consists of a 100 mW Continuous Wave Laser operating at 645 nm (-5/+7 nm), a set of beam expanding optics, collection optics, and a camera. A digital pulse generator controls the occurrence and the duration of the CMOS exposure, together with the shutter opening. A frame grabber records the frames, captured by the camera, on a computer.The Photon Focus MV1 12-bit progressive scan monochrome camera is mounted with a conventional lens, equipped with a narrow band filter centred at 645 nm (± 2 nm), as well as with a band width, at one half of the trans-

  haviour, especially on the Santoro's configuration. A single ethylene flow rate of 3.83 ± 0.07 cm 3 /s is established, corresponding to a mean fuel injection velocity of V F = 3.98 ± 0.09 cm/s, i.e to a fuel exit Reynolds number of Re F = 33. The oxidizer mixture flow rate is set at a single value of 717 cm 3 /s, corresponding to a mean oxidizer velocity of V ox = 8.9 ± 0.09 cm/s.The oxidizer Reynolds number Re ox is then about 300. These flow rates are chosen so that the flames presented in this study are contained within the domain of almost uniform ∇(B 2 ). Furthermore, all flames are laminar, steady, and non-smoking.The study consists of investigating flames that are burning in oxidizer streams of four different O 2 contents, i.e. X O 2 = 0.21, 0.30, 0.40, and 0.50 (± 0.01). Every flame experiences each of the three upward ∇(B 2 ) reported in Tab. 1.

  e. 0.87, into the experimental profile for X O 2 = 50 % and MagF 0. At this non-dimensional height, F v /F v,max is equal to 0.97. Intending to reproduce the profile for X O 2 = 50 % and MagF 0, we solved Eq.(3) for I/I max with z/z max = 0.87, F v /F v,max = 0.97, and X O 2 = 30 %. The computed current was I/I max = 0.167, i.e. I=10.0 amperes. This current induces a magnetic gradient of about 1.1 T 2 /m. We finally established the flame for X O 2 = 30 % and I = 10.0 Amperes and measured the soot volume fraction field. Figure 5 exhibits the targeted profile F v /F v,max measured for X O 2 = 50 % and I = 0, together with the one measured for X O 2 = 30 % and I = 10.0 Amperes. The close up view in the upper right corner shows that the specific point targeted is fairly reproduced in the flame for X O 2 = 30 % and I = 10.0 Amperes. This strategy assessment looks here especially relevant in the lower part of the profiles, i.e. in a region

Figure 1 :

 1 Figure 1: Schematic view of the experimental setup. This allows soot volume fraction measurements in an axisymmetric ethylene non-premixed flame experiencing different magnitudes of ∇(B 2 ). The insert in the upper left corner exhibits the different profiles of ∇(B 2 ) as a function of HAB. As an illustration, the flame established for X O2 = 50 % and MagF 0 is shown in reversed intensity scale.

Figure 2 :

 2 Figure 2: Maps of the raw integrated extinction coefficient κ int λ for increasing ∇(B 2 ) magnitudes, starting from MagF 0 to MagF 3. The ethylene non-premixed flame burns here in air.

Figure 3 :Figure 4 :

 34 Figure 3: Non-dimensional soot volume fraction maps for increasing magnetic gradient from the left to the right, and for increasing oxygen content in the oxidizer stream from the top the the bottom. For every oxygen content, the reference f v,max is reported on the right of the corresponding row.

Figure 5 :

 5 Figure 5: Profiles as a function of non-dimensional height z/z max and non-dimensional integrated soot volume fraction F v /F v,max . The experimental profile for X O2 = 30 % and I = 10.0 Amperes intends to reproduce the experimental profile for X O2 = 50 % and I = 0 Ampere at the specific location z/z max = 0.87 and F v /F v,max = 0.97. The matching between both profiles around this location is shown in the close up view in the insert.

Table 1 :

 1 ∇(B 2 ) magnitudes.

	Magnetic field	MagF 1 MagF 2 MagF 3
	Average peak (T 2 /m)	1.8	7.1	18.2
	Spatial fluctuation (T 2 /m) ± 0.01	± 0.1	± 0.3
	Current (A)	13	26	60

Acknowledgments

The authors acknowledge the Université franco-allemande/Deutsch-Französische Hochschule for their financial support. The authors feel very grateful to M. Kashif, J. Bonnety, H. Dutilleul, and J.M. Citerne for their technical support.