S Anthoine
email: anthoine@cmi.univ-mrs.fr

J.-F Aujol

Y Boursier
email: boursier@cppm.in2p3.fr

C Mélot
email: melot@cmi.univ-mrs.fr

Some proximal methods for CBCT and PET tomography

Cone-Beam Computerized Tomography (CBCT) and Positron Emission Tomography (PET) are two complementary medical imaging modalities providing respectively anatomic and metabolic information of a patient. In the context of public health, one must address the problem of dose reduction of the potentially harmful quantities related to each exam protocol : X-rays for CBCT and radiotracer for PET. Two demonstrators based on a technological breakthrough (acquisition devices work in photon-counting mode) have been developed and we investigate in this paper the two related tomographic reconstruction problems. We formulate separately the CBCT and the PET problems in two general frameworks that encompass the physics of the acquisition devices and the specific discretization of the object to reconstruct. These objects may be observed from a limited number of angles of views and we take into account the specificity of the Poisson noise. We propose various fast numerical schemes based on proximal methods to compute the solution of each problem. In particular, we show that primal-dual approaches are well suited in the PET case when considering non differentiable regularizations such as Total Variation. Experiments on numerical simulations and real data are in favor of the proposed algorithms when compared with the well-established methods.

1

).

We will describe in the following the general model of acquisition of a CT-scan.

The PIXSCAN II demonstrator

. Hybrid pixel detectors were developed originally for particle tracking in high energy physics experiments. They are actually very promising X-ray detectors for

Introduction

Cone Beam Computerized Tomography (CBCT) and Positron Emission Tomography (PET) scans are medical imaging devices that provide respectively anatomical and metabolic complementary information about the patient.

From the point of view of modelization they are key-examples of ill-posed inverse problems. Indeed the observations in both cases are obtained by a given number of projections of the original data corrupted by noise. The projection operator is generally highly ill-conditionned and non invertible. The noise takes into account additive dark noise (gaussian noise), from the electronic of the measurement device, and Poisson noise from the counting process to get the measurement.

Moreover on the one hand CBCT scans use the intrinsic physical property of biological tissues, since they absorb X rays. One would like to lower the dose delivery necessary to get an image, and this could be done by improving the detection efficiency of the scanner. On the other hand a PET scan is based on the principle of injection of active radiotracer to the patient and one would like also to reduce this dose. Remark furthermore that the informations provided by a CBCT scan are needed in the reconstruction of data from a PET scan (see explanations below). This is why the CCPM in Marseille, France, developed recently a new prototype of small animal scanner ClearPET/XPAD. It is designed to record images for both modalities simultaneously. Moreover it includes the use of a new generation of detectors which are not affected by dark noise, and thus lead to model the measurements with pure Poisson noise.

Taking into account this new setting, the goal of this work is to apply reconstruction algorithms both in CBCT and TEP cases for pure Poisson noise. With this model building numerical algorithms to reconstruct the data can be anyway challenging. Recent developments in the field of convex optimization allow to use proximal methods, and in particular primal-dual algorithms [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], which solve in a very accurate way these inverse problems in tomography. To our knowledge some of them were never applied to tomography. We will test their robustness when increasing the noise and lowering the number of projections used to reconstruct.

The contributions of the papers are the following:

• We propose a complete analysis of both a new TEP and a new CBCT scanner. Due to their new technology (use of hybrid pixels, as detailed in section 2), this lead us to some new functionals to minimize.

• By making use of recent developments in convex optimization, we can minimize the functional we propose without the need to smooth them (as far as we know, this is new in the case of the TEP functional involving a logarithm term).

• We present extensive numerical results, so that the reader can judge of the efficiency of the proposed methods.

• Since we have prototypes of the new scanners, we can test the modeling and the algorithms introduced in the paper on real datas (and not just on synthetic examples). Our analysis of the scanners really starts from the modeling of the physics of the acquisition device, and ends at the testing on real datas.

The outline of this paper is the following. In Section 2 we describe the physics and explain the mathematical modeling in both cases. In Section 3 we review quickly the state of the art in inverse problems related to tomography. From Section 4 to Section 6 we explain recent developments about proximal, and especially primal-dual algorithms in the convex optimization community. We detail the algorithms we are going to use. In section 7, we present numerical results on simulations and real datas.

2 Physics of CBCT and PET imaging 2.1 Cone-Beam Computerized Tomography

Introduction

Computerized Tomography (CT) is a medical imaging modality that provides anatomical information on contrast images. A CT-scan basically consists in a X-ray source and a X-ray camera with the imaged object placed in between. For a given position of the source and the detector relative to the object, an image is acquired and it consists in a projection along the line-of-sight defined by the source-detector direction. A full set of acquisitions is then obtained by stacking successive images taken for different positions during the rotation of the block source-detector around the object. The data thus acquired, that are the input of the tomographic reconstruction, express the value of the projection as a function of the projection angle, and it is called a sinogram.

The function to reconstruct in the case of the CT-scan is the absorption coefficient, denoted by µ, of the different tissues that constitute the object. This coefficient is energy-dependent and its values are tabulated (see figure 1 which gives some examples). For instance, the main organs of the human body are made of a mixture of tissues, and each organ is characterized by a single coefficient µ. The X-ray source delivers a cone-beam and rotates around the patient lying on a bed. The X-ray detectors are located on a ring all around the patient. material science and micro-computed tomography [START_REF] Debarbieux | Repeated imaging of lung cancer development using pixscan, a low dose micro-ct scanner based on xpad hybrid pixel detectors[END_REF]. They are based on the silicon hybrid pixel tracker developed for the Atlas collaboration at CERN. The use of hybrid pixels is expected to help to improve the contrast for soft tissues and to reduce both the scan duration and the dose absorbed by the animal.

Hybrid pixel detectors are used to count X-ray photons above a given energy threshold. The XPAD3 camera is composed by 130 µm 2 pixels with a readout every microsecond without dead time. The readout is independent between pixels. There is no charge amplification in the readout process, so that the XPAD3 camera is not subject to dark noise like usual charge integration detectors. Notice that dark noise that usually affects any CCD-like integration charge detector is classically modeled by an additive gaussian noise.

Thus, hybrid pixel detectors allow acquisition of low statistics data with a pure Poisson noise due to unavoidable photon noise. We therefore expect that photon counting CBCT will permit to lower drastically the exposition dose.

Let us now quickly describe the PIXSCAN II microCT system, displayed on Figure 3. It is basically composed of three blocks which are :

1. The first block contains a Oxford Instruments X-ray source which has a W target and a 13 to 40 µm focal spot size, and it can be operated up to 90 kV and 2 mA. A filter wheel is placed in front of the point of emission in order to modify the spectrum of the source according to the filter used. This source delivers a cone-beam characterized by a 33 degrees apex angle.

2.

The second block consists in a rotating platform placed between the X-ray block and the detector block. The animal that is imaged stands on the platform (fixed to an animal holder), while the platform rotates step by step (typically 360 deg., 1deg/step). Three additional motors allow for translation in the three main directions.

3. The third block supports the X-ray hybrid pixels XPAD3/Si camera. The detectors implement a 500 µm thick silicon sensors and it is made of 8 barrettes of 7 chips each, assembled in a tile geometry. Each chip contains 9600 hybrid pixels (80 × 120) so that the full detector contains a half million 130x130 µm 2 pixels. The full camera size is 125x75 mm 2 .

General and monochromatic acquisition models

Let us modelize the reconstruction of images from CBCT scan.

The probability p that a photon at a given energy E is not absorbed by the material along a line-of-sight (los) L is expressed thanks to the Beer-Lambert law by : This probability only depends on the µ coefficient that reflects anatomical properties of the human body.

p = exp - L µ(l, E)dl . (1)
In the general case, let us define an object by its absorption coefficient µ = µ(r; E), where r stands for a position vector in R 3 , and y L by the measurement of a CT-scanner.

The measurement is made counting all the photons that have not been absorbed by the imaged object, thus by integrating all the µ coefficients along the los L and for all energies E according to the Beer-Lambert law in equation 1 :

y L = ∞ 0 z L (E) exp - L µ(l, E)dl dE (2)
where z L (E) = I(E) • S(E) with :

• S(E) : the detector efficiency at energy E.

• I(E) : the number of photons emitted by the source at energy E along the los L.

In this paper, we will exclusively focus on a monochromatic model, i.e. we assume that the X-ray source delivers a monochromatic beam at a given energy E 0 , which implies the following equations:

z L (E) = z L (E 0 ) • δ(E -E 0 ) (3)
The acquisition model is then simplified:

y L = z L (E 0 ) exp - L µ(l, E 0 )dl (4)

Discretization of the problem

Let us consider the discrete framework that is relative to the real-world aspects : objects to reconstruct are discretized into voxels and CT cameras are pixellized. First, if a CT camera contains N pixels (indexed by n = 1, . . . N ), a tomographic set of measurements y is obtained with Θ angles of projection (indexed by θ = 1, . . . Θ) so that y ∈ R J 1 is indexed by j = 1, . . . , J 1 with J 1 = N Θ.

Secondly, the continuous object µ to reconstruct is discretized according to a decomposition on a finite dimension basis {b i } i∈I 1 such that µ ≈ i∈I 1 µ i b i . The choice of {b i } i∈I 1 is critical and may be done by minimizing the approximation error or by offering an adapted representation to the objects to reconstruct (see [START_REF] Hanson | Local basis-function approach to computed tomography[END_REF] for more results).

However, for numerical reasons, the family function of voxels is most of the time chosen by the users.

The object to reconstruct is then a vector µ ∈ R I 1 indexed by i = 1, . . . , I 1 . Equation ( 4) then becomes :

y j = z j (E 0 ) exp - I 1 i=1 a ij µ i (E 0 ) (5)
where z j (E 0 ) represents the number of photons emitted by the source into the solid angle defined by the pixel j seen from the source point. In the following, we will write z j (E 0 ) = z j for simplicity. In this problem, (a i,j ) I 1 ×J 1 , (z j ) I 1 , (y j ) J 1 and the unknown vector (µ i ) I 1 are nonnegative constants. More precisely, a ij = 0 when the ray j does not intersect the voxel i, and 0 < a ij ≤ 1 when the ray j intersect the voxel i, depending the definition of the linear operator A = (a ij ) I 1 ×J 1 ∈ M(R I 1 , R J 1 ), called the system matrix. This operator A is a numerical implementation of the operators of projection that fully describes the geometry of the acquisition system. The coefficient a i,j of A typically characterizes the probability that any event occurring on a photon in pixel i will be detected on pixel j.

Incorporating noise

We model the measurements as independently distributed Poisson random variables and we can assume the general model of measurements :

y j ∼ P z j exp -[Aµ] j + r i (6)
where P(λ) describes a Poisson distribution with λ parameter. Additive classical detector read-out noise can be accounted for in several ways. Since new generation photon-counting detectors are not affected by dark noise classically mod-eled by additive Gaussian noise, there is no additive classical detector read-out noise and we have ∀i r i = 0:

y j ∼ P z j exp -[Aµ] j . (7)
We recall that the Poisson likelihood reads:

P (Y = y|X = x) = x y e -x y! . ( 8
)
From Bayes rule, we then obtain that P (µ|y) = P (y|µ)P (µ) P (y) [START_REF] Brunner | Imaging performance of the hybrid pixel detectors xpad3-s[END_REF] so that by keeping the µ-dependent part, we obtain that the negative log-likelihood writes for our photon-counting detectors:

L CT (µ) = - J 1 j=1 {y j log z j exp -[Aµ] j -z j exp -[Aµ] j } (10)
Since we want to minimize this negative log likelihood we can work with

L CT (µ) = J 1 j=1 {y j [Aµ] j + z j exp -[Aµ] j } (11)
This problem is ill-posed, and we add a regularization term J(µ) (described later) to the data fidelity term L CT . Thus we consider the following problem:

μ = arg min µ≥0 L CT (µ) + λJ(µ). (12)
In this problem, λ is a regularization parameter that makes the balance between the datafidelity term and J(µ) which is a cost function that introduces some apriori information on the object to reconstruct.

Positron Emission Tomography

Introduction

Unlike the CBCT, Positron Emission Tomography (PET) is a medical imaging modality that provides a measurement of the metabolic activity of an organ. The exam consists in injecting to the patient and further during the exam a radiotracer attached to a molecule, for instance an analog of glucose like the [START_REF] Dewaraja | Regularized reconstruction in quantitative spect using ct side information from hybrid imaging[END_REF]-FDG, that will be absorbed preferentially by some organs depending on their function. The measurements made in PET modality are the results of a radioactive decay that emits a positron. This positron will annihilate with an electron after a very short time of flight in the human body (the mean free path is around 1mm). During the positronelectron annihilation, two gamma rays of 511keV are emitted in opposite directions, and several rings of detectors around the patient are detecting these gamma rays.

The ClearPET/XPAD demonstrator

The ClearPET/XPAD small animal PET/CT scanner design is based on the combination of the ClearPET scanner [START_REF] Lartizien | The clearpet project: Development of a 2nd generation high performance small animal pet scanner[END_REF], built at EPFL, Lausanne within the Crystal Clear Collaboration with the PIXSCAN microCT prototype, built at CPPM, Marseille (see Figure 4). The ClearPET main characteristics are its detectors constituted of LSO/LuYAP phoswich crystal matrices read-out by multi-anode photomultiplier tubes and its partial ring geometry, which allows for inserting an X-ray tube and the hybrid pixel camera which together constitute a micro-CT scanner [START_REF] Mosset | Development of an optimized lso/luyap phoswich detector head for the lausanne clearpet demonstrator[END_REF][START_REF] Rey | Count rate performance study of the lausanne clearpet scanner demonstrator[END_REF]. Both modalities can therefore share the same field of view. The positron tomograph is a modified version of the ClearPET camera built in Lausanne. Indeed the new prototype is built with a new partial ring geometry and it uses an upgrade of the data acquisition software. A 150mm lead shielding has been added in order to protect the PET modules against scattered X-rays. The PET partial ring geometry consists of 3 rings of 7 detector modules as shown. This geometry suppresses axial gaps between rings and allows acquisition of a complete set of 3D sinograms while optimizing the sensitivity. This configuration has a 55 mm axial field of view and 112 mm transverse field of view. The combined PET/CT system inserts modification of both scans to use a better arrangement of the detectors and to allow for the integration of the micro-CT elements on the same rotating gantry. It has been used to demonstrate the feasibility of simultaneous bi-modal acquisitions [START_REF] Khodaverdi | Design study for the ClearPET/XPAD small animal PET/CT scanner[END_REF][START_REF] Nicol | Design and construction of the clearpet/xpad small animal pet/ct scanner[END_REF].

Acquisition model

Let us define the concentration of activity v = v(r) that we want to reconstruct, i.e. the concentration of radioactive decay. We also define w L as a measurement acquired by the sensors of the PET-scan along a line of sight L, i.e. the number of gamma rays detected.

The absorption coefficient µ at 511keV of the material that are lying along the los is taken into account to model the possible absorption of photons emitted at an energy of 511keV. For the PET reconstruction problem, we will consider that µ is known. We obtain the following acquisition model in a continuous framework :

w L = L v(l) exp (-µ(l, E 511 ))dl, (13)

Discretization of the problem

The discretization of the PET acquisition can be deduced from Equation ( 13) by assuming that the concentration activity is discretized into a number of voxels equal to I 2 , so that the vector to reconstruct is v ∈ R I 2 . In the same way as for the CT-scan, we assume that J 2 measurements, indexed by j = 1, . . . , J 2 , are acquired by Θ cameras uniformly located in angles around the imaged object, and each composed of N pixels so that J 2 = N Θ. We finally obtain :

w j = I 2 i=1 c ij v i exp (-c ij µ i (E 511 )), (14)
where the matrix C = (c ij ) I 2 ×J 2 ∈ M(R I 2 , R J 2 ) incorporates the parallel-beam geometry of the acquisition, similarly to the system matrix A described for the CT-scan problem. Depending on the numerical implementation of the projection and backprojection, the coefficient c ij may typically describe the length of intersection between the voxel i and the los joining the source of the gamma emission to the pixel j. If we define the final system matrix B for the PET-scan problem that incorporates the µ(E 511 ) correction by ∀i, j, b ij = c ij exp (-c ij µ ij (E 511 )), we deduce the final discrete acquisition framework :

w j = I 2 i=1 b ij v i , (15)

Incorporating noise

As the PET-scan measurements are a counting process, we model the data acquisition of PET-Scan in a similar way as the CBCT problem by incorporating a pure Poisson noise :

w j ∼ P([Bv] j ), (16)
where v ∈ R I 2 denotes the radioactive concentration vector to reconstruct, w ∈ R J 2 the vector of measurements and B the system matrix which describes the full properties of the PET-Scan (J 2 < I 2 ). The negative log-likelihood that we want to minimize then reads :

L P ET (v) = J 2 j=1 [Bv] j -w j log ǫ ([Bv] j ) (17)
with the notation: log ǫ (x) = log(x + ǫ). Since this problem is also ill-posed, we add a regularization term J to the data fidelity term L and we finally consider the following PET-reconstruction problem :

v = arg min v≥0 L P ET (v) + λJ(v) (18)
3

State of the art

The two problems we address are actually different, but algorithms to solve them come from the same seminal family.

Some reconstruction algorithms start by applying a variance-stabilizing transform such that the Anscombe's transform [START_REF] Andscombe | The transformation of poisson, binomial and non negativebinomial data[END_REF] or the Haar-Fisz transform [START_REF] Fisz | The limiting distribution of a function of two independant random variables and its statistical application[END_REF] to go closer to the hypothesis of gaussian noise [START_REF] Dupé | A proximal iteration for deconvolving poisson noisy images using sparse representations[END_REF] and [START_REF] Fryzelewicz | A haar-fisz algorithm for poisson intensity estimation[END_REF]. Since for small values of the parameter in the Poisson noise, these transforms do not yield exactly gaussian noise, in the following, we will keep the hypothesis of Poisson noise and we will deal directly with the data fidelity terms [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] and [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] derived from the physic of the problem.

In both of the problems of CBCT and PET a Radon transform is involved. Thus a set of method deals with the use of filtered back projection (FBP) to recover the unknown datas [START_REF] Feldkamp | Practical cone-beam algorithm[END_REF]. The idea of filtered backprojection is to use the Fourier slice theorem and the definition of the dual Radon transform (the so called "backprojection") to get an inversion formula. This amounts in 2 dimensions to filter the Radon transform of the data by a ramp filter and to backproject the data (see [START_REF] Helgason | Groups and Geometric Analysis: Integral Geometry, Invariant Differential Operators, and Spherical Functions[END_REF] for more details and mathematical proofs). The ramp filter has the property to filter out low frequencies and to emphasize high frequencies, which let recover well sharp edges but in the same time accentue the high frequencies of the noise. To deal with this issue, this is a classical method to use filters with some cut off frequency. In [START_REF] Mariano-Goulart | A priori selection of the regularization parameters in emission tomography by fourier synthesis[END_REF], the authors focus in the setting of PET on this specific point. They propose to solve it by minimization of a least square functional, penalized with a quadratic term related to the cut off frequency of the filter. This amounts in their case to solve a linear system. Their algorithm yields low pass filtered solutions.

Another set of methods is based on the Expectation-Maximization algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], the so called EM or ML-EM algorithm (see [START_REF] Shepp | Maximum likelihood reconstruction in positron emission tomography[END_REF] for PET and [START_REF] Lange | EM reconstruction algorithms for emission and transmission tomography[END_REF] for both CT and PET). This iterative algorithm seeks the maximum of the log likelihood function. Its main drawback is its numerical cost. Furthermore, since the problem is anyway ill posed, the estimation is known to be unstable even if the algorithm garantees that it will always satisfy positivity constraints. One often prefers to stop iterations before reaching convergence in order to get an acceptable result. To improve and stabilize convergence, works focus on acceleration procedures and on penalizing the data fidelity terms with an additional term. Following the approach of De Pierro [START_REF] De Pierro | A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography[END_REF], Fessler and all [START_REF] Erdoǧan | Monotonic algorithms for transmission tomography[END_REF] develop an algorithm for the PET problem using quadratic surrogates functions to compute the minimum of the -log likelihood penalized by a quadratic constraint. Iterations are accelerated by using Ordered Subsets algorithms [START_REF] Hudson | Accelerated image reconstruction using ordered subsets of projection data[END_REF], or the so called block iterative-incremental gradient methods. Unfortunately these algorithms may cycle and one needs to build relaxed versions or use surrogates functions to insure convergence [START_REF] Ahn | Globally convergent image reconstruc-tion for emission tomography using relaxed ordered subsets algorithms[END_REF]. Among regularizations Huber regularization [START_REF] Huber | Robust Statistics[END_REF] is sometimes preferred [START_REF] Dewaraja | Regularized reconstruction in quantitative spect using ct side information from hybrid imaging[END_REF] since it is differentiable. Popular is also the median root prior, which leads generally to modified EM iterations [START_REF] Alenius | Bayesian image reconstruction for emission tomography based on median root prior[END_REF] but, written as a penalization term added to the data fidelity term, is non convex. Total variation (TV) as a regularization has the main advantage to preserve well edges but leads to a convex non differentiable penalization term, which is challenging from the point of view of algorithms and convergence proofs. It has been used in the context of the CT problem changing the non quadratic log likelihood into a quadratic one [START_REF] Sidky | Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[END_REF]. This approximation is actually plausible as soon as the approximation of the Poisson's noise by gaussian noise is feasible. Making use of quadratic surrogates functions together with TV penalization or also ℓ 1 regularization (ℓ 1 ), the authors in [START_REF] Harmany | This is SPIRAL-TAP: Sparse poisson intensity reconstruction algorithms theory and practice[END_REF] develop algorithms which are actually very close to proximal algorithms developed in convex optimization. They take ε > 0 in the data fidelity term [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] in order to make it differentiable.

During the last decade the community of convex optimization improved significantly the performances of algorithms for optimization of the sum of eventually non differentiable functions [START_REF] Chen | A proximal-based decomposition method for convex minimization problems[END_REF], [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF], [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF], [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF].

The proximal operator [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF] is one of the key tool in these works to derive convergent algorithms in the case where one of the function is differentiable. The forward-backward splitting algorithm [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] focus on the minimization of the sum of a convex differentiable function f • A + g with f satisfying a L-Lipschitz continuous gradient, A a linear operator and a convex eventually non smooth function g. This algorithm was used in the context of image deconvolution under a Poisson noise when a Anscombe's transform is applied in [START_REF] Dupé | A proximal iteration for deconvolving poisson noisy images using sparse representations[END_REF]. In [START_REF] Figueiredo | Restoration of poissonian images using alternating direction optimization[END_REF] the authors study the Alternative Direction Method of Multipliers (ADMM) described in the setting of Poisson noise, a method which doesn't require any differentiability of f or g. The algorithm PPXA which deals with the sum of possibly more than 2 functions (the functions may even be non differentiable) was applied in the setting of dynamical PET [START_REF] Pustelnik | Parallel algorithm and hybrid regularization for dynamic pet[END_REF], [START_REF] Pustelnik | Parallel proximal algorithm for image restoration using hybrid regularization[END_REF]. Mathematical proofs insure convergence but each iteration requires several steps. FISTA is an acceleration procedure for non constraint algorithms [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], and in the special case of the TV functional and ℓ 1 , additional constraints such as the positivity of the solution can be addressed [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF]. The choice of the parameter λ in the regularization term was studied in [START_REF] Bertero | A discrepancy principle for poisson data[END_REF] for data fidelity terms issued from a Poisson noise. In [START_REF] Chen | A proximal-based decomposition method for convex minimization problems[END_REF] an algorithm which belongs to primal-dual methods is described since primal and dual variables are successively computed in the algorithm. Then recently in [START_REF] Bonettini | An alternating extragradient method for total variation based image restoration from poisson data[END_REF] the authors developed an algorithm using also a primal dual scheme to miminize a data fidelity term based on Poisson noisy data penalized with TV. They apply it to denoising and deblurring problems in images. A new algorithm [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] was presented in order to minimize the sum of two convex functions without additional assumptions on their regularity but some others on the Legendre-Frechel transform of one of the function. It turns out that this setting fits perfectly with the goal of minimizing the -log likelihood penalized by non differentiable regularizations either in the case of the PET problem or in the one of the CT. In all these works theoretical convergence to the unique solution as well as numerical performances are adressed.

In the following we will present several algorithms to minimize the -log likelihood penalized by the TV and l 1 norm. We will compare their performance with well established algorithms.

Numerical algorithms

Let X and Y be two finite-dimensional real vector spaces embedded with an inner product ., . and the associated norm . = ., . .

Discrete setting

We take here the same notations as in [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. The image is a two dimensional vector of size N × N . We denote by X the Euclidean space IR N ×N , and Y = X × X. The space X will be endowed with the inner product u, v = 1≤i,j≤N u i,j v i,j and the norm u = u, u . We introduce a discrete version of the gradient operator. If u ∈ X, the gradient ∇u is a vector in Y given by: (∇u) i,j = ((∇u) 1 i,j , (∇u) 2 i,j ). with

(∇u) 1 i,j = u i+1,j -u i,j if i < N 0 if i = N (19) and (∇u) 2 i,j = u i,j+1 -u i,j if j < N 0 if j = N (20)
We also introduce a discrete version of the divergence operator. We define it by analogy with the continuous setting by div = -∇ * where ∇ * is the adjoint of ∇: that is, for every p ∈ Y and u ∈ X, (-divp, u) X = (p, ∇u) Y . It is easy to check that:

(div(p)) i,j =    p 1 i,j -p 1 i-1,j if 1 < i < N p 1 i,j if i=1 -p 1 i-1,j if i=N +    p 2 i,j -p 2 i,j-1 if 1 < j < N p 2 i,j if j=1 -p 2 i,j-1 if j=N (21)
We can define a discrete version of the Laplacian operator by setting ∆u = div(∇u) if u ∈ X.

Regularization choice

In this paper, we will consider the following regularization operators:

• The discrete total variation of u is then defined by:

J T V (u) = 1≤i,j≤N |(∇u) i,j | (22)
• A regularized version of the total variation:

J reg T V (u) = α 2 + |∇u| 2 , 1 = 1≤i,j≤N α 2 + |(∇u) i,j | 2 (23)
• A sparsity-inducing norm on a frame expansion (wavelets, curvelets, . . . ):

J l 1 ,ϕ (u) = λ∈Λ | u, ϕ λ | = R ϕ (u) l 1 (24)
where {ϕ λ } λ∈Λ is a tight frame, i.e. a family of elements of X such that

∀u ∈ X, u = λ∈Λ u, ϕ λ ϕ λ & u 2 = 1≤i,j≤N |u i,j | 2 = λ∈Λ | u, ϕ λ | 2 . ( 25
)
and R ϕ is the frame analysis operator:

u ∈ X → R ϕ (u) = { u, ϕ λ } λ∈Λ ∈ l 2 (R).
Frames are a generalization of orthornormal bases that include unions of several orthonormal bases and other redondant systems.

Sparsity-inducing norms on frames such as wavelet frames or curvelets allow to recover efficiently images that have sharp edges and thus are widely used in image processing.

4.3 Forward-backward splitting and FISTA

Proximal operator

We will use some mathematical concepts coming from convex analysis and refer the reader to [START_REF] Rockafellar | Convex Analysis[END_REF] for a complete introduction to this subject. Let F be a convex proper function. We recall that the subgradient of F , which is denoted by ∂F , is defined by

∂F (x) = {p ∈ X such that F (y) ≥ F (x) + p, y -x ∀y} (26)
One can show (see [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF] for more details, or [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF] for a complete proof) that for any h > 0 the following problem always has a unique solution :

min y hF (y) + 1 2 x -y 2 (27)
This solution is given by:

y = (I + h∂F ) -1 (x) = prox F h (x) (28)
The mapping (I + h∂F ) -1 is called the proximity operator. We refer to [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] for examples of proximal operator computations.

Forward-backward splitting

Assume one wants to solve :

min x F (x) + G(x) ( 29
)
where F is a convex C 1,1 function, with ∇F L Lipschitz, and G a simple convex function (simple means that the proximity operator of G is easy to compute). The Forward-Backward algorithm [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF][START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] reads in this case:

x 0 ∈ X x k+1 = (I + h∂G) -1 (x k -h∇F (x k )) ( 30
)
This algorithm is known to converge provided h ≤ 1/L. In terms of objective functions, the convergence speed is of order 1/k.

Acceleration

It has been shown by Nesterov in [START_REF] Nesterov | Introductory lectures on convex optimization: a basic course[END_REF][START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF][START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] , and by Beck and Teboule in [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], that the previous algorithm could be modified so that a convergence speed of order 1/k 2 is obtained.

The algorithm proposed by Beck and Teboule in [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] is the following:

         x 0 ∈ X ; y 1 = x 0 ; t 1 = 1; x k = (I + h∂G) -1 (y k -h∇F (y k )) t k+1 = 1+ √ 1+4t 2 k 2 y k+1 = x k + t k -1 t k+1 (x k -x k-1 ) (31)
This algorithm converges provided h ≤ 1/L.

Projected gradient algorithm

An interesting particular case for problem [START_REF] Herzen | Quantitative phase-contrast tomography of a liquid phantom using a conventional x-ray tube source[END_REF] arises when G is the indicator function of some closed convex set C, i.e.:

G(x) = 0 if x ∈ C +∞ otherwise (32)
Then Algorithm ( 30) is nothing but the classical projected gradient method, and Algorithm (31) an accelerated projected gradient method. In particular, (30) now reads:

x 0 ∈ X x k+1 = P C (x k -h∇F (x k )) (33)
with P C is the orthogonal projection operator onto C. (31) then reads:

         x 0 ∈ X y 1 = x 0 t 1 = 1 x k = P C (y k -h∇F (y k )) t k+1 = 1+ √ 1+4t 2 k 2 y k+1 = x k + t k -1 t k+1 (x k -x k-1 ) (34)
This algorithm converges provided h ≤ 1/L.

FISTA and constrained total variation

Beck and Teboule have shown in [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF] that FISTA could be used to solve the constrained total variation problem.

The considered problem is the following:

min u∈C J T V (u) + 1 2λ f -u 2 (35)
with C a closed non empty convex set.

In [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF], it is shown that the dual problem is a constrained smooth minimization problem. It can be solved with the previous projection algorithm with the FISTA acceleration.

Notice that when C = X (i.e. non constraint), then derivation of the dual problem was first done in [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF], and the acceleration using Nesterov ideas in [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in im[END_REF].

More precisely, the result shown in [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF] is the following:

Proposition 1 Let us set:

h(v) = -H C (f -λdivv) 2 + f -λdivv 2 (36)
where H C (u) = u -P C (u) and P C (u) is the orthogonal projection of u onto C. Let us define:

ṽ = arg min v ≤1 h(v) (37)
Then the solution of problem ( 35) is given by:

u = P C (f -λdivṽ) (38)
Hence it suffices to compute the minimizer of (37) to get the solution of Problem [START_REF] Lange | EM reconstruction algorithms for emission and transmission tomography[END_REF]. The minimization Problem [START_REF] Mariano-Goulart | A priori selection of the regularization parameters in emission tomography by fourier synthesis[END_REF] is the constrained minimization of a smooth function, and this can be done with the FISTA algorithm [START_REF] Huber | Robust Statistics[END_REF].

Indeed, it is proven in [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF] that h given by ( 36) satisfies: Lemma 1 The function h defined by ( 36) is continuously differentiable, and its gradient is given by ∇h

(v) = -2λ∇(P C (f -λdivv)) (39)

FISTA and constrained regularization

We show here how the results of Beck and Teboule in [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF] can be adapted to some general L 1 regularization. We consider the following problem (notice that when K is the gradient operator, then Problem [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF] is the constrained total variation regularization Problem [START_REF] Lange | EM reconstruction algorithms for emission and transmission tomography[END_REF]):

min u∈C Ku 1 + 1 2λ f -u 2 (40)
with C a closed non empty convex set. K is a continuous linear operator from X to Y (two finite-dimensional real vector spaces). Proposition 2 Let us set:

h K (v) = -H C (f + λK * v) 2 + f + λK * v 2 (41)
where

H C (u) = u -P C (u) and P C (u) is the orthogonal projection of u onto C. Let us define: ṽ = arg min v ≤1 h K (v) (42)
Then the solution of problem ( 40) is given by:

u = P C (f + λK * ṽ) (43)
Proof : We are interested in finding the minimizer of the following problem:

min u∈C 2λ Ku 1 + f -u 2 (44)
We start with the remark that the following equality holds:

Ku 1 = max |v|≤1 Ku, v = max |v|≤1 u, -K * v (45)
Hence Problem (44) can be written

min u∈C 2λ Ku 1 + f -u 2 = min u∈C max |v|≤1 2λ u, -K * v + f -u 2 (46)
We can swap here the min and the max (see e.g. [START_REF] Ekeland | Analyse convexe et problmes variationnels[END_REF]), and thus :

min u∈C 2λ Ku 1 + f -u 2 = max |v|≤1 min u∈C 2λ u, -K * v + f -u 2 (47)
We have:

2λ u, -K * v + f -u 2 = u -(f + λK * v) 2 -λ 2 K * v 2 -2λ f, K * v = u -(f + λK * v) 2 -f + λK * v 2 + f 2
Hence Problem (44) can be written:

max |v|≤1 min u∈C u -(f + λK * v) 2 -f + λK * v 2 (48)
The solution of the minimization with respect to u ∈ C is simply given by u = P C (f + λK * v) where P C is the orthogonal projection onto C.

Hence Problem (44) can be written:

max |v|≤1 P C (f + λK * v) -(f + λK * v) 2 -f + λK * v 2 (49)
And therefore Problem (44) amounts to solving

max |v|≤1 (-h K (v)) (50)
with

u = P C (f + λK * v).
Hence it suffices to compute the minimizer of (42) to get the solution of Problem [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF]. The minimization Problem [START_REF] Nicol | Design and construction of the clearpet/xpad small animal pet/ct scanner[END_REF] is the constrained minimization of a smooth function, and this can be done with the FISTA algorithm [START_REF] Huber | Robust Statistics[END_REF].

Indeed, one can show that the function h K given by (41) satisfies: Lemma 2 The function h K defined by ( 41) is continuously differentiable, and its gradient is given by ∇h

K (v) = -2λK(P C (f + λK * v)) (51)
Proof : Let us consider the function s(x) = H C (x) 2 = P C (x) -x 2 . It is standard (see e.g [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF]) that s ′ (x) = 2(x -

P C (x)). Now let s(v) = s(f + λK * v). We have ∇s(v) = λK(s ′ (f + λK * v)) = 2λK(f + λK * v -P C (f + λK * v)). Let r(v) = f + λK * v 2 . Then ∇r(v) = 2λK(f + λK * v). But h K (v) = s(v) -r(v).
And we get the result.

Chambolle Pock Algorithm 4.4.1 Introduction

Recently, A. Chambolle and T. Pock [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] have introduced an algorithm that will be useful to solve our minimization problem. The framework is the following. As before, X and Y are two finite-dimensional real vector spaces embedded with an inner product ., . and the associated norm . = ., . . We also introduce a continuous linear operator K : X → Y with respect to the induced norm

K = max { Kx | x ∈ X and x ≤ 1} (52)
Let F : X → [0, +∞) and G : X → [0, +∞) be two proper, convex, lower semi continuous (l.s.c.) functions.

We define the Legendre-Fenchel conjugate of F by :

F * (y) = max x∈X ( x, y -F (x)) (53)
The general minimization problem considered is the following saddle-point problem:

min x∈X max y∈Y ( Kx, y + G(x) -F * (y)) (54)
This saddle point problem is a primal-dual formulation for the nonlinear primal problem :

min x∈X (F (Kx) + G(x)) (55)
or of the corresponding dual problem max y∈Y (-G * (-K * y) + F * (y)) (56)

Algorithm

The algorithm proposed by A. Chambolle and T. Pock is the following:

• Initialization: Choose τ, σ > 0, (x 0 , y 0 ) ∈ X × Y ), and set x0 = x 0 .

• Iterations (n ≥ 0): Update x n , y n , xn as follows:

   y n+1 = (I + σ∂F * ) -1 (y n + σK xn ) x n+1 = (I + τ ∂G) -1 (x n -τ K * y n+1 ) xn+1 = 2x n+1 -x n (57)
The following theorem is proved in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]: Theorem 1 Let L = K , and assume Problem (54) has a saddle point. Choose τ σL 2 < 1, and let (x n , xn , y n ) be defined by (57). Then there exists a saddle point (x * , y * ) such that x n → x * and y n → y * .

Notice that Algorithm (54) can be used even when both F and G are non smooth. One should also be aware that similar approaches had been previously proposed in [START_REF] Chen | A proximal-based decomposition method for convex minimization problems[END_REF][START_REF] Briceo-Arias | Convex variational formulation with smooth coupling for multicomponent signal decomposition and recovery[END_REF].

Application to CBCT

Functional

We recall that the data term is the following :

L CT (µ) = y, Aµ + z exp(-Aµ), u (58)
y is the data, A the operator, µ the unknown, u the vector whose coordinates are all 1, and z is the number of photons.

Notice that we have to take into account the fact/constraint that µ ≥ 0. The considered minimization problem then reads:

min µ≥0 λJ(µ) + L CT (µ) ( 59
)
where J is some regularization operator.

Solving Problem (59) in the case of J reg T V

In this case, Problem (59) can be solved with the accelerated projected gradient algorithm [START_REF] Khoury | A geometrical calibration method for the pixscan micro-ct scanner[END_REF] with F (µ) = λJ(µ) + L CT (µ), and C = {µ ≥ 0}. It can easily be computed that:

∇J reg T V (µ) = -div ∇µ α 2 + |∇µ| 2 (60)
and

∇L CT (µ) = A * (y -z exp(-Aµ)) (61)
We denote by TVreg this algorithm.

Solving Problem (59) in the case of J l 1 ,ϕ

In this case, Problem (59) can be solved with the FISTA algorithm [START_REF] Huber | Robust Statistics[END_REF] with F = L CT in [START_REF] Herzen | Quantitative phase-contrast tomography of a liquid phantom using a conventional x-ray tube source[END_REF], and G(µ) = λJ l 1 ,ϕ (µ) + χ C (µ) We recall that J l 1 ,ϕ (µ) = λ | µ, ϕ λ |, and that χ C is the indicator function of C, i.e. χ C (µ) = 0 if µ ∈ C and +∞ otherwise. The proximal point prox λJ l 1 ,ϕ +χ C (µ) is computed with the fast algorithm of Section 4.3.6. We denote by FBwav this algorithm.

Solving Problem (59) in the case of J T V

In this case, Problem (59) can be solved with the FISTA algorithm [START_REF] Huber | Robust Statistics[END_REF] with F = L CT in [START_REF] Herzen | Quantitative phase-contrast tomography of a liquid phantom using a conventional x-ray tube source[END_REF], and G(µ) = λJ T V (µ) + χ C (µ). The proximal point prox J T V +χ C (µ) is computed with the algorithm of [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF] (see Section 4.3.5). We denote by FB-TV this algorithm.

6 Application to TEP

Functional

We recall that the data term is the following :

L T EP (v) = Bv -w log ǫ (Bv), u (62)
with the notation:

log ǫ (x) = log(x + ǫ) (63)
w is the data, B the operator, v the unknown and u the vector whose coordinates are all 1.

Notice that we have to take into account the fact/constraint that v ≥ 0. The considered minimization problem then reads:

min v≥0 λJ(v) + L T EP (v) (64)
6.2 Solving Problem (64) in the case of J reg T V

In this case, Problem (59) can be solved with the projected gradient algorithm [START_REF] Khodaverdi | Design study for the ClearPET/XPAD small animal PET/CT scanner[END_REF] with

F (v) = λJ(v) + L T EP (v), and C = {v ≥ 0}.
It can be easily computed that:

∇L T EP (v) = B * 1 - w Bv + ǫ (65)
We denote by TVreg this algorithm.

Lemma 3 Let F (x) = x -w log(x), u + χ C (x) = i x i -w i log(x i ) + χ IR + (x i ) (68)
Let y = (I + σ∂F * ) -1 (x). Then:

y i =      1 2 x i + 1 -(x i -1) 2 + 4σw i if w i > 0 1 if w i = 0 and x i ≥ 1 x i if w i = 0 and x i ≤ 1 (69)
Proof :

Let us set the function on IR:

f i (x i ) = x i -w i log(x i ) + χ IR + (x i ) ( 70
)
f i is a convex lower semi continuous function on IR (since

w i ≥ 0). Hence F (x) = i f i (x i
) is a convex lower semi continous on IR N . We have:

x = (I + σ∂F ) -1 (X) = arg min x x -X 2 2σ + F (x) = arg min x x -X 2 2σ + x -w log ǫ (x), u + χ C (x) = arg min x i 1 2σ (x i -X i ) 2 + f i (x i ) + χ IR + (x i )
with the obvious notations x = (x i ) and X = (X i ). Solving the previous minimization problem with respect to x amounts to solving each problem with respect to x i (since it is separable). We are thus led to consider the minimization problem: arg min

x i 1 2σ (x i -X i ) 2 + f i (x i ) + χ IR + (x i ) (71)
We now split the proof into two cases.

• First case : w i > 0

If x > 0, we have :

f ′ i (x) = 1 - w i x (72)
The solution of Problem (71) is thus given by:

xi = 1 2 X i -σ + (X i -σ) 2 + 4w i σ (73)
• Second case : w i = 0 Problem (71) then becomes: arg min

x i 1 2σ (x i -X i ) 2 + x i + χ IR + (x i ) (74)
The solution of Problem ( 71) is thus given by: xi = max(X i -σ, 0)

We deduce that if x = (I + σ∂F ) -1 (X), then xi is given by:

xi = 1 2 X i -σ + (X i -σ) 2 + 4w i σ if w i > 0 max(X i -σ, 0) if w i = 0 (76)
We use the identity:

(I + σ∂F * ) -1 (x) = x -σ I + 1 σ ∂F -1 x σ (77)
And if

y = I + 1 σ ∂F -1 x
σ then y i is given by:

y i = 1 2σ x i -1 + (x i -1) 2 + 4σw i if w i > 0 max( x i -1 σ , 0) if w i = 0 (78)
Thus we deduce the result.

Second approach

It is possible to avoid the inner loop of the previous approach by writing Problem (64) in another way. This is more elegant, but since it involves considering a higher dimensional operator, this will be less efficient in term of computation speed (as we will see in the numerical section of the paper). Problem (64) can be rewritten into:

min v F (Kv) + G(Bv) + χ C (v) (79) with K = ∇ for J T V and K = R ϕ for J l 1 ,ϕ and C = {v ≥ 0}, F (p) = λ p L 1 , G(x) = x -w log(x), u + χ C (x) for both J l 1 ,ϕ or J T V .
The associated saddle point problem is:

min v max y,z ( Kv, y + Bv, z -F * (y) -G * (z) + χ C (v)) (80)
The 3 optimality conditions are:

   0 ∈ K * y + A * z + ∂χ C (v) 0 ∈ Kv -∂F * (y) 0 ∈ Bv -∂G * (z) (81)
Chambolle-Pock algorithm thus becomes in this case (gradient descent in v, gradient ascent in y and z):

       v k+1 = (I + α∂χ C ) -1 (v k -αK * y k+1 -αB * z k+1 ) y k+1 = (I + β∂F * ) -1 (y k + βK vk ) z k+1 = (I + γ∂G * ) -1 (z k + γBv k ) vk+1 = 2v k+1 -v k (82)
There are 3 proximity operators to compute. It is straightforward to check that

(I + α∂χ C ) -1 (v) = P C (v) (83)
with P C the orthogonal projection onto C.

It is easy to see that

(I + β∂F * ) -1 (y) = P C (y) (84)
with P C the orthogonal projection onto C = {y / y ∞ ≤ 1}. As for (I + γ∂G * ) -1 , its value is given by Lemma 3 (see formula (69)). We denote by CP-TV this algorithm.

7 Results

Experimental setup

In this section, we present the results of our numerical experiments. We analyze the performances of the different algorithms we have presented so far, as well as three stateof-the-art algorithms. In the case of the TEP problem, we also ran experiments using the very similar spiral [START_REF] Harmany | This is SPIRAL-TAP: Sparse poisson intensity reconstruction algorithms theory and practice[END_REF] algorithm.

The experiments are made on simulated observations for both PET and CT. The simulations were done using the open-source Matlab toolbox Image Reconstruction Toolbox (irt) 1 for several levels of noise. The true objects that we seek to reconstruct are displayed in Figure 5 (these are 128x128 images for Zubal's phantom, 256x256 images otherwise). The three phantoms used here are standard tomography phantoms. Zubal's phantom displays tissues of different absorption and is used here to study the global performances of the algorithms. For this phantom, the performances are evaluated in terms of the standard Signal-To-Noise Ratio (snr), the Structural Similarity ssim [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF] (see Eq. ( 85) and ( 86)), as well as the computation time. The contrast and resolution phantoms are used to study the precision of the algorithms in terms of contrast or resolution i.e. to evaluate which is the lowest contrast and the smallest object size that can be resolved. As is well-known, global measures such as the snr are not accurate enough to evaluate fine contrast or resolution. Therefore, we also compute the better suited Contrast-to-Noise Ratio (cnr) [START_REF] Herzen | Quantitative phase-contrast tomography of a liquid phantom using a conventional x-ray tube source[END_REF] defined in Eq. ( 87).

When comparing a reconstructed image I to the true object T , the snr, ssim and cnr are defined as:

SN R(I, T ) =10 log 10 mean(I 2 ) mean(|I -T | 2 ) , (85)
) 87
For the ssim, w is a window sliding through the image and a and b are constants. For the cnr, the subscript in refers to the inside of the test tubes and out to the outside of the test tubes within the phantom. This measure is straightforward in the case of the resolution phantom since all test tubes have the same intensity. In the case of the contrast phantom, we average the cnr measures for the 6 test tubes since they have different contrasts.

Figure 5: The images we wish to reconstruct: (from left to right) Zubal's phantom in CBCT and in PET; the contrast phantom in CBCT and in TEP; the resolution phantom (both CBCT and PET).

CT

For the CT problem, the level of noise in the observations is driven by the source power measured by the photon counts. We started from a standard level (1e5 photons) and studied how the quality degrades at lower levels (1e4 and 1e3 photons).

We compared three algorithms we proposed, namely

• the gradient descent with a regularized version of the tv-norm (referred to as TVreg),

• the Forward-Backward algorithm with a l 1 penalization in Wavelet space using the Haar wavelet, 3 levels of decomposition (referred to as FB-Wav ),

• the Forward-Backward algorithm with a tv penalization (referred to as FB-TV ), to three state-of-the-art algorithms implemented in the irt toolbox:

• the filtered back-projection (referred to as FBP ),

• the MLEM algorithm (referred to as MLEM and described in [START_REF] Lange | EM reconstruction algorithms for emission and transmission tomography[END_REF]),

• the MLEM algorithm penalized by a Huber function (referred to as MLEM-Huber or MLEM-H and described in [START_REF] Erdoǧan | Monotonic algorithms for transmission tomography[END_REF]).

TEP

For the TEP problem, the level of noise in the observations is driven by the efficiency of the detector. We studied several levels of fcount: 5e5, 2e5, 1e5.

In addition to the six algorithms studied in the CT problem, we also ran expriments in TEP with three versions of the Chambolle-Pock algorithm we described in Section 4:

• Chambolle-Pock's algorithm computed using the first approach (see 6.4.1), with a tv penalization computed using FISTA (referred to as CP-TV-BT ),

• Chambolle-Pock's algorithm computed using the second approach (see 6.5), with a tv penalization (referred to as CP-TV ),

• Chambolle-Pock's algorithm computed using the second approach with a Ll 1 penalization in Wavelet space using the Haar wavelet, 3 levels of decomposition (referred to as CP-Wav ), as well as the spiral [START_REF] Harmany | This is SPIRAL-TAP: Sparse poisson intensity reconstruction algorithms theory and practice[END_REF] algorithm (referred to as SPIRAL).

In the next two subsections we give the results of our experiments.

CT results

Zubal's phantom

For this phantom, we show in Fig. 6 the reconstructions obtained when the level of noise increases from 1e4 to 1e2 photons count. Table 1 synthesizes the results. This table gives the mean snr, ssim and computation time for a number of realizations nb. repet of the noise indicated in the last column. If relevant, we also give the number of iterations nb. iter. and the regularization(s) parameter(s) (λ), which is tuned to yield the best snr. We indicate in bold the settings yielding the best snr and ssim.

Both snr and ssim indicate that the proximal algorithms yield the best results, and that they are more robust to noise than the state-of-the-art algorithms. The Filtered-Back-Projection clearly suffers many artifacts due to the-ill-posedness of the problem. Plain MLEM is able to recover good quality results at high photon count; however it is clear that it is also unstable as the noise increases. All other methods are able to stabilize at higher noise level since they incorporate a regularization. 1: CT reconstruction, photon count from z=10000 to z=100.

In terms of snr and ssim, MLEM-Huber yields comparable results to the proximal algorithms we proposed at high photon counts, however those indicators drop faster as noise increases for MLEM-Huber. Finally, the visual aspect of the forward-backward reconstructions is quite different from the MLEM-Huber ones. The flat part in Zubal's phantom tend to be reconstructed as flat parts with the forward-backward algorithms while they seem more "cloudy" with MLEM-Huber. This is precisely the advantage of using a L 1 -type regularization (as in the Forward-Backward algorithms) versus a quadratic one (as in MLEM-Huber). Notice that the blocky aspects in the wavelet-based Forward-Backward reconstructions is due to the fact that we use the Haar wavelet.

Effects of noise and number of projections on contrast and resolution

We now study how the quality degrades with noise and when the number of angles of projections decreases for both the resolution and contrast phantom. The noise level considered are again 1e4, 1e3 and 1e2 photons count, while the number of angle of projections is 90, 60 and 30. The performances are evaluated in terms of cnr, ssim and snr. But for each algorithm, the hyperparameters are tuned to yield the best cnr. Since the algorithms we proposed yields similar results and Filtered-Back-Projection is clearly not efficient, we summarize in Tables 2 and3 the numerical results we obtained for the best proximal algorithm (FB-TV ) we proposed and two state-of-the-art algorithms only. Images of reconstructions are displayed in Figures 7 and9. Furthermore, profiles are extracted from both phantoms and their reconstructions to assess the precision in terms of resolution and contrast recovery for each method. These are shown in Fig. 8 For the proximal approaches and the MLEM approaches we also noticed that the quality of the reconstructions does not degrade much when the number of projections decays from 90 to 30. On the other hand, the quality does decay with the photon count, as may be seen in Figure 7 and Figure 9. These figures show the reconstructions for the two MLEM approaches and the proximal FB-TV approach, for 60 projections as the photon count decreases from 1e4 to 1e2. The snr, ssim and cnr of each reconstruction are given below it. These three indicators all show that FB-TV performs the best, then MLEM-Huber and finally MLEM for a photon count and number of projection fixed. One may also notice that the indicators degrade faster for MLEM, a little slower for MLEM-Huber and yet even slower for FB-TV as the photon count decreases.

As we saw for Zubal's phantom , the visual characteristics of the MLEM-Huber and FB-TV are quite different, the former yielding a "cloudy" reconstructions, the latter flat-by-part ones. This visual sensation at low photon count is also verified at higher counts. To show this in more details, we plot in Figure 8 a profile cut through the test tubes for both the contrast and resolution reconstructions for both algorithms.

As far as the contrast phantom is concerned, there is a clear advantage to using the FB-TV algorithm over the MLEM-Huber, since the latter will not reconstruct the two lowest contrasted tubes from 1000 photons count to counts below, while the FB-TV algorithm reconstructs all contrasts at 1000 photons and only fails at reconstructing the faintest one at 100 photons. As far as the resolution phantom is concerned, the differences are not as clear to the eye, since one may distinguish objects of the same size for both approaches. However one may argue that since the FB-TV are flat by parts while the ones obtained with MLEM-Huber are fuzzier, when feeding such reconstructions to automatic programs that segment different parts, those will yield better performances when using the FB-TV reconstructions than the MLEM-Huber ones since they rely on edge detectors.

TEP results

Zubal's phantom

As in the CT case, for Zubal's phantom, we study the performances of the different algorithms in terms of snr and ssim (the hyperparameters are tuned to yield the best snr). Table 4 synthesizes the results (see Section 7.2.1 for the legend). We show in Fig. 10 the reconstructions obtained for the levels of noise fcount = 5e5 and 1e5 (note that for the algorithms we proposed, the figure displays only the most efficient ones: Chambolle-Pock's algorithm).

The conclusions here are very similar to the CT case, the algorithm we proposed yielding constantly better snr and ssim, and being more robust to increasing noise levels. MLEM-Huber does yield interesting results though not as high quality as the proximal algorithms. The "cloudy" and "flat-by-part" aspects again differentiate the L 1 versus L 2 regularizations.

Effects of noise and number of projections on contrast and resolution

As in the CT case, we study how the quality degrades with noise and when the number of angles of projections decreases for both the resolution and contrast phantom. The noise level considered are again fcount = 5e5, 2e5 and 1e5, while the number of angle of projections is 90, 60 and 30. The performances evaluated in terms of snr and ssim displayed in Tables 5 and6.

Compared to the CT problem, the PET problem seems to be a much harder problem to solve. The quality of the reconstructions for all algorithms compared to the CT ones is much degraded. One notes however that as happens for CT, the Filtered Back-Projection is fast left behind by the other methods. Again all proximal algorithms perform similarly, although the tv based algorithms give slightly better results than the wavelet-based ones. We focus on the Chambolle-Pock algorithm, which deals with the exact fit-to-data term (ǫ = 0) as opposed to the Forward-Backward ones.

In general, the influence of the noise and number of projections is more pronounced in PET than in CT and this time, the number of projections has more influence than the level of noise. Therefore, we rather present images of reconstructions at a given level of noise (f count = 2e5) when the number of projections decreases in Figures 11 and12.

Although MLEM-Huber yields better numerical numerical results for low noise ( f count = 5e5), at low fcounts (2e5 and 1e5) i.e. when the dose decreases, all three indicators show that the proximal tv-based algorithm CP-TV-BT performs better than MLEM-Huber, which itself performs better than a plain MLEM. The plain MLEM also degrades faster than the regularized algorithms (see Figures 11 and12). Here however, one notices that the differences between the indicators is not as marked as for CT.

Concerning the contrast phantom, visual inspection of Figure 11 seems to show that the lower right tube has more chances to be recovered by the MLEM-Huber while it is the opposite for the top right tube. Examining the profile in Figure 13, the explanation is that the CP-TV-BT slightly underestimates the intensity of the flat part, which is a well known artifact of L 1 estimation. The results obtained with the resolution phantom are similarly not totally conclusive since it is not clear whether one method or the other fails at recovering the finest resolution objects. However in both cases, it remains clear that CP-TV-BT will lead to reconstructions that are flat by part, thus making it more likely to yield high quality automatic segmentations of the different imaged tissues than the MLEM-Huber approach.

Results on real data

In this section, we present the results obtained on real CBCT acquisitions obtained with the PIXSCAN II demonstrator developed at CPPM and incorporating the XPAD3 hybrid pixels detectors. As described in Section 2.1.2, we used a wide aperture cone beam, Tungsten target, X-ray tube (90kV, 60W) whose spectrum was hardened by a 2.5mm Aluminium filter. In this configuration, four complete acquisitions have been realized where 15000, 10000, 1000 and 600 photons were counted on average per pixel for flat irradiation conditions (in absence of a mouse) after respective 1.5s, 1s, 100ms and 50ms exposure times. For each acquisition, 360 projections of the same healthy mouse were acquired every degree in about 15 minutes per acquisition.

Image reconstruction was then performed using the different algorithms described in 7.4.1. For each of the four acquisitions, we performed the 2D reconstruction of a slice, chosen to be the medium slice lying in the plane containing the optical axis, which joins the source point and the detector and is perpendicular to the detector plane. This slice corresponds to the abdominal region of the mouse where soft tissues and spine bones are visible. Given the pixel size of 130 × 130µm 2 and a geometrical zoom factor of two, a slice of 614 × 614 pixels of an average resolution of 65 × 65µm 2 and of thickness 65µm was reconstructed.

Experimental setup

The slice reconstruction was performed by the Forward-Backward algorithm with a tv penalization (referred to as FB-TV ). The results have been compared with a standard filtered back-projection implemented in the IRT toolbox. In particular, the system matrix was computed with a strip-integral model incorporating the geometry of the acquisition obtained as described in [START_REF] Khoury | A geometrical calibration method for the pixscan micro-ct scanner[END_REF][START_REF] Nicol | Design and construction of the clearpet/xpad small animal pet/ct scanner[END_REF].

We performed three different experiments for each acquisition: the reconstruction has been performed with 90, 60 and 36 projections obtained by extracting data respectively every 4, 6 and 10 degrees from the complete acquisition with 360 projections. These experiments aim at evaluating the robustness of our algorithms to the reduction of projection angles, one of the current hot topic in the field of CBCT reconstruction.

Results and discussion

Results obtained with a decimated dataset of 60 projections over the 360 that have been acquired are displayed on Figures 14 to 17. They respectively correspond to photon counts of 15000, 10000, 1000 and 600 in flat irradiation conditions. Note that between these two extreme experimental conditions, the X-ray dose absorbed by the imaged mouse is reduced by a factor 25 for a full acquisition of 360 projection angles. For each experiment, we compare the results obtained with four different values of the parameter λ, i.e. λ ∈ {5, 15, 25, 40} with the results of the Filtered Back-Projection (FBP) reconstructed from the full acquired data (i.e. 360 projections) and from the decimated acquired data (i.e. 60 projections). The FBP reconstruction from the full data can be considered as the reference reconstruction for these tests, although it suffers from artifacts introduced by the erratic behavior of a very small proportion of pixels that can not be fully corrected by pre-processing. We can typically notice circular-shaped artifacts in the mouse body onto the reconstructed slices with FBP at 15000 and 10000 photon counts, and several rays of light are visible onto the reconstructed slices for number of angles lower than 360 whatever the photon count. Moreover, the Filtered Back-Projection is known to produce negative values without any physical meaning. For this reason, it does not appear relevant to use a snr or a ssim criteria using the FBP reconstruction with full data as a reference.

The very good quality of results obtained for usual flat irradiation (10000 photon per pixels) highlights the robustness of the FB-TV algorithm with respect to the number of projection angles. With a strong enough regularization weight, the reconstruction from only 60 projection angles is not subject to any initial artifacts induced by some dead or failing pixels : ringing artifacts are eliminated as well as spurious rays of light. Moreover, the essential structures of the image are preserved, which allows to identify the main biological features : pawns (outside the main body), spine bone (bright part at the center of the body), a very tiny part of a rib bone (small bright part around the bottom of the body), some air or gas areas (black holes at the top and top-left of the body), and the cylindrical plastic support of the mouse (on the top-left of the images, outside the body). The Total Variation prior seems well adapted to the presented images, but these first results on real data have to be completed by additional tests onto images from other anatomical parts that can be significantly more complex and less adapted to the TV prior.

An optimal choice of the regularization parameter is obviously highly subjective on real data. A value of 25 or 40 for the λ parameter strongly regularizes the images whatever the photon count whereas a value of 5 or 15 will mainly reduce the noise and preserve the finest details that might be of interest from a medical point of view. The images acquired with a 15000 and 10000 photon count and reconstructed with λ ∈ {25, 40} are accurate and edges are sharp. However, the same values of λ for lower photon counts introduce a blur which in particular slightly degrades the accuracy of the paws or the spine of the mouse. A lower value of λ leads to reconstructions that are satisfactory for a diagnosis, but that may be noisy and then less adapted to be the input of a segmentation or classification algorithm, or to a statistical analysis. The trade-off has then to be made jointly with medical experts, and it strongly depends on the future use of the images.

Figure 18 displays a comparison of reconstructed images as a function of the number of projection angles used as measurements. Reconstructions with 90, 60 and 36 projection angles are showed together with the reference reconstruction obtained from the full set of 360 projections by Filtered Back-Projection. These reconstructions have been obtained with the same value λ = 20. Even if the quality of images is logically decreasing when the number of projections is reduced, the overall quality of reconstruction with only 36 projection angles and a 600 photon count is high enough to visually identify the main features. This worst tested case corresponds to a theoretical reduction of X-ray dose of 250 when compared to a reconstruction with 360 projection angles with a photon count of 15000. However, these images suffer from noise, due to the low level of the acquired signal which is consequently strongly corrupted by Poisson noise. This noise can only be completely eliminated by a strong regularization which may introduce a blur and make a diagnosis more difficult to establish.

Figure 1 :

 1 Figure 1: Several mass attenuation coefficients curves as a function of energy. Notice that the absorption coefficient is decreasing with the energy excepted for Iodine and Barium, that are characterized by a strong discontinuity due to their intrinsic K-edge energy.

Figure 2 :

 2 Figure 2: Sketch of a classical Cone-Beam Computerized Tomography (CBCT) Scan. The X-ray source delivers a cone-beam and rotates around the patient lying on a bed. The X-ray detectors are located on a ring all around the patient.

Figure 3 :

 3 Figure 3: Picture of the PIXSCAN II scanner with the three distinct blocks. From left to right : the XPAD3 detector block, the imaged object block, and the source block.

Figure 4 :

 4 Figure 4: Left: sketch of the ClearPET/XPAD. One can see the partial ring detectors (in green) dedicated to the PET scan, and the X-ray source and the XPAD3 camera for the CT scan that form a line and fill the gap in the partial ring of PET detectors. Right: picture of the ClearPET/XPAD.

SSIM

  (I, T ) =mean w (2mean(I w )mean(T w ) + a)(2cov(I w , T w ) + b) mean(I w ) 2 + mean(T 2 w + a)(var(I w ) + var(T w ) + b) , (86) CN R(I) = |mean(I in ) -mean(I out )| var(I in ) + var(I out ) . (

Figure 6 :

 6 Figure 6: CT reconstructions for Zubal's phantom, two top rows: photon count z=10000, two bottom rows: photon count z=100.

Figure 7 :

 7 Figure 7: CT contrast phantom reconstructions for 60 projections. Left to right: MLEM, MLEM-Huber and FB-TV algorithms. Top to bottom: 1e4 to 1e2 photon count.

Figure 8 :

 8 Figure 8: Profile of the CBCT contrast (left) and resolution (right) phantom reconstruction results for 60 angles of projections a noise level of 1e3 photons.

Figure 9 :

 9 Figure 9: CT resolution phantom reconstructions for 60 projections. Left to right: MLEM, and FB-TV algorithms. Top to bottom: 1e4 to 1e2 photon count.

Figure 11 :

 11 Figure 11: PET contrast phantom reconstruction results for a noise level of f count = 2e5. Leftmost to rightmost columns: MLEM, MLEM-Huber and CP-TV-BT algorithms. Top to bottom row: from 90 to 30 projections.

Figure 12 :

 12 Figure 12: PET resolution phantom reconstruction results for 60 angles of projections and a noise level of fcount = 2e5. Leftmost to rightmost columns: MLEM, MLEM-Huber and CP-TV-BT algorithms.

Figure 13 :

 13 Figure 13: Profile of the PET contrast (left) and resolution (right) phantom reconstruction results for 60 angles of projections and a noise level of fcount = 2e5.

Figure 14 :

 14 Figure 14: Results obtained on real CBCT acquisitions with flat irradiation of 15000 photons per pixel. Top-left: Filtered back-projection with 360 projection angles (full data), top-middle: Filtered back-projection with 60 projection angles (decimated data). Top-right: FB-TV reconstruction with 60 projection angles and λ = 5. Bottom, from left to right: FB-TV reconstruction with 60 projection angles with λ = 15, λ = 25 and λ = 40.

Figure 15 :

 15 Figure 15: Results obtained on real CBCT acquisitions with flat irradiation of 10000 photons per pixel. Top-left: Filtered back-projection with 360 projection angles (full data), top-middle: Filtered back-projection with 60 projection angles (decimated data). Top-right: FB-TV reconstruction with 60 projection angles and λ = 5. Bottom, from left to right: FB-TV reconstruction with 60 projection angles with λ = 15, λ = 25 and λ = 40.

Figure 16 :

 16 Figure 16: Results obtained on real CBCT acquisitions with flat irradiation of 1000 photons per pixel. Top-left: Filtered back-projection with 360 projection angles (full data), top-middle: Filtered back-projection with 60 projection angles (decimated data). Top-right: FB-TV reconstruction with 60 projection angles and λ = 5. Bottom, from left to right: FB-TV reconstruction with 60 projection angles with λ = 15, λ = 25 and λ = 40.

Figure 17 :

 17 Figure 17: Results obtained on real CBCT acquisitions with flat irradiation of 600 photons per pixel. Top-left: Filtered back-projection with 360 projection angles (full data), top-middle: Filtered back-projection with 60 projection angles (decimated data). Top-right: FB-TV reconstruction with 60 projection angles and λ = 5. Bottom, from left to right: FB-TV reconstruction with 60 projection angles with λ = 15, λ = 25 and λ = 40.

Figure 18 :

 18 Figure 18: Results obtained on real CBCT acquisitions. From left to right : Filtered Back-Projection reconstruction from the full data set (360 projection angles), FB-TV reconstruction with 90, 60 and 36 projection angles with λ = 20. From top to bottom : Reconstructions performed from real data obtained with a photon count of 15000, 10000, 1000 and 600 in flat irradiation conditions.

Table

	Photon	Algorithm snr	ssim	λ	nb. iter. time (s) nb. repet.
	count
		TVreg	19.57 0.900	300	300	44	100
		FB-Wav	18.63 0.938	75	300	110	100
	1e4	FB-TV FBP	19.61 0.916 13.40 0.395	300 -	300 -	85 0.12	100 25
		MLEM	15.60 0.711	-	120	46	25
		MLEM-H 19.61 0.856	2e5	1000		25
		TVreg	15.06 0.808	200	300	36	100
		FB-Wav	14.06 0.826	25	300	110	100
	1e3	FB-TV FBP	15.10 0.845 200 9.08 0.201 -	300 -	85 0.09	100 25
		MLEM	11.86 0.462	-	43	14	25
		MLEM-H 14.52 0.680	7e5	752		25
		TVreg	11.34 0.625	80	300	32	100
		FB-Wav	10.62 0.695	10	300	110	100
	1e2	FB-TV FBP	11.35 0.690 0.44 0.076	80 -	300 -	78 0.07	100 25
		MLEM	7.90	0.200	-	17	5.67	25
		MLEM-H 10.78 0.489 3.5e4	605		25

Table 2 :

 2 . CT contrast phantom reconstruction results.

		Nb. angles	90		60		30
	Photon	Algorithm cnr ssim	snr	cnr ssim	snr	cnr ssim	snr
	count
		TVreg	5.18 0.866 21.86 4.69 0.867 21.04 4.08 0.855 19.54
		FB-Wav	4.91 0.911 21.53 4.18 0.911 20.09 3.57 0.884 18.28
	1e4	FB-TV FBP	6.05 0.944 23.03 5.33 0.914 21.57 4.74 0.912 20.78 2.23 0.201 8.98 1.84 0.155 8.33 1.35 0.100 6.79
		MLEM	2.88 0.425 15.13 2.69 0.393 14.71 2.34 0.353 13.39
		MLEM-H 3.33 0.492 17.52 3.23 0.521 17.50 3.00 0.591 17.13
		TVreg	3.17 0.793 17.15 3.00 0.743 16.95 2.78 0.728 15.46
		FB-Wav	3.18 0.851 17.68 2.96 0.839 17.01 2.61 0.802 15.39
	1e3	FB-TV FBP	3.93 0.831 17.65 3.50 0.817 17.42 3.36 0.756 15.25 0.82 0.046 3.57 0.65 0.033 2.09 0.44 0.017 -0.79
		MLEM	1.95 0.274 11.87 1.77 0.253 11.36 1.53 0.218 10.33
		MLEM-H 2.27 0.337 13.49 2.11 0.351 13.33 2.04 0.393 13.47
		TVreg	2.17 0.686 13.30 2.10 0.699 13.28 1.80 0.679 12.21
		FB-Wav	2.03 0.728 13.83 2.08 0.779 12.93 1.65 0.736 11.48
	1e2	FB-TV FBP	2.30 0.706 12.90 2.34 0.779 13.90 2.07 0.704 12.30 0.25 0.005 -5.46 0.21 0.003 -7.32 0.14 0.001 -10.28
		MLEM	1.17 0.187	8.77	1.06 0.170	8.13	0.87 0.164	7.36
		MLEM-H 1.53 0.282 10.78 1.46 0.294 10.28 1.36 0.334 10.07

Table 3 :

 3 CT resolution phantom reconstruction results.

		Nb. angles	90		60		30
	Photon	Algorithm cnr ssim	snr	cnr ssim	snr	cnr ssim	snr
	count
		TVreg	6.11 0.928 21.21 5.50 0.912 20.35 3.00 0.719 15.60
		FB-Wav	0.12 0.176	1.20	4.74 0.875 19.15 2.76 0.721 15.19
	1e4	FB-TV FBP	6.27 0.938 21.38 5.66 0.925 20.54 3.03 0.722 15.67 3.31 0.349 8.26 2.73 0.300 7.69 1.78 0.219 6.17
		MLEM	3.64 0.556 15.27 3.27 0.526 14.60 2.19 0.402 12.04
		MLEM-H 5.53 0.852 20.32 4.69 0.818 18.85 2.45 0.623 14.27
		TVreg	3.37 0.811 16.65 2.95 0.753 15.75 1.97 0.516 12.50
		FB-Wav	2.96 0.779 15.63 2.64 0.586 14.56 1.89 0.495 12.36
	1e3	FB-TV FBP	3.41 0.825 16.74 2.98 0.766 15.84 1.97 0.520 12.53 2.05 0.218 6.91 1.68 0.182 5.92 1.13 0.119 3.84
		MLEM	2.36 0.355 12.59 2.07 0.322 11.88 1.48 0.245 10.20
		MLEM-H 3.08 0.683 15.79 2.50 0.615 14.26 1.68 0.476 11.85
		TVreg	1.71 0.521 12.22 1.48 0.505 11.56 1.11 0.484 10.49
		FB-Wav	1.54 0.437 11.50 1.37 0.490 11.27 0.99 0.459 10.13
	1e2	FB-TV FBP	1.70 0.526 12.17 1.48 0.508 11.52 1.08 0.490 10.43 0.75 0.063 0.99 0.62 0.045 -0.50 0.42 0.019 -3.36
		MLEM	1.14 0.170	9.38	0.96 0.149	8.94	0.67 0.113	8.11
		MLEM-H 1.57 0.451 11.46 1.30 0.418 10.71 0.99 0.364	9.55

Table 5 :

 5 TEP contrast phantom reconstruction results.

	Photon	Algorithm	snr	ssim	λ	nb. iter. time (s) nb. repet.
	count
		TVreg	15.33 0.902	0.70	200	10	15
		FB-Wav	14.77 0.889	0.10	150	89	25
		FB-TV	15.37 0.905	0.70	100	62	15
		CP-Wav	14.68 0.885	0.10	80	63	25
	5e5	CP-TV-BT 15.32 0.905 CP-TV 14.84 0.860	0.70 0.70	80 400	63 266	15 15
		SPIRAL	15.17 0.905	0.70	100	76	10
		FBP	11.59 0.429	-	-	0.04	25
		MLEM	13.38 0.819	-	17	2	25
		MLEM-H	15.22 0.866 0.9/0.25	267	46	25
		TVreg	13.68 0.866	0.50	200	11	100
		FBwav	13.11 0.867	0.0875	150	89	50
		FB-TV	13.71 0.875	0.50	100	65	100
		CPwav	13.46 0.870	0.0875	100	79	50
	2e5	CP-TV-BT 14.17 0.882 CP-TV 14.01 0.867	0.50 0.50	100 150	82 107	100 100
		SPIRAL	13.40 0.872	0.50	100	65	100
		FBP	9.03	0.322	-	-	0.07	25
		MLEM	12.08 0.774	-	12	2	25
		MLEM-H	13.97 0.853 0.9./0.25	274	53	25
		TVreg	12.12 0.841	0.40	200	13	10
		FBwav	11.55 0.834	0.0625	150	89	50
		FB-TV	12.14 0.847	0.40	100	68	10
		CPwav	11.65 0.835	0.0625	50	40	50
	1e5	CP-TV-BT 13.13 0.862 CP-TV 12.86 0.823	0.40 0.40	50 100	46 78	10 10
		SPIRAL	11.77 0.841	0.40	100	86	10
		FBP	6.66	0.254	-	-	0.08	25
		MLEM	11.06 0.731	-	10	2	25
		MLEM-H	12.92 0.837 0.8/0.25	278	58	25
	Table 4: TEP Zubal phantom reconstruction, detector efficiency fcount=500000 to
	100000.

Table 6 :

 6 TEP resolution phantom reconstruction results.

http://www.eecs.umich.edu/ fessler/irt/irt

In this case, Problem (64) can be solved with the classical Forward-Backward splitting algorithm [START_REF] Hudson | Accelerated image reconstruction using ordered subsets of projection data[END_REF] with F = L T EP in [START_REF] Herzen | Quantitative phase-contrast tomography of a liquid phantom using a conventional x-ray tube source[END_REF], and In this case, Problem (64) can be solved with the classical Forward-Backward splitting algorithm [START_REF] Hudson | Accelerated image reconstruction using ordered subsets of projection data[END_REF] with F = L T EP in [START_REF] Herzen | Quantitative phase-contrast tomography of a liquid phantom using a conventional x-ray tube source[END_REF], and

is computed with the algorithm of [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF] (see Section 4.3.5).

Solving Problem (64) with Chambolle-Pock algorithm (57)

In this case, Problem (64) can be solved with Chambolle-Pock algorithm (57) and even with ǫ = 0

First approach

We have to work a little to adapt the framework to our case. Problem (64) can be rewritten into:

with C = {v ≥ 0}, F (x) = x-w log ǫ (x), u +χ C (x), and G(x) = λJ(x)+χ C (x). (recall that u = (1, .., 1)). Then Chambolle-Pock algorithm becomes in this case:

(I + σ∂F * ) -1 is given by (69) (see Lemma 3 here-after).

In the case when wavelet regularization is considered, i.e. J = J l 1 ,ϕ , then (I +τ ∂G) -1 is computed with the fast algorithm of Section 4.3.6.

We denote by CPwav this algorithm.

In the case of total variation regularization, i.e. J = J T V , then (I + τ ∂G) -1 is computed with the algorithm of [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF] (see Section 4.3.5).

We denote by CP-TV-BT this algorithm.