Brest Latim U1101

Dominique Thers

A collaboration between

i

Declaration of Authorship

I, Debora Giovagnoli, declare that this thesis titled, "Image Reconstruction for Three-gamma PET Imaging" and the work presented in it are my own. I confirm that:

• This work was done wholly or mainly while in candidature for a research degree at this University.

• Where any part of this thesis has previously been submitted for a degree or any other qualification at this University or any other institution, this has been clearly stated.

• Where I have consulted the published work of others, this is always clearly attributed.

• Where I have quoted from the work of others, the source is always given.

With the exception of such quotations, this thesis is entirely my own work.

• I have acknowledged all main sources of help.

• Where the thesis is based on work done by myself jointly with others,

PSF plot and

Gaussian fit for the bone sphere in (0, 0, -60). . .

PSF plot and

Gaussian fit for the water sphere in (0, 0, 0). . . .

PSF plot and

Gaussian fit for the lung sphere in (0, 0, 60). . . . in numerical order sphere of 2-mm, 4-mm and 8-mm radius.

With "conv" we refer to variant kernel convolution, while "stat conv" refers to the invariant kernel convolution correction. . . L'algorithme LM-MLEM pseudo-TOF peut être exprimé comme suit:

l (k+1) j = l (k) j  i,b p i,j,b  i p i,j,b  j 0 p i,j 0 ,b l (k) j 0 . ( 1
)
où i fait référence au bin de la LOR et b au bin pseudo-temporel correspondant à l'événement considéré; l (k) j est la valeur de l'image dans le voxel j pour l'estimation k-th. On désigne comme pseudo-time-bin le bin du temps correspondant à l'équivalent en ps de la position LCI.

Resume

Nous montrons des résultats liés au fantôme de type NEMA obtenu par reconstruction de 3-g et 2-g MLEM. Pour la reconstruction pseudo-TOF, nous avons utilisé deux approches pseudo-TOF différentes: dans une reconstruction nous avons utilisé un écart-type gaussien fixe de 70 ps, qui est la valeur correspondant à l'incertitude moyenne sur le calcul LCI. Dans la seconde approche, la résolution était basée sur les événements, donc variable, calculée pour chaque événement comme l'équivalent en temps de la différence entre le LCI et la projection du point d'annihilation sur la LOR, comme le montre la Les travaux présentés dans cette thèse ont porté sur le développement d'une nouvelle technique de reconstruction pour l'imagerie 3-gammas, qui a été évaluée sur la géométrie du système d'imagerie XEMIS2. L'étude XEMIS2 est la partie préclinique du projet plus général XEMIS, qui a déjà développé un premier prototype, le XEMIS1. La caractérisation de la caméra Compton XEMIS1 a donné la démonstration expérimentale de la faisabilité de l'imagerie trois gammas, on retrouve les études précédentes dans Gallego-Manzano [START_REF] Gallego | Optimization of a single-phase liquid xenon Compton camera for 3g medical imaging[END_REF],

Oger [START_REF] Oger | Développement expérimental d'un télescope Compton au xenon liquide pour l'imagerie médicale fonctionnelle[END_REF] et Grignon [START_REF] Grignon | Étude et développement d'un télescope Compton au xénon liquide dédié à l'imagerie médicale fonctionnelle[END_REF].

Dans ce travail, nous utilisons les interactions du troisième-gamma dans le LXe pour déterminer un cône Compton, dont le point d'intersection avec le LOR, obtenu à partir des photons de coïncidence, permet de restreindre la localisation des LOR section impliquée dans l'événement combiné annihilation + émission. L'événement est ensuite reconstruit à l'aide d'un pseudo-TOF/MLEM algorithme où l'intersection ligne-cône est prise comme centre d'un PDF gaussien utilisé dans la reconstruction, de la même manière que dans TOF-PET, d'où la dénomination de la méthode proposée comme pseudo-TOF. Un avantage principal associé à l'approche proposée est de permettre l'intégration de l'information du troisième gamma dans des cadres de reconstruction existants avec une capacité TOF, tels que CASTOR.

Resume

Introduction

Clinical nuclear medicine can be considered to have started as early as 1927 when Blumgart and Weiss [START_REF] Herrmann L Blumgart | STUDIES ON THE VELOC-ITY OF BLOOD FLOW: IV. The Velocity of Blood Flow and Its Relation to Other Aspects of the Circulation in Patients with Arteriosclerosis and in Patients with Arterial Hypertension 1[END_REF] analyzed the normal and abnormal circulation of blood in patients using the radioisotope Bi 214 . The first medical detectors were mainly borrowed from experimental and nuclear physics, made of scintillators coupled to photomultipliers (PMTs) through a light guide to read the signal.

In the 1970s computed tomography (CT) scanning systems were invented in conjunction with tomographic reconstruction. In the same years magnetic resonance imaging (MRI) was developed, and later on positron emission tomography (PET) and single photon emission tomography (SPECT) imaging systems [START_REF] James | From Roentgen to magnetic resonance imaging: the history of medical imaging[END_REF].

Nowadays, nuclear medicine allows for the study of in vivo biological and metabolic processes in the body using mainly two functional molecular imaging modalities: SPECT and PET imaging. Both imaging modalities are based on the use of radiopharmaceuticals (or radiotracer), whose biodistribution in the body can be evaluated using external detection and tomographic reconstruction. The radiotracer is composed of a substance that follows a physiological process labeled with a radioactive isotope that emits radiation detectable from outside the body. For SPECT diagnostics the radiopharmaceutical is a single photon emitter, detected through a gamma camera rotating around the patient, while in PET studies positron emitters are employed and the detection occurs in a ring-shaped detector surrounding the body.

While nuclear imaging allows for improved screening, diagnosis and follow-up of several diseases, it also introduces some risks related to radiation exposure [START_REF] Stauss | Guidelines for 18 F-FDG PET and PET-CT imaging in paediatric oncology[END_REF]. Nowadays the two main factors that guide research and development in nuclear imaging are arguably the reduction of the administered dose and of the acquisition time, for both safety and practical purposes.

These two factors have a critical influence on image quality as the signal to noise ratio (SNR) is mostly dependent on the number of photon counts in the detectors [START_REF] Brix | Radiation exposure of patients undergoing wholebody dual-modality 18F-FDG PET/CT examinations[END_REF][START_REF] Xu | 200x Low-dose PET Reconstruction using Deep Learning[END_REF][START_REF] Baxter | Determination of signal-to-noise ratios and spectral SNRs in cryo-EM low-dose imaging of molecules[END_REF]. To maintain sufficient image quality standards required in modern medicine while decreasing the total activity, image processing and software-based methods have been implemented [START_REF] Shtok | Sparsity-based sinogram denoising for low-dose computed tomography[END_REF][START_REF] Chen | Artifact Suppressed Dictionary Learning for Low-Dose CT Image Processing[END_REF], together with the use of new imaging systems [START_REF] Romo-Luque | PETALO: Time-of-Flight PET with liquid xenon[END_REF][START_REF] Salvador | Design of a high performances small animal PET system with axial oriented crystals and DOI capability[END_REF][START_REF] Aggarwal | Impact of preclinical PET scanner characteristics on the overall image quality[END_REF] New systems, known as time of flight (TOF)-PET, were first conceived in the 1970s. Their goal is to measure the time difference between the two Algorithms used to reconstruct tomographic images can be classified into two main groups: analytical and iterative methods. In the first case, we deal with an inversion problem between a function and its line integral. In the second case we model all the possible aspects of data acquisition and all the physical interactions that can occur and we look for the image that best fits the data in an iterative fashion. The data on which we base our reconstruction are acquired detecting 511-keV coincidence photons in PET and single photons in SPECT. In both cases, the limited amount of events and the particle stochastic behavior result in noisy outputs. Historically, the rationale for changing from analytical to iterative image reconstruction was related to better noise-modeling of the latest iterative algorithms. Shepp and Vardi [START_REF] Shepp | Maximum Likelihood Reconstruction for Emission Tomography[END_REF] were the first to propose an maximum likelihood expectation maximization (MLEM) method that showed significant improvement over the previous analytical filtered backprojection (FBP). Today the most used image reconstruction algorithm is based on their first model and its variations (list mode (LM)-MLEM, ordered subset expectation maximization (OSEM)).

In this thesis we present three-gamma (3-g) imaging, where the acquisition system relies on a b+ and g emitter. The reasoning behind 3-g imaging is that the third gamma detection information may help to provide better localization of the annihilation point, thus enabling higher image quality of reduced dose delivered to the patient. Detectors taking advantage of a third gamma emission have been proposed in a number of studies [START_REF] Kacperski | Three-gamma annihilation imaging in positron emission tomography[END_REF][START_REF] Lang | Sub-millimeter nuclear medical imaging with high sensitivity in positron emission tomography using b + g coincidences[END_REF][START_REF] Masełek | Towards 2+ 1 photon tomography: Energy-based selection of two 511 keV photons and a prompt photon with the J-PET scanner[END_REF][START_REF] Moskal | Feasibility study of the positronium imaging with the J-PET tomograph[END_REF], [START_REF]J-PET: A New Technology for the Whole-body PET Imaging[END_REF][START_REF] Krzemie Ń | J-PET analysis framework for the prototype TOF-PET detector[END_REF]. Among these efforts, in this work we present the xenon medical imaging system 2 (XEMIS2) camera, a preclinical 3-g imaging scanner developed in Nantes, France. Two main characteristics of XEMIS2 make it suitable for 3-g imaging: the use of liquid xenon (LXe) as scintillation medium, that offers a continuous geometry, and the use of Scandium-44 (Sc-44), a b + g emitter, which enables the 3-g detection.

The aim of this work is to present a new image reconstruction method based on 3-g detection, which could allow better image quality and thus a reduction of the administered dose. The principle on which 3-g image reconstruction is based is the intersection between a LOR and a Compton cone belonging to the same event. From the two coincidence photons we obtain a LOR, while from the third gamma, emitted as a prompt gamma, we determine a cone. The intersection between the LOR and the cone provides the LOR/cone intersection (LCI) coordinates that locate the most probable annihilation position on the LOR, as for the time difference in TOF, whence the pseudo-TOF denomination of the proposed reconstruction method. In this case, we proposed to carry on a study on the distribution sigma, which depends not only on the scanner resolution but on multiple factors (positron range, energy resolution, scattering, etc). This work is composed of five main chapters in addition to general introduction and conclusion. In Chapter 1 we introduce nuclear imaging, starting with the physics behind nuclear medicine: positron decay, positron-electron annihilation and other positron interactions with matter, positron range, photon attenuation and interactions. The second part of the chapter is about PET imaging, we illustrate how the coincidence photons are detected, we present the challenges concerning scattering and random coincidences, and also the sensitivity and the depth of interaction (DOI) problem in PET. In Chapter 2 we introduce Image Reconstruction in a historical way: first, discussing data acquisition and then talking about analytical image reconstruction algorithms and model-based ones. At the end of the chapter we present the system matrix and all the effects that can be included in the modelization, as the detector geometry, the attenuation, the normalization, the scattering, and the positron range. In Chapter 3 we present the proposed 3-g reconstruction method. We begin addressing our attention to the XEMIS2 project, we discuss the detector and its geometry, the used radioisotope, the data acquisition, and detection process. Then we expose the proposed image reconstruction technique, based on 3-g detection, including the LCI computation, and all the needed requirements as a matter of spatial and energy resolution and, even more important, as angular resolution. The second part of the chapter deals with the Geant4 Application for Emission Tomography (GATE) Monte Carlo simulation used to obtain data. We show the simulated phantoms and talk about customizable and advanced software for tomographic reconstruction (CASToR), the reconstruction software, in which we implemented the 3-g reconstruction algorithm. In Chapter 4 we show the results obtained from the simulation, and the image analysis that was carried on. First, we present the figures of merit used in the analysis, and then the results: the reconstructed image using 3-g and conventional PET method and some metrics and profile plots. In Chapter 5 we discuss the positron range correction issue related to the use of Sc-44. We talk about the state-of-the-art techniques to correct the effect and we show what has been done in our case using variant and invariant kernel convolution, implemented directly in CASToR reconstruction software. At the end of this work, in the conclusion 5.6, we discuss all that has already been done for 3-g image reconstruction and possible future developments among which is the use of Neural Networks to correctly compute the LCI point.

Notation: vectors will be denoted by bold italic letters (e.g. v), matrices will be written in upper-case bold italic letters (e.g. A), with a ij being the entry at row i, column j.

Chapter 1

General principles of nuclear imaging Summary

In this chapter, we introduce the main subject of this work and we gather all the necessary information to fully understand the 3-g image reconstruction method, starting with the introduction to PET imaging. We first present the physics of PET, the radioactive b + decay and the interaction of particles with matter, in particular positron and photon's interaction at low energies. We then introduce the basics of PET imaging, as a matter of detection, coincidence discrimination, and all the challenges encountered with detector sensitivity and DOI. At the end of the chapter we mention TOF-PET introducing briefly the theory and the basic notions that will be helpful in the following of this work.

Introduction

Nuclear imaging is based on the use of radioactive compounds, composed of a molecule labeled with a radioisotope, that follow physiologic or pathological processes linked to a variety of diseases. There exist more than 100 different types of diagnostic exams for nuclear medicine, that can detect early stages of oncological or neurological diseases, when cures are more effective. Unlike other medical imaging procedures, which map the anatomy of the body with almost no information about the metabolism, nuclear imaging, and in particular PET imaging, is a valuable mean for providing disease-related quantitative and qualitative information.

Physical principles

In this section, we discuss some of the basic principles of nuclear physics, especially the ones in which we are most interested, as radioactive decay and photon and positron's interactions with matter.

Radioactive decay

There exist six different modes of radioactive decay [START_REF] Pfützner | Radioactive decays at limits of nuclear stability[END_REF], but the only one interesting in PET imaging is positron decay, also called b + decay. More generally, in PET, radiotracers are labeled with positron emitters (see Table 1.2) that are unstable elements decaying through b + emission. Positron decay consists of the conversion of a proton into a neutron with the emission of a positron and a neutrino, according to the following expression:

p ⇤ ! n + e + + n + energy. (1.1)
where the ⇤ indicates that the proton is not a free particle but belongs to a nucleus.

And furthermore through this equation involving the parent and daughter nuclide:

A Z X b + ! A Z 1 Y. (1.2)
A minimum transition energy requirement is expected for b + decay, since the nucleus reduces its atomic number by one and after losing a positron, the daughter-atom loses an electron to reach its ground state. The mass of both electron and positron being respectively particle and anti-particle, is 511 keV, from which the total minimum of 1.022-MeV energy is required [START_REF] Robley | The atomic nucleus[END_REF].

As the positron travels in the surrounding tissue, it loses its kinetic energy in Coulomb interactions with the atomic electrons. As the rest mass approaches the electron mass, the trajectory deviation may be accentuated, giving an intricate and twisted path to the particle. When the positron reaches thermal energy (few keV), it interacts with an electron by annihilation or forming a positronium atom with a lifetime of about 0.12 ns [START_REF] Louis | Scattering of Positrons by Hydrogen Atoms and Formation of Positronium[END_REF]. Two forms are possible in its ground state: ortho-positronium (parallel particles' spins)

or para-positronium (anti-parallel particles' spins), the latter decays in selfannihilation producing two 511-keV photons as for positron-electron annihilation, while the former annihilates producing three gammas. The back-to-back emission of free annihilation is required for the momentum conservation in the electron-positron pair. Since the photons emission takes place in a frame of reference that moves with a different velocity than the one of the detection system, the photons can be emitted in a slightly different direction from the ideal 180 by a few tenths of degree (s ⇡ 0, 5 ) [START_REF] Dryzek | An undergraduate experiment to test relativistic kinematics using in flight positron annihilation[END_REF].

FIGURE 1.1: Positron emission and annihilation scheme for F-18 [START_REF] Bailey | Physics and Instrumentation in PET[END_REF].

Mathematics of radioactive decay

Nuclear radioactive decay obeys the following equation:

dN dt = lN(t). (1.3)
in which dN is the number of nuclei decaying in a very short amount of time dt, l is the decay constant, and N(t) the number of nuclei at time t. The decay constant does not depend on external factors, but only on the effective decay of the nuclei and it is not tabulated. Instead, we refer to the half-life, defined as the amount of time it takes to reduce the original activity by a factor two:

T 1/2 = ln(2) l .
(1.4)

Charged particle's energy loss

A charged particle of mass M 1 and charge Z 1 penetrating a material of atomic number Z 2 and atomic mass M 2 , slows down losing energy to the medium's atoms. It interacts with the matter via electromagnetic force with the electrons and protons and via strong nuclear force with the nuclei. The particle's stopping process involves complex interactions, although we can measure the particle's average energy loss per unit pass length, called stopping power, dE dx [START_REF] Saadi | Development of a new code for stopping power and CSDA range calculation of incident charged particles, part A: Electron and positron[END_REF]. The computation of this quantity requires a complete knowledge about the particle energy loss mechanisms and it changes drastically from particle to particle. In particular between ions and electrons/positrons, due to the difference in mass.

In this work, we are mostly interested in the energy loss of the positron, which is mainly caused by two mechanisms: Coulomb interaction with orbital electrons of the nuclei or collision stopping power (CSP) and Bremsstrahlung emission or radiation stopping power (RSP) [START_REF] Ça | A study on the calculation of stopping power and CSDA Range for incident positrons[END_REF]:

✓ dE dx ◆ tot = ✓ dE dx ◆ c + ✓ dE dx ◆ r . (1.5)
The first term is the CSP and the second term is the RSP. CSP was firstly studied by Bethe and Bloch [START_REF] William R Leo | Techniques for nuclear and particle physics experiments: a how-to approach[END_REF] and later on by Bhabha [START_REF] Hj Bhabha | The scattering of positrons by electrons with exchange on Dirac's theory of the positron[END_REF], which formalized it as follows:

✓ dE dx ◆ c = KZ 2Ab 2 m e c 2 b 2 g 2 (m e c 2 g 1 2 I 2 ! + 2 ln 2 b 2 12 ✓ 23 + 14 g + 1 + 10 (g + 1) 2 + 4 (g + 1) 3 ◆ d ! . (1.6)
where m e is the electron/positron mass, b = v c , with v the particle velocity and c the speed of light; g = 1 p 1 b 2 is the Lorentz factor. The variable I represent the mean excitation energy, Z the atomic number, A the atomic mass of the material, K = 4pe 4 N 0 m e c 2 = 0.307075 is a coefficient, being N 0 the Avogadro number and e the electron charge, and d is the density correction effect.

On the other hand, the RSP was studied by Lindhard and Ritchie [START_REF] Lindhard | On the properties of a gas of charged particles[END_REF][START_REF] Rh Ritchie | Interaction of charged particles with a degenerate Fermi-Dirac electron gas[END_REF].

In the positron case, due to the small mass, the particle deviates from its trajectory when in the nucleus' electric field, with a radial acceleration causing the Bremsstrahlung emission. When the particle reaches the critical energy, that is when the radiative energy loss equals the collision energy loss, the Bremmstrahlung becomes dominant. The computation of the energy loss due to radiation emission is rather complex and we define it through a parameter, called radiation length X 0 , that corresponds to the distance over which the positron energy is reduced by 1 e due to radiative loss only [START_REF] Donald | Passage of particles through matter[END_REF]:

✓ dE dx ◆ r = E X 0 .
(1.7)

Positron range

A positron is emitted during nuclear decay with an energy that spans from 0 to the maximum available energy, following an asymmetrical distribution centered around half of the E max . In Table 1.1 we see some reference values for the most common PET radionuclides. Due to its initial kinetic energy, the positron path can be not negligible for high energy positrons. The positron range depends also on the material in which it propagates since the interactions probability is calculated on the atomic number and atomic mass.

Positron range is one of the main sources of blur in PET imaging: for F-18, the most used PET radionuclide, we deal with 0.5-mm range correction in the human body, while for Sc-44, the radioisotope chosen for the image reconstruction method proposed in this work, the correction is around 2 mm [START_REF] Ferguson | Comparison of scandium-44 g with other PET radionuclides in pre-clinical PET phantom imaging[END_REF].

We show in Fig. 1.2 a comparison between the energy spectra of the two isotopes. The energy distribution function used to obtain the histograms can be approximated following Levin and Hoffman's formula [START_REF] Craig | Calculation of positron range and its effect on the fundamental limit of positron emission tomography system spatial resolution[END_REF]:

N(E)dE = pF(Z 1, E) ✓ 1 + E 0.511 ◆ (E max E) 2 dE. (1.8)
where Z is the atomic number of the nucleus that undergoes b + decay, E is to the kinetic energy of the positron in MeV, E max is the maximum kinetic energy in MeV, p is the momentum of the positron p =

q 1 + E 0.511 2 1
and

F(Z, E) is the Fermi function: F(Z, E) = 2ph 1 e 2ph .
(1.9)

in which h = Za p 1 + E 0.511 , a = 1 137 being the fine structure constant. Due to the small correction, positron range for F-18 can be neglected for all those cases in which the scanner resolution is higher than the correction itself [START_REF] Emond | Effect of positron range on PET quantification in diseased and normal lungs[END_REF].

Isotope Mean positron range Isotope Mean positron range [mm]

[mm] 18 F 0.5 15 O 2.5 64 Cu 0.57 60 Cu 3.09 13 N 1.31 76 Br 3.09 61 Cu 1.32 62 Cu 4.39 44 Sc 2.4 66 Ga 6.13 The blurring introduced in the reconstructed images by the positron range effect is due to the displacement of the LOR. In

Photon interactions with matter

High-energy photons are secondary ionizing radiation, meaning that they interact with atoms, nuclei, and electrons without causing direct ionization. The interaction can be seen as a collision, that results in the ejection of an electron, which can cause ionization. There are four significant photon-matter interactions in nuclear medicine: photoelectric effect, Compton scattering, pair production, and Rayleigh scattering.

The photoelectric effect consists of a collision between a photon and an atom. The incident photon is absorbed by the atom and all its energy is passed over to an orbital electron. If the transferred energy E 0 is higher than the binding energy E b , this particle, also called photoelectron, can be ejected with an energy given by the following equation:

E pe = E 0 E b .
(1. [START_REF] Brix | Radiation exposure of patients undergoing wholebody dual-modality 18F-FDG PET/CT examinations[END_REF] this interaction can in turn create a vacancy in an orbital shell, which leads to When the collision occurs between a photon and an atomic external electron, we deal with Compton scattering. If the photon's deposited energy is higher than the electron's binding energy, then the electron is ejected as a recoil electron. In any case, the photon is only deflected, not absorbed, with an angle q C and an energy that is calculated as in the following equation:

E 0 = E 0 1 + ⇣ E 0 m e c 2 ⌘ (1 cosq C ) . (1.11)
where E 0 is the incident photon's energy, m e is the mass of the electron and c is the speed of light. Thus the Compton scattering angle can then be calculated as follows:

cos q C = 1 + m e c 2 ✓ 1 E 0 1 E 0 ◆ . (1.12)
The transferred energy range goes from 0 keV to a maximum value E max .

The extreme cases occur when there is a back-scattering event, see Fig. 1.4.

Compton scattering is not equally probable at all energies or for all scattering angles, the interaction probability follows the Klein-Nishina's formula:

ds dW = Zr 2 0 ✓ 1 1 + (1 cosq C ) ◆ 2 ✓ 1 + cos 2 q C 2 ◆ ✓ 1 + (1 cosq C ) 2 (1 + cos 2 q C )(1 + (1 cosq C )) ◆ .
(1.13)

where ds dW is the differential cross-section, Z the atomic number of the material, r 0 the electron radius. For energies higher than 1.022 MeV we can deal with pair production.

This consists of the interaction between a photon and a charged particle's electric field, in which the photon is converted into an electron-positron pair.

The requirement of minimum energy is due to the sum of the electron and positron mass, which accounts for 1.022 MeV. The difference between the Chapter 1. General principles of nuclear imaging incident photon's energy and the sum of the two masses is equally shared between the new particles pair as kinetic energy.

Another kind of scattering named coherent or Rayleigh scattering involves a photon and an atom: due to the great atomic mass almost no energy is absorbed in the collision as recoil energy, and the photon is scattered with very few energy loss.

These four main photon interactions do not occur at all energies and in all materials with the same probability. We can see in Fig. 1.6 how for lower energies in low Z materials the photoelectric effect is the most probable interaction, while for energies higher than 1.022 keV pair production is predominating.

Compton scattering is most likely to happen for 511-keV photons in an organic body, however in case of small diffusion angles the two gammas can still be recorded as coincidence photons in PET. Scattered events represent between 40% and 60% of the total events in a scanner acquisition [START_REF] Dj Thompson | Calibration of the energetic gamma-ray experiment telescope (EGRET) for the Compton gamma-ray observatory[END_REF].

FIGURE 1.6: Main photon interactions versus photon energies for different atomic number [START_REF] Sorenson "simon | Physics in Nuclear Imaging[END_REF].

Photon attenuation

When a photon beam passes through an object, some of the photons interact with the material and the beam intensity is attenuated. This effect depends on four factors: thickness, density, atomic number, and the photons energy. The beam intensity follows this equation widely known in literature [START_REF] Heinrich | Photometria sive de mensura et gradibus luminis, colorum et umbrae[END_REF][START_REF] Beer | Bestimmung der absorption des rothen lichts in farbigen flussigkeiten[END_REF][START_REF] Donald | The beer-lambert law[END_REF]: where I 0 is the initial intensity of the beam, µ is the linear attenuation coefficient depending on the density r, on the atomic number Z, and on the photons energy, and x is the thickness of the object. The linear coefficient can be seen as a measure of the probability of photon attenuation in a material unit length [µ]=cm 1 .

I x = I 0 e µ(

PET imaging

PET imaging is used to obtain the uptake distribution of a radiopharmaceutical in the patient's body, mostly in oncological studies, before and after treatment and during radiotherapy. PET scans are also important diagnostic tools in neurology for degenerative disease, for example Alzheimer's. See Table 1.2 for some common radiotracers and their conventional applications. The radiation reaching the detector is converted in electric pulses, which are collected in the acquisition system to be analyzed and reconstruct the image.

The detection is based on the knowledge of these two photons' characteristics: co-linearity and 511-keV of initial energy. Several factors can influence the correct event detection, among them we find scattering in the body and/or in the detector and positron range. The scintillating material needs to be coupled to an electronic system equipped with a light detector to collect the radiation light yield. There are two principal ways: the one-to-one coupling technique (one crystal for one photo-detector) and the block detector one (multiple crystals coupled to the same photo-detector through a light guide). The coupling system to adopt is chosen upon a study on the photodetector and on the scintillator to employ.

If a light guide needs to be used, the properties of the transparent material ought to match the scintillator: similar refraction index, almost no reflection and so on.

The event detection and acquisition process occurs as follows: a photon scintillate in the detector producing light, that is converted into a pulse from the PMT, giving information about the position of the incident photon. Two pulses belonging to the same time window are recorded together as a coincidence. The line connecting the two detected events is called LOR and it joins two sides of the detector, crossing the field of view (FOV), where the patient is located.

Scintillation in LXe

The way photons are produced depends on the type of scintillator we are using; in LXe this can happen via atomic ionization or atomic excitation [START_REF] Miyajima | Numbers of scintillation photons produced in NaI(Tl) and plastic scintillator by gamma-rays[END_REF][START_REF] Doke | Let dependence of scintillation yields in liquid argon[END_REF]. The atomic excitation consists of the excitation of an atom in the xenon by a photoelectron, that combines with another xenon atom creating an excited dimer. After few picoseconds, the dimer de-excites and emits an ultraviolet photon. For the atomic ionization case, the photoelectron creates an electronion pair with one of the xenon atoms, that combines with another atom producign an ionized dimer. Eventually, a visible photon will be emitted in the same way as in the atomic excitation case [START_REF] Gallego | Optimization of a single-phase liquid xenon Compton camera for 3g medical imaging[END_REF].

Photodetectors

The conversion between radiation and visible light is performed by a photodetector, typically a PMT. The total number of electrons that reach the end of the PMT for each photoelectron defines the gain of the PMT, which is normally around 10 5 -10 8 [START_REF] Glenn | Radiation detection and measurement[END_REF]. The current pulse contains both time and energy information about the arrival of the photon on the PMT photocatode.

Other photodetectors based on semiconductors (i.e. Silicon PMTs) are also widely used in medical imaging due to their insensitivity to the magnetic field.

Coincidences

In PET imaging we are interested in the collection of coincidences, since we aim at detecting the two 511-keV photons produced by the same annihilation.

Due to several factors (e.g. spatial and energy resolution of the detector, scattering in the phantom, positron range, etc.) the acquired coincidences do not always cross the real annihilation point.

In a classical PET acquisition we can find three types of events: true, random, and scatter coincidences. With true coincidence we denominate the ideal case in which two photons that do not undergo significant interactions before being acquired, track a correct LOR. In the other two cases, we have an effect of LOR mispositioning, that carries an error on the activity distribution.

A scatter coincidence occurs when one or both photons are scattered before being detected and it is more likely to happen in large objects, see Fig. 1.9, example 1. We refer to a random event if the coincidence photons do not belong to the same annihilation event, see Fig. 1.9, example 2 and 3. Random events can be due to a large time window or to multiple annihilations occurring in the same time window.

These events add bias to the reconstructed image if not accounted for properly: for the randoms it increases proportionally to the amount of activity A as 2tA 2 , where t is the coincidence resolving time of the system. We can limit their acquisition by working with a narrow time window and with a good time resolution scintillator. Moreover, for scattered events we can work with an energy window (> 400 keV), rejecting events that lost already a conspicuous part of their energy, and use detectors with a high energy resolution.

FIGURE 1.9: Event 1 represents a double scattered event, events 2 and 3 show a random coincidence. [START_REF] Bailey | Physics and Instrumentation in PET[END_REF]

Sensitivity and Depth Of Interaction

The term sensitivity referred to a PET system indicates the ability to detect photons and it mainly depends on two characteristics: the detector stopping FIGURE 1.10: Parallax error due to the unknown event DOI [START_REF] Bailey | Physics and Instrumentation in PET[END_REF].

power and the scanner geometry. High sensitivity is often related to a large solid-angle coverage in the scanner, as in the case of narrow and long axial FOV, and of a material with high atomic number.

The ability of efficiently contain and detect photons is useful not only to improve the sensitivity but also to account for another phenomenon, called the DOI effect. When a photon enters the detector, its interaction position is determined on the detector entrance surface and it is not the actual interaction point. If the detector is sufficiently deep for these coordinates to not match the actual interaction point, we have an effect of LOR mispositioning, that increases with the depth, as seen in Fig. 1.10. This effect is important especially for deep detectors and in the oblique coincidence case, for which every recorded event carries a systematic error that needs correction. The idea of using the arrival time of the two photons to locate the annihilation point developed in the 80s, when the first TOF-PET systems were built.

Time of Flight PET

Such systems were designed to improve the SNR of the reconstructed images and reduce the random rate acquisition.

In TOF-PET we impose a time window on the coincidence detection (few ns), in order to accept photons belonging to the same annihilation, and we record the time difference between the two. If we assume to have an ideal detector with a perfect time resolution, we can determine the location x of the event on the LOR as follows:

x = cdT 2 .
(1.15) In the ideal case, using only the LOR coordinates and the time difference, we can determine the annihilation point in a 3D volume, without passing through image reconstruction [START_REF] Tom K Lewellen | Time-of-flight PET[END_REF].

In reality we need to consider the actual time resolution of the detector, which gives us the uncertainty on the position computed through Eq. (1.15).

We can see in Fig. 1.12 how the localization through TOF is considered the center of a Gaussian distribution, whose FWHM depends on the time resolution. For example, with a 600 ps timing resolution our uncertainty on the LOR would be around 9 cm.

The main difference compared to conventional PET is that in classic reconstruction the probability associated with the LOR voxels is the same for each of them, while in the TOF case we give a Gaussian-shaped distribution to the LOR.

Looking at the benefits carried by this technology, we could list noise reduction, especially for little objects [START_REF] Thomas | Time-of-flight positron emission tomography: status relative to conventional PET[END_REF], the earlier convergence of the iteration method and at a clinical level, the reduction of the acquisition time due to the improvement in image quality. As for the SNR reduction compared to non-TOF-PET, the gain factor can be seen from the following equation from [START_REF] Thomas | Time-of-flight positron emission tomography: status relative to conventional PET[END_REF]:

SNR TOF SNR non TOF = r D Dx . (1.16)
where D is the size of the patient body and Dx the spatial equivalent to the TOF measurement. The gain factor depends on the size of the emitting object, thus we expect a higher value for heavier bodies. For a 20-cm-diameter object

(

Introduction

Tomographic image reconstruction aims to form images of the radiotracer distribution in the patient's body starting from raw PET data. In the data acquisition process we collect LORs through projections or coordinates pairs and the image reconstruction process can be seen as the estimation of the function f (x): R n ! R, with n = 2, 3 for a 2D or a 3D image. f (x) gives the value of the image in the point x and it corresponds to the activity distribution, following this equation:

E det.photons/second = Z Z VOR s(x) f (x)dx. (2.1)
where the integral is on the volume of response (VOR), the parallelepiped that joins the pair of detector in which the coincidence occurred, and s(x) is the scanner sensitivity.

If we consider the ideal case without attenuation, random or scattered coincidences, and detector geometry effect, we can assume that the total number of events detected in this volume is proportional to the amount of activity in the same.

In image reconstruction, there are two types of algorithms to obtain images:

Data acquisition

In PET imaging, data acquisition consists of collecting the information on the coincidence photons, looking for the detection of colinear 511-keV gammas in the time and energy acceptance windows.

In order to collect a high number of events and to maximize detector sensitivity, a PET system is typically composed of multiple ring-shaped detectors with a wide-volume coverage around the patient, positioned in the center. All events are recorded, including scattered and random coincidences, and raw data are stored in sinograms [START_REF] Fahey | Data acquisition in PET imaging[END_REF][START_REF] Leahy | Recent developments in iterative image reconstruction for PET and SPECT[END_REF].

In the data acquisition process detectors produce coincidence pairs with opposed detectors (axially and transaxially), which form a LOR identified by the rotational angle f and the radial displacement s, see Fig. During the reconstruction process, we need to consider a correction factor for LORs' spacing, since they tend to be closer together as we approach the edge of the ring; this effect is not considered during the acquisition.

In TOF-PET the difference in time between the coincidence photons arrival is also measured and stored. For a LM acquisition the dt-value is stored event by event, while in sinogram mode the sinogram is divided into several TOF-bin, each of them corresponding to a time difference.

Deadtime

All detectors present a limit due to the data treatment speed. When a PET system starts the event acquisition, there is a time interval in which other events cannot be detected, due to the system saturation while treating the first impulse. This amount of time is defined as the deadtime of the system. It is mostly due to the time the electronics take to register the first signal, but it can also depend on the scintillation time of the detector's material. This time gap can be seen as a blind time for the detector, during which we are unable to treat new interactions, and it is characteristic for each scanner.

In most high count rate situations another effect can affect measurements, known as pile-up. This happens when two signals arrive too close to each other and are summed up as one signal only, leading to energy resolution degradation and a count rate loss [START_REF] Usman | Radiation detector deadtime and pile up: A review of the status of science[END_REF].

Analytic Image Reconstruction

Analytical approaches are based on the hypothesis that for every registered event we have a measure of the line integral associated that describes the radioactive distribution on the line. The first to solve this problem was Radon, in 1917 [START_REF] Radon | On the determination of functions from their integrals along certain manifolds[END_REF]. He proposed what we call the Radon transform, illustrated in Fig. 2.2:

R f (s, f) = Z L s,f f (x, y)dxdy = Z • • f (s cos f l sin f, s sin f + l cos f)dl.
(2.2) in which L s,f is the unit sphere, s and f are the radial displacement and the angle that determine the LOR and f (x, y) is the function to be determined.

Actually, the integral is constrained by the FOV size and it occurs only between R FOV and R FOV .

In TOF-PET systems the same equation becomes:

R TOF f (s, f, t) = Z • • f ( s sin f + l cos f, s cos f + l sin f)h(t l)dl. (2.3)
where t is the time difference and h(t) is the Gaussian kernel used to modulate the probability on the LOR. R f ( s , ⌃ ) Analytical image reconstruction is considered to give direct results and linearity, which allows to better control noise correlation during the reconstruction process; and for its speed, mandatory in large 3D data sets [START_REF] Alessio | PET image reconstruction[END_REF].

2D analytic image reconstruction

Some properties of the Radon transform need to be stressed out: first, the operator R is linear

R[ f + lh] = R f + lRh, ( 2.4)
It is also continuous and invariant for rotations [START_REF] Bailey | Physics and Instrumentation in PET[END_REF]:

R f ( s, f + p) = R f (s, f). (2.5)
This is valid for scanners that collect all the LORs crossing the object, corrections to raw data need to be applied in case of incomplete scanner ring, to cover for the missing part of data, for example through interpolation.

The main objective of analytical reconstruction is to reconstruct an image f , given a sinogram g, for which:

g(s, f) = R f (s, f). (2.6)
meaning that we need to find R 1 the inverse operator of R.

Let us introduce the central slice theorem [START_REF] Natterer | Mathematical methods in image reconstruction[END_REF][START_REF] Bracewell | The projection-slice theorem[END_REF], which ties the 2D Fourier For all (s, f) 2 Z:

g f (s) = R f (s, f), ( 2.7)
we can write then:

ĝf (s) = f (sw(f)). (2.8)
where s = r • w(f) with r 2 R 2 and w(f) = (cos f, sin f).

For a fixed angle f we can write:

ĝf (s) = Z R e 2ipss g f (s)ds = (2.9) = Z R e 2ipss Z R f (sw ? (f) + tw(f))dtds. (2.10)
Let r be equal to sw ? (f) + tw(f), dr = dtds and r • w(f) = s:

ĝf (s) = Z 2 R e 2ipsr•w(f) f (r)dr = f (sw(f)). ( 2

.11)

We can now find the values of f (x, y) relative to the image through an inverse Fourier transform, this algorithm is also called direct Fourier reconstruction:

f (r) = 1 2 Z 2p 0 Z • • |s| ĝf (s)e 2ipsr•w(f) dsdf. (2.12)
where |s| is a ramp filter. Ramp filtering can be seen as a convolution in the spatial domain or as a multiplication in the Fourier domain.

Back-projection

Two of the main operators in image reconstruction are back-projection and its adjoint, forward projection. They can be seen as the Radon transform (forward projector) and its adjoint, the back-projection operator R ⇤ , defined as:

R ⇤ p(s, f) = Z 2p 0 p(x cos f + y sin f, f)df. (2.13)
and for a fixed angle f is also valid:

b 1 (x, y; f) = p(x cos f + y sin f, f). (2.14)
where b 1 (x, y; f) is the 2D function formed through the back-projection of the single projection indexed by f.

In reality, we do not observe directly g = R f , but a noisy version of it:

g e = R f + e. (2.15)
Noise is amplified by the back-projection operation, in fact if we take the norm:

||R ⇤ [g + e] R ⇤ g|| X = ||R ⇤ e|| X . (2.16)
where X is the image space, R ⇤ g = f and ||R ⇤ e|| X does not tend to zero when e ! 0. A small perturbation on g can disrupt the reconstruction of f .

A solution to this problem is to apply a ramp filter |s|, the Jacobien of the transformation, which amplifies high frequencies corresponding to the image details, since high spatial resolution requires high frequency information. The ramp filter in 2D is obtained through apodization and discretization, using the fact that we can sample |s| resulting in a discrete ramp-filter [START_REF] Gengsheng | Can the backprojection filtering algorithm be as accurate as the filtered backprojection algorithm?[END_REF].

This method is called back-projection filtering or FBP [START_REF] Avinash | Principles of computerized tomographic imaging[END_REF], depending on when the filter is applied. The ramp filter could not be sufficient due to its behavior that ensures high spatial resolution at the expense of noise, especially since it carries the FBP is the most traditional image reconstruction approach. It is based on an idealized model that does not consider noise, scatter, attenuation, and other real data important features. In fact, it assumes that the number of detected events in a specific direction can give an approximation of the radiotracer distribution along the same direction, that is one of the projections introduced earlier. The only way to consider noise in FBP is smoothing the projection before image reconstruction.

f (x, y) = Z p 0 R F f (s, f)df. (2.17)
where R F f (s, f) = F 1 1 {|s|F 1 {p f (s)}} is the filtered projection, |s| being the ramp filter. The FBP method works as follows: first, the data are collected through a projection for a specific angle and a sinogram is created; then the algorithm exploits the back-projection operation to smear each projection back to the region along with the acquired direction (or angle) of acquisition.

Adding the back-projections for all the angles, we will obtain a blurred approximation of the original object. Sharpening the projections is also possible, using a Fourier filter [START_REF] Paul E Kinahan | Analytic image reconstruction methods[END_REF]; although analytical methods like FBP are mostly used to produce fast and practical solutions that highlight issues about data acquisition. Techniques which are used to implement real statistical estimations are of the iterative kind and they include the corrections already in the reconstruction process [START_REF] David | Iterative image reconstruction[END_REF].

3D analytic image reconstruction

One of the main 3D reconstruction's challenges is the time-consuming computation. Due to this problem, the first reconstructed 3D volumes were ensembles of independent 2D reconstructions [START_REF] Lee A Feldkamp | Practical conebeam algorithm[END_REF]. Other methods were also used to avoid long computation time as for example transforming 3D projections into 2D ones through rebinning, single-slice rebinning [START_REF] Margaret E Daube-Witherspoon | Treatment of axial data in three-dimensional PET[END_REF], multi-slice rebinning [START_REF] Robert M Lewitt | Three-dimensional image reconstruction for PET by multi-slice rebinning and axial image filtering[END_REF],

direct Fourier [START_REF] Stark | Direct Fourier reconstruction in computer tomography[END_REF] and Fourier rebinning [START_REF] Defrise | Exact and approximate rebinning algorithms for 3-D PET data[END_REF].

Image reconstruction algorithms using directly 3D sinograms are called fully 3D and the methods seen for 2D can be extended to this category. In order to generalize them, we need to introduce a representation for the fully-3D data, the Radon transform operator in 3D can be written as follows, see Fig 2 .4:

R f (t, ô) = p(t, f, q) = Z Z Z R 3 f (x)d(x • ô(f, q) t)dx.
(2.18)

where d(x) is the Dirac function and the vector ô(f, q) locates the plane related to t 2 R, signed distance between the origin and the perpendicular plane to ô(f, q) 2 S 3 , where S 3 is a 3D unit sphere. Before approaching the analytic algorithms for 3D data, we need to introduce the 3D generalized central slice theorem, that states the equivalency between the 2D Fourier transform of a projection of arbitrary direction ẑ(f, q)

and the central section of the 3D Fourier transform of the object at the same angle. This implies that the image can be reconstructed from a set of nontruncated projections, if the set has an intersection with every equatorial circle on the sphere S 2 . What is known as Orlov's condition assesses that the image can be reconstructed from the set of non-truncated projections if there is no great circle on the unit sphere that does not intersect the set of non-truncated 2D parallel projections. We can imagine that the acquisition is repeated for a large range of z and once we get the stack of images from different planes we can form a 3D image. This is a way of considering 3D imaging, although the fully 3D imaging includes the acquisition on all the planes perpendicular to the scanner axis, called direct planes, and on all the oblique planes that cross the direct planes as well, in Fig. 2.5 some examples of different 3D reconstructions.

Back-projection

For the 3D back-projection operation, it is often more efficient to back-project keeping a fixed direction (fixed angles f and q) and placing back the values of the projection into a 3D array along the corresponding LORs.

b 1 (x, y, z; f, q) = p(r • ô(f, q), f, q). (2.19)
where r = (x, y, z). The complete back-projection can be obtained integrating in df and dq.

To better understand the parallel with the 2D case, we can derive the FBP equation to be used in reconstruction, for which the back-projection operator can be written as:

R ⇤ g(r) = Z S 2 g(r • ô, f)d ô. (2.20)
Filtering in 3D becomes more complex: we can see the 2D filter as the crosssection of the 3D one. For symmetry the filter will be independent of the projection angle f, but dependent on q and on the maximum acceptance angle [START_REF] Colsher | Fully-three-dimensional positron emission tomography[END_REF].

If we write g ô(t) = R f (t, ô), we obtain:

f (r) = 1 2 Z S 2 Z • • |s| 2 ĝô (s)e 2ipsr• ôdsd ô. (2.21)
Reconstruction with analytic techniques remains more efficient for very large 3D data sets, especially when multiple data sets are acquired in wholebody or dynamic studies, due to the amount of required time. However the speed is not the only FBP strength: analytic algorithms are linear and thereby allow easier control of the spatial resolution and noise correlations in the reconstruction.

Model-based Image Reconstruction

Tomographic image reconstruction can be seen as the operation of inversion of the Radon transform. This simplistic approach can be solved through FBP, without considering random detection, scattering in the tissues and the acquired data's noise properties. Unlike analytic approaches, statistical methods incorporate all known information in a discrete framework: through A the SM, detailed in section 2.6, we are able to discretize the problem and consider attenuation, noise, scattering, detector response, and a priori information about the images.

For analytic methods we required some common steps to reconstruct (i.e.

filtering and back-projection), for iterative methods, we need five elements:

-Basis function: the whole image is discretized and composed of basis functions, for example, voxels. This way the image can be represented as a vector of basis function coefficients (i.e. the voxel's uptake), and we deal with a discretized problem in the data and image domain.

-System model: a model including the scanner geometry effects and the physical effects on data acquisition.

-Noise model: this model concerns the data deviation from their expected values based on a Poisson model.

-Objective function: the objective function measures the fit between the observed data p and the modeled data p, that need to be optimized.

-Numerical optimizer: this is the mathematical algorithm that maximizes (or minimizes) the objective function.

All iterative methods are based on a first forward modeling step followed by a correction step. The activity is discretized voxel by voxel, every element a ij of the SM corresponds to the probability of a positron emitted in the voxel j to be detected by the crystals' couple i, connected by a LOR.

p = Af + n = p(f) + n. (2.22)
where p is the measured projection, f the radioactive distribution, and n represents the projections' noise, that can be seen as a shifted Poisson variable, since p ⇠ Poisson(p; f ).

As for maximizing (or minimizing) the discrepancy between measured and expected data p, we can take as objective function the least-square (LS)

function [START_REF] Kaufman | Maximum likelihood, least squares, and penalized least squares for PET[END_REF]:

O LS = I Â i=1 (p i pi (f)) 2 . (2.23)
where the sum is on the projection bins' values and p i represents the measured data, while pi ( f ) is the expected LOR data. In case we have some knowledge about the projection data, we can weight each element in the sum differently and obtain the weighted-least squares (WLS) function, in which the weights can be estimated directly from the projections of the current image estimate.

LOR projection data are typically based on Poisson statistics, due to the positron emission following the same probability, thus we can introduce the Poisson distributed objective function from Shepp and Vardi [START_REF] Shepp | Maximum Likelihood Reconstruction for Emission Tomography[END_REF] known as the maximum-likelihood (ML) function:

O ML = L(f) = I ' i=1 Pr(p i | pi (f)) = I ' i=1 ( pi (f)) p i e pi (f) p i ! . (2.24)
where Pr(p i | pi (f)) is the Poisson distribution and it gives the probability of obtaining a measured value equal to p i in the LOR i while pi (f) is the expectation. The multiplication is due to the assumption that all the LOR counts are mutually independent and the probability of obtaining the vector p is the product of the Poisson probabilities L(f), also called likelihood [START_REF] Andrew | Image reconstruction and correction techniques for positron volume imaging with rotating planar detectors[END_REF].

Starting from Eq. (2.24), many algorithms have been proposed to maximize the logarithm of the likelihood:

l(f) = ln L(f) = I Â i=1 (p i ln pi (f) pi (f) ln p i !). (2.25)
where ln p i ! is a constant and it can be omitted. Since the log-likelihood l(f) depends on p, which are unknown measurements, we cannot directly calculate the ML estimate for f, altough there are algorithms to maximize the expected value of l(f):

f(k+1) = argmax f E h l(f)|p, f(k) i . ( 2.26)
ML estimators have interesting properties that make them useful in many situations. First of all they are asymptotically unbiased, their bias is reduced with the growth of the observations. Second, they are asymptotically efficient as they yield the minimum variance for large numbers. Although these properties make them not particularly susceptible to noise, their variance is still high and the images are often filtered to introduce spatial smoothing.

Iterative algorithms

The process of iterative reconstruction begins with an estimated f (k) of the values in the image. The current image estimate is then projected, giving a set of values that would be expected if this was the true image. The predicted projections are then compared to the measured ones to elaborate a set of error values in the projection space. These are back-projected to obtain image-space error values that are used to update the image estimate f (k+1) . And this whole process is reproduced for each iteration.

The main difference with analytical back-projection-based algorithms is that there is no feedback about the image estimate in the direct reconstruction.

Another way to estimate f k is to find the image that satisfies all the constraints imposed by the measured data and the prior knowledge. ART [START_REF] Gordon | Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography[END_REF] is one of the methods using this process: the difference between the measured value and the expected value is calculated for each LOR and through the transpose SM the correction is applied in the image space. Neglecting the noise parameter from Eq. (2.22) we have:

f (k) = f (k 1) a i f (k 1) p i ||a 2 i ||. a i (2.27)
where a i is the SM line associated to the i-th LOR. This algorithm corresponds to what is known as the Kaczmarz method [START_REF] Haller | Kaczmarz algorithm in Hilbert space[END_REF], which allows the iterative solution of a system consisting of N equations in N unknown variables.

Taking Eq. (2.22) we can easily see how it represents a system of I equations of hyperplanes in the J dimensional space. The algorithm projects the current image vector f (k) orthogonal on the hyperplanes and updates the image accordingly, as done in ART. Several algorithms originated from ART, i.e.

SART, SIRT [START_REF] Landweber | An iteration formula for Fredholm integral equations of the first kind[END_REF], MART, all neglecting the noise factor, thus not very efficient in PET imaging.

The most famous iterative algorithm for PET and SPECT is MLEM. It was studied by Shepp and Vardi [START_REF] Shepp | Maximum Likelihood Reconstruction for Emission Tomography[END_REF], Dempster [START_REF] Arthur P Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], Lange and Carson [START_REF] Lange | EM reconstruction algorithms for emission and transmission tomography[END_REF] and others. The basic problem consists in the estimation of an image vector f from the counts measurements stored in a vector g. f i is the number of events taking place in the i-th voxel, and g j is the number of events occured in the j-th LOR connecting two detector elements.

Let A be the SM, for which each element a ij gives the probability of an event detected in the i-th LOR to have occurred in LOR bin j. The relation between the average number of events detected on a LOR and the matrix is then:

E[p] = Af. (2.28)
where f is the image we want to reconstruct inverting this equation.

Let p be the incomplete observed data, incomplete because we only know that an event occurred on a LOR, we do not know exactly from where it originates. P is the matrix having as elements the number of events detected on the i-th LOR that originated from the j-th bin. Being able to observe P directly we could solve our problem as:

fj = Â i p ij . ( 2

.29)

Since the likelihood depends on the elements p ij , which refer to unknown measurements, we cannot directly calculate the ML estimate and obtain f, thus we are going to maximize its expected value. From Eq. (2.26) we can rewrite the expected value inserting the log-likelihood as:

E "  i  j f i a ij + p ij ln f j a ij p, f(k) # , (2.30)
due to linearity we can rewrite:

 i  j ⇣ f i a ij + E h p ij p, f(k) i ln f j a ij ⌘ . (2.31)
For the probability theory we know that for independent Poisson random variables, as p ij and p, the conditional probability distribution given the sum of its values, its a binomial distribution with parameters

⇣ Â j p ij , E[p ij ] Â j E[p ij ] ⌘ [91]
.

With E[p ij ] = f j a ij we obtain:

E h p ij |p, f(k) i = p i f (k) j a ij  n f (k) n a in . (2.32)
Once the expectation value is developed, we can proceed with the maxi-

mization step: ∂ ∂ f l E h l(f) p, f(k) i = 0, (2.33)
 i a il +  i E h p il p, f(k) i 1 f l = 0, (2.34)
f l =  i E h p il p, f(k) i  i a il , (2.35)
f l = f (k) l  i a il  i a il p i  n f (k) n a in , (2.36)
where f l is the l-th component of our estimation. The MLEM equation is the following:

f k+1 l = f (k) l  i a il  i a il p i  n f (k) n a in . (2.37)
In Fig. 2.6 we show the convergence properties of MLEM algorithm in different situations.

One of the characteristics of MLEM is the slow convergence: we need around 30-50 iterations to reach convergence, with a forward projection and a back-projection in each of them, meaning that it is twice as slow as FBP. However, in general, MLEM performs better than any analytic method, due to the possibility of including nonuniform constraints.

Other algorithms can be obtained from similar logic, especially envisioning improvements of the MLEM algorithm, as for example to accelerate convergence and reduce the iteration number. An algorithm following this logic is the OSEM, in which the projections are divided in ordered subsets to reach convergence earlier [START_REF] Malcolm | Accelerated image reconstruction using ordered subsets of projection data[END_REF][START_REF] Verhaeghe | AB-OSEM reconstruction for improved Patlak kinetic parameter estimation: a simulation study[END_REF]:

f (k,m+1) l = f (k,m) l  i2S a il  i2S a il p i  l 0 a il 0 f (k) k,m . (2.38)
where m is the subset index, ans S is the group of subsets [START_REF] Florent | Impact of image-space resolution modeling for studies with the high-resolution research tomograph[END_REF][START_REF] Verhaeghe | Quantification task-optimized estimates from OSEM and FBP reconstructions in single-and multi-subject studies[END_REF][START_REF] Reader | Advances in PET Image Reconstruction[END_REF]. To reach the one-iteration step all the subsets need to be treated, if only one subset is used, we go back to the MLEM algorithm. This modification of the MLEM algorithm leads to a more efficient update, due to the computation being a factor m smaller [START_REF] Hsiao | An accelerated convergent ordered subsets algorithm for emission tomography[END_REF]. The main problem with OSEM is that we face an interations limit corresponding to the number of subsets, thus convergence to an ML estimated is lost.

In all these algorithms the projections p i are considered as the sum of the detected events in a bin, without considering the data format (sinogram or histogram). A variation of the MLEM algorithm was proposed in 1998 by Parra and Barrett [START_REF] Parra | List-mode likelihood: EM algorithm and image quality estimation demonstrated on 2-D PET[END_REF] and later by others [START_REF] Ronald H Huesman | List-mode maximum-likelihood reconstruction applied to positron emission mammography (PEM) with irregular sampling[END_REF][START_REF] Thomas | Spatiotemporal reconstruction of list-mode PET data[END_REF][START_REF] Lamare | List-mode-based reconstruction for respiratory motion correction in PET using non-rigid body transformations[END_REF], in order to consider LM data with no prior conversion: .39) this algorithm is known as list-mode expectation maximization (LMEM) and ML represents the events index.

f (k+1) l = f (k) l  i a il  i2ML a il  l 0 a il 0 f (k) l 0 . ( 2
In Fig. 2.7 we can see some images obtained through OSEM reconstruction for different subsets number at different iterations. It seems that the OSEM algorithm has all the properties of the MLEM algorithm but requires fewer iterations to achieve a stable result, nevertheless, it adds the difficulty of choosing the right number of subsets. All the methods shown in this section are ray-driven, although there exist some voxel-driven techniques. Instead of updating multiple voxel values based on some LOR values, they try to minimize (or maximize) the objective function on a single voxel basis [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF][START_REF] Jeffrey | Penalized weighted least-squares image reconstruction for positron emission tomography[END_REF]. The minimum (or maximum) is found for a single value in the image and all the other values are constant. The procedure needs to be repeated for each voxel in order to end the iteration.

System modelization and corrections

To achieve a quantitative reconstruction in which we obtain the real radiopharmaceutical's uptake in each voxel, we should model all the effects happening during data acquisition: the positron range, the scattering and the attenuation, etc. Evidences from multiple studies [START_REF] Edward J Hoffman | Quantitation in positron emission computed tomography: 1. Effect of object size[END_REF][START_REF] Edward | Quantitation in positron emission computed tomography: 4. Effect of accidental coincidences[END_REF][START_REF] Huang | Quantitation in positron emission computed tomography: 2. Effects of inaccurate attenuation correction[END_REF][START_REF] John C Mazziotta | Quantitation in positron emission computed tomography: 5. Physical-anatomical effects[END_REF][START_REF] Simon | Effects of scatter on model parameter estimates in 3D PET studies of the human brain[END_REF][START_REF] Montandon | Atlas-guided non-uniform attenuation correction in cerebral 3D PET imaging[END_REF] show that all these non-linear corrections can make a difference in the reconstructed image quality.

System matrix

Iterative algorithms model the system through a SM, which is applied to the forward and back-projection operators and its elements can include the corrections for the bias-carrying effects. Due to the matrix size (e.g. 10 14 for a

Philips PET system) and to the computation time required to calculate all the elements, the SM is not entirely stored and often computed on the fly. There exist several models to make the system matrix computation easier, one of which is known as the factored system model [START_REF] Qi | High-resolution 3D Bayesian image reconstruction using the microPET small-animal scanner[END_REF][START_REF] Qi | Effect of errors in the system matrix on maximum a posteriori image reconstruction[END_REF], which considers the system matrix as a multiplication of matrices.

A 2 R M⇥N is the SM in which the element a ij defines the probability of detecting an emission from voxel j, j = 1, • • • , N in the detector i with i = 1, • • • , M, and we can express it as:

A = A sens A blur A att A geom A positron .
(2.40)

where A att 2 R M⇥M is a diagonal matrix containing the attenuation factors, expressed in cm 1 ; A sens 2 R M⇥M is the normalization matrix, a diagonal matrix measured through the acquisition of a homogeneous cylinder;

A positron 2 R N⇥N is the matrix considering the positron range correction;

A geom 2 R M⇥N considers the geometrical sensitivity of the detector and it represents the probability that a photon pair produced in voxel j reaches the detector i in absence of other effects (attenuation and non-colinearity). A blur is the sinogram blurring matrix that considers non-colinearity, scattering in the detector, and DOI effect. In this section we present some of the overmentioned effects and the correction approaches that can be adopted during the image reconstruction process.

Detector geometry

The SM detector geometry correction factor is mainly related to the parallax error, which can be worsened by the scattering in the detector [START_REF] Braem | Feasibility of a novel design of high resolution parallax-free Compton enhanced PET scanner dedicated to brain research[END_REF], and it needs to be accounted for.

In the factored system model the correction is carried by the A geom matrix, which can be determined in advance through empirical [START_REF] Vv Selivanov | Detector response models for statistical iterative image reconstruction in high resolution PET[END_REF][START_REF] Zhou | Fast and efficient fully 3D PET image reconstruction using sparse system matrix factorization with GPU acceleration[END_REF], analytical [START_REF] Staelens | A three-dimensional theoretical model incorporating spatial detection uncertainty in continuous detector PET[END_REF] or Monte Carlo based studies [START_REF] Alessio | PET image reconstruction[END_REF]; or on the fly to avoid large memory storage [START_REF] Pratx | Online detector response calculations for high-resolution PET image reconstruction[END_REF][START_REF] Bert | A fast CPU/GPU ray projector for fully 3d list-mode PET reconstruction[END_REF]. The A geom matrix controls the projection operation between the image and the projection space: the forward projection is the multiplication of the image vector and the A geom matrix, while the back-projection is the multiplication of the projection vector and the A T geom (transposed matrix). This operation can be carried on through two different strategies: the LOR-driven or the voxel-driven approach. The difference consists in associating the A geom 's coefficient to a LOR or to a voxel. In LM is normally preferred the LOR-driven method since each event is associated with a LOR, for the voxeldriven approach all the events have to be read in order to determine which LORs cross the considered voxel, with an increase of the needed computation time.

Attenuation and normalization

Before reaching the detector each photon is attenuated, first in the patient body and then in every detector layer crossed. The attenuation effect needs to be considered for each LOR separately due to the different trajectories. A CT scan of the imaged object is often used to estimate the attenuation coefficients [START_REF] Paul E Kinahan | Attenuation correction for a combined 3D PET/CT scanner[END_REF].

Dual systems PET-MRI have also been developed, with the great advantage over PET-CT that the patient does not receive an additional radiation dose for the attenuation study [START_REF] Martinez-Möller | Tissue classification as a potential approach for attenuation correction in whole-body PET/MRI: evaluation with PET/CT data[END_REF][START_REF] Hofmann | MRI-based attenuation correction for wholebody PET/MRI: quantitative evaluation of segmentation-and atlasbased methods[END_REF]. The main disadvantage of using MRI images is that the conversion from gray scale to attenuation coefficients is not a straightforward operation as different materials can be seen under the same color.

In the reconstruction process, images are normalized through the sensitivity coefficients that vary for each scanner and each acquisition. Usually the normalization values are obtained before reconstruction and stored in a sensitivity map, sized as the image. A proper calculation of the sensitivity coefficients takes in consideration the attenuation in the body and it is based on a hour-long acquisition (real or simulated) of a homogeneous phantom covering entirely the FOV volume.

Positron range

Positron range correction is mandatory only for those acquisitions in which we employ a radiopharmaceutical emitting an energetic positron for which the range is relevant. In most PET imaging acquisitions, the employed radioisotope is F-18 which carries a submillimetrical correction, visible only if the detector resolution is higher, thus the A positron matrix is often equal to the identity matrix I.

There are different ways to deal with the positron range effect, some methods to reduce it through the use of magnetic fields were developed in the past [START_REF] Nl Christensen | Positron emission tomography within a magnetic field using photomultiplier tubes and lightguides[END_REF]. More recent approaches operate a correction on the projection or on the image. We can blur sinograms to simulate the average positron range [START_REF] Sf Haber | Application of mathematical removal of positron range blurring in positron emission tomography[END_REF], with the disadvantage that this method is valid only in sinogram mode. As an alternative, we can convolve (or deconvolve) the image during or after reconstruction. Different kinds of convolution can be applied: on the estimated image at each iteration, or directly on the reconstructed image. The convolution kernel needs to be studied in advance, especially for dishomogeneous objects: it strictly depends on the imaged object density, thus on the materials.

Another valid method is a Monte Carlo-based correction through the simulation of the positron propagation in the object. This method allows to reach a good correction factor, however, the use of GPUs is required [START_REF] Autret | Amelioration qualitative et quantitative de reconstruction TEP sur plate-forme graphique[END_REF].

Scattering and random coincidences

In clinical PET acquisition between 30 and 60% of the events undergo scattering before detection, consequently, the image quality is affected by a consistent amount of noise. A first approach to the correction of this phenomenon is the random coincidences estimation and subtraction. Hoffman in 1981 [START_REF] Edward | Quantitation in positron emission computed tomography: 4. Effect of accidental coincidences[END_REF] showed that the random coincidences number increases as the square of the activity: the injected dose needs to be tuned in order to avoid a high number of scattered events and to have a reasonable amount of counts.

Other methods include the use of one or multiple energy windows [START_REF] Shao | Triple energy window scatter correction technique in PET[END_REF] and the single scatter simulation by Watson [START_REF] Charles C Watson | A single scatter simulation technique for scatter correction in 3D PET[END_REF]. This approach is based on the estimation of the coincidences that scattered only once, calculated through the Klein-Nishina's formula (1.13). We can write the estimated number of coincidences detected by the AB pair that scattered only once, S AB , as:

S AB = Z V s dV s s AS s BS 4pR 2 AS R 2 BS ! µ s c ds c dW [I A + I B ], (2.41)
where V s represents the scattering volume, s the scattering position, s AS,BS the geometrical cross-sections of the two detectors and ds c dW the total cross section from the Klein-Nishina's formula.

I A,B = e AS e 0 BS e ( R A S µds+ R B S µ 0 ds) Z A,B S f ds. (2.42)
where e AS,BS is the detector efficiency for 511-keV photons and R AS,BS are respectively the distances from the detectors A and B to the scattering point S;

at last, µ is the attenuation coefficient for 511-keV photons and f the emission density.

This method known as single scatter simulation (SSS) is widely used in PET for the estimation of the scatter contribution. When using SSS the modeled correction can be scaled to match the data through tail-fitting techniques [START_REF] Ye | Scatter correction with combined single-scatter simulation and Monte Carlo simulation for 3D PET[END_REF],

although they result in artifacts in the reconstructed images when the tail is too noisy or small. The SSS method is fast and it reaches good results for homogeneous phantoms, nevertheless it does not consider multiple scattered events, and corrections based on Monte Carlo simulations have shown better results [START_REF] Sang | Ultra-fast hybrid CPU-GPU multiple scatter simulation for 3-D PET[END_REF]. Through a simulation of the photon propagation we can obtain the estimation of the scattered events and random coincidences distribution.

The main advantage of the Monte Carlo method is that all the events can be simulated in a realistic way, going from no scattering to multiple scattering;

the main disadvantage is related to the elevated computational cost in terms of time, that can be partly avoided via GPUs [START_REF] Sang | Fully 3D iterative scatter-corrected OSEM for HRRT PET using a GPU[END_REF].

Introduction

In recent times the use of imaging procedures has been growing as they develop and cover new and different diagnostics (oncological, neurological, etc.), and concerns about administered dose have been raised, especially regarding pediatric patients [START_REF] Stauss | Guidelines for 18 F-FDG PET and PET-CT imaging in paediatric oncology[END_REF]. A typical amount of total activity for a

PET study can span between 25 and 290 MBq [START_REF] Ramsey | First human imaging studies with the EX-PLORER total-body PET scanner[END_REF] for clinical scanners and between 4 and 40 MBq for small animal systems [START_REF] Simon R Cherry | Essentials of in vivo biomedical imaging[END_REF]. In any case, we deal with around 10 mSv of administered dose for a total data acquisition time between 1 and 20 minutes. A dose reduction can cause image quality to degrade [START_REF] Brix | Radiation exposure of patients undergoing wholebody dual-modality 18F-FDG PET/CT examinations[END_REF][START_REF] Xu | 200x Low-dose PET Reconstruction using Deep Learning[END_REF], thus post-acquisition image processing and software-based methods may be used [START_REF] Shtok | Sparsity-based sinogram denoising for low-dose computed tomography[END_REF][START_REF] Chen | Artifact Suppressed Dictionary Learning for Low-Dose CT Image Processing[END_REF] to recover SNR and contrast to noise ratio (CNR) in the images. New imaging devices have also been studied [START_REF] Romo-Luque | PETALO: Time-of-Flight PET with liquid xenon[END_REF], in order to tackle this same problem, especially detectors and systems taking advantage of three-gamma detection [START_REF] Lang | Sub-millimeter nuclear medical imaging with high sensitivity in positron emission tomography using b + g coincidences[END_REF][START_REF] Masełek | Towards 2+ 1 photon tomography: Energy-based selection of two 511 keV photons and a prompt photon with the J-PET scanner[END_REF][START_REF]J-PET: A New Technology for the Whole-body PET Imaging[END_REF][START_REF] Krzemie Ń | J-PET analysis framework for the prototype TOF-PET detector[END_REF].

The notion of triple coincidence was first introduced by Liang in 1987 [START_REF] Liang | Triple gamma coincidence tomographic imaging without image processing[END_REF] and resumed in 2001 by Kurfess [START_REF] James | Coincident Compton nuclear medical imager[END_REF] for a solid-detector study. The J-PET [START_REF] Moskal | Feasibility study of the positronium imaging with the J-PET tomograph[END_REF][START_REF] Kacperski | Three-gamma annihilation imaging in positron emission tomography[END_REF] project is based on the use of the three 511-keV photons originated from the ortho-positronium annihilation to reconstruct the time and position of the annihilation point on an event-by-event basis, using a TOF framework and reaching an average time resolution of around 40 ps.

In 2006 the team of Thers [START_REF] Grignon | Étude et développement d'un télescope Compton au xénon liquide dédié à l'imagerie médicale fonctionnelle[END_REF][START_REF] Thers | Nuclear medical imaging using beta (+) gamma coincidences from Sc-44 radio-nuclide with liquid xenon as detection medium[END_REF] proposed a monolithic LXe detector used as a Compton camera. A first research project called xenon medical imaging system (XEMIS) was launched to study the feasibility of the novel LXe technology to detect and reconstruct a three-gamma signal generated from Sc-44, a positron and prompt gamma emitter. A small-dimension single-phase LXe time-projection chamber (TPC) known as xenon medical imaging system 1 (XEMIS1) was developed and tested [START_REF] Gallego | XEMIS: A liquid xenon detector for medical imaging[END_REF] as a prototype. Promising results from this preliminary study brought the project to the actual preclinical state with the XEMIS2 detector, a small animal LXe detector for 3-g data acquisition,

whose goal is to exploit the intersection of the LOR with the third-gamma Compton cone to narrow the annihilation position on the LOR. In recent years, other research teams started exploring the subject of 3-g direct reconstruction,

in which through the LCI we can obtain directly the annihilation position, with no need for tomographic reconstruction [START_REF] Yoshida | Whole gamma imaging: a new concept of PET combined with Compton imaging[END_REF].

In this chapter we explore the proposed image reconstruction technique, we begin describing the LXe small animal imaging system, XEMIS2, suitable for 3-g detection; we discuss the logic behind 3-g image reconstruction, and we show the implementation and the challenges of the new technique.

The XEMIS2 project

The XEMIS2 project concerns a new small animal imaging system, whose technology has been developing at the Subatech laboratory in Nantes, France [START_REF] Manzano | XEMIS2: A liquid xenon detector for small animal medical imaging[END_REF]. The camera uses LXe as detection mean and Sc-44, produced on-site by the ARRONAX cyclotron [START_REF] Haddad | ARRONAX, a high-energy and high-intensity cyclotron for nuclear medicine[END_REF][START_REF] Roesch | Scandium-44: benefits of a long-lived PET radionuclide available from the 44Ti/44Sc generator system[END_REF][START_REF] Eppard | Pre-Therapeutic Dosimetry Employing Scandium-44 for Radiolabeling PSMA-617[END_REF], as 3-g (b+ and g) emitter. The goal is to produce good quality images with less administered dose, using the third-gamma information to narrow the annihilation position on the LOR.

Scanner geometry

The XEMIS2 is a whole-body camera for small animals that presents a particular geometry: an axially long FOV in a monolithic detector, in order to maximize the scanner sensitivity. In the XEMIS2 system we do not deal with crystals or block detectors, all the volume surrounding the small animal is filled with the detection medium, LXe. The scanner is a cylinder composed of two TPC filled with LXe; they are placed back to back and separated by a shared cathode. The active zone of detection goes from the inner 7-cm-radius cylinder to the outer 19-cm-radius cylinder and it is 24 cm long, 12 cm for each TPC, see Fig. 3.1.

All around the detector cylinder, there are 380 PMTs used to detect the scintillation light generated during photon-LXe interaction. A homogeneous electric field is applied between the cathode and the anodes, in order to drift the ionization charges and measure them on the anodes.

Radio-isotope

The positron-and-g emitter chosen for the XEMIS2 project is Sc-44. This radioisotope is a good candidate for our study due to its 4-hour lifetime, ideal for medical applications. The most probable interaction is through b-decay, branching ratio (BR)=94.27% [START_REF] Sitarz | Radionuclide candidates for b+g coincidence PET: An overview[END_REF], with the emission of a positron and a neutrino it becomes Calcium*-44 (Ca*-44). The emitted positron travels in the body and annihilates with an electron in the surroundings, while the Ca*-44 de-excites emitting an 1157-keV photon with an isotropically distributed probability in space, see Fig. 3.4. Thus, the radiation that we can detect is composed of the two back-to-back photons generated in the annihilation, and the third gamma emitted during the de-excitation process. Due to the energy difference between the coincidence photons (511 keV) and the third gamma (1157 keV), we are able to discriminate between the interactions in the acquisition process, and we can also exclude three-gamma cases deriving from the ortho-positronium state during positron decay.

An important characteristic of Sc-44 as positron emitter is the positron range. The particle is emitted in a range of energies that goes from 0 to 1474 keV, which leads to a mean positron range of 2.4 mm [START_REF] Ferguson | Comparison of scandium-44 g with other PET radionuclides in pre-clinical PET phantom imaging[END_REF], versus the 0.5 mm of the most common PET radioisotope, F-18.

Detection

The main characteristic that makes LXe a suitable detection material for 3-g detection is the excellent scintillation property, with 68 photons/keV and 2.2 ns of fast scintillation decay time. Due to its liquid state, the density is 3.100 g/mL and characterizes the material also as a high stopping power medium.

The photon detection process occurs as follows: due to the photon-LXe interaction in the active area, scintillation light is produced [START_REF] Aprile | Liquid xenon detectors for particle physics and astrophysics[END_REF] and detected in the PMTs. This signal is used to trigger the event acquisition, giving the interaction's time and position on the z-axis [START_REF] Zhu | Scintillation Signal in XEMIS2, a Liquid Xenon Compton Camera with 3g Imaging Technique[END_REF][START_REF] Xing | XEMIS: Liquid Xenon Compton Camera for 3g Imaging[END_REF][START_REF] Aprile | Noble gas detectors[END_REF]. All the other information regarding the event, as x-y position and released energy, are gathered through the ionization signal's analysis. The ionization process produces electron-ion pairs and, to avoid recombination between pairs, an electric field is applied to drift the charges towards the segmented anodes.

The signal is recorded on the anodes and it provides the x y position and the arrival time. The z position is computed through a trivial operation:

z = v drift • (t arrival t trigger )
, for which we need to know the drift velocity and the time difference between the trigger and the arrival time. The amplitude of the signal is proportional to the amount of produced charge, thus we can obtain the deposited energy from the signal recorded on the anodes.

Preliminary studies regarding the spatial and energy resolution have been carried on for the previous prototype XEMIS1 [START_REF] Gallego | XEMIS: A liquid xenon detector for medical imaging[END_REF]. A whole calibration study of the Compton camera was performed with a Sodium-22 (Na-22) source in order to characterize the properties of the prototype. From the DOI profile of the photoelectric 511-keV events in the detector, we obtained the longitudinal resolution as the result of the Gaussian fit on the peak, with a sigma of 100 µm.

A study on the ionization signal in LXe was also carried on to obtain the energy resolution value, which is expected to be around 5 % for 511-keV photons in an electric field of 1 kV/cm. While the angular resolution evaluation obtained a value of 4 for an electric field of 0.75 kV/cm.

Proposed 3-g image reconstruction technique

The principle of 3-g image reconstruction is based on the use of the two coincidence photons to determine the LOR and on the detection of the third-gamma interactions to obtain a Compton cone to use for a better localization of the annihilation position on the line. The third photon is most likely to undergo

Compton scattering in LXe due to the higher initial energy (1.157 MeV). The aim is to identify a Compton cone to delineate the direction from where the photon was emitted. In order to do so, information about the position and the energy of the first two interactions in LXe are needed. Once the cone is built, we intersect it with the LOR obtaining the LCI coordinates. This point can be used as the center of a Gaussian probability distribution function (PDF)

that gives a non-constant probability of annihilation along the LOR in image reconstruction.

LOR/cone intersection

For an 1157-keV photon the probability of scattering in LXe is 79% in 12 cm of detection mean [1]. Due to this high probability, we can obtain a Compton cone for most of the recorded events. The Compton cone axis corresponds to the line connecting the first two detection points, A and C in Fig. 3.6, and through the Klein-Nishina's Compton scattering formula in Eq. (1.13),

knowing the energy before the interaction E 0 and the deposited energy E 1 , we can compute the cone aperture angle q C as:

cos q C = 1 mc 2 E 1 E 0 (E 0 E 1 ) . (3.1)
To identify the LCI point(s), results of the intersection of the LOR and the Compton cone, we compute the system in Eq. (3.2) considering an infinite line and cone as in Fig. 3.5 [START_REF] Guertault | Intersection of a ray and a cone[END_REF]:

O D P V Q C A FIGURE 3
.5: Vectorial diagram of the cone/LOR intersection: C and A are respectively the first and second third-gamma interaction point, q = q C is the scattering angle, the vector D represents the LOR direction and the points P and Q the two intersections between the cone and the infinite line.

<

:

P = O + tD P C ||P C|| V = cos q C . (3.2)
in which C and A are respectively the third-gamma first and second interaction point; D is the LOR vector, V is the cone normal vector and P and Q are the two LCI points. We can rewrite the system as:

<

:

P = O + tD (P C•V) 2 (P C)•(P C) V = cos 2 q C , (3.3) 8 <
:

P = O + tD (P C • V) 2 (P C) • (P C) cos 2 q C = 0, (3.4)
Replacing P = O + tD we get a quadratic function:

t 2 ((D • V) 2 cos 2 q C ) + 2t((D • V)(CO • V) D • CO cos 2 q C )+ (CO • V) 2 CO • CO cos 2 q C = 0, (3.5)
We can easily solve it, if we write it in the form at 2 + bt + c = 0:

8 > > > > > > < > > > > > > : a = (D • V) 2 cos 2 q C b = 2((D • V)(CO • V) D • CO cos 2 q C ) c = (CO • V) 2 CO • CO cos 2 q C . (3.6)
To solve the system in Eq. (3.6) we must look at the determinant first:

D = b 2 4ac (3.7)
in order to discriminate if there are solutions and how many of them:

-if D < 0, there is no intersection;

-if D > 0 there is a unique intersection point;

-if D = 0 there are two points of intersection.

The computation is carried on for infinite cone and line, thus we add the requirement that only solutions within the FOV are accepted. We also decided to discard all the events presenting a double intersection (2%), as well as events with no common point between LOR and cone (4%). These two types of event mostly occur due to Compton scattering before reaching the detector. As shown in Section 3.2.1, there are several metallic layers between the FOV and the LXe volume, that increase the probability of scattering before detection, which needs to be summed up with the already existing probability of scattering in the patient's body. This can lead to LCI miscalculation, either due to an error in the LOR or in the cone computation. In Fig. 3.6 we see how a small error on the Compton angle can derive in a much bigger uncertainty on the LOR.

In our case, where we deal with three photons, events with no scattering are rare and even small-angle deviations can be an issue in the intersection determination. In Fig. 3.7 we show only a part of the possible cases: (a) a true event, (b) a coincidence photon that undergoes scattering in the metallic layers around the FOV, and (c) a 1157-keV photon that scatters in the phantom.

Other cases can be related to the scattering in the phantom for one, two or three photons at the same time, etc. All these situations can affect the LCI correct determination, and the uncertainty on the coordinates computation needs to be included in the width of the Gaussian PDF applied on the LOR during image reconstruction.

Requirements for 3-g imaging

The main deal in 3-g imaging is the use of the information carried by the third gamma to improve the image reconstruction. In order to do so, we need to determine a Compton cone and the crossing point between this object and the coincidence LOR. The precision on the LCI coordinates depends on several factors, among which the energy and spatial resolution of the third-gamma interactions, the distance between the first two detections and the distance LOR-cone vertex. In the next paragraphs we explore some of the factors that influence the LCI calculation and especially the angular resolution, which has the most influence on the cone.

Energy resolution

In Eq. (3.1) we wrote the Compton angle formula, using two variables corresponding to the transferred energy in the scattering E 1 , and the incident photon energy E 0 . We can also write the formula as:

cos q C = 1 + m e c 2 ✓ 1 E 0 1 E 1 ◆ = 1 + m e c 2 ✓ 1 E 0 E e ◆ . ( 3.8)
where E e is the energy transferred to the ejected electron.

Any error in the energy measurement can deeply affect the Compton angle value, if we consider E 0 as a known variable and we apply error propagation on the previous equation, we obtain:

s 2 E = ✓ ∂q ∂E e s E e ◆ 2 . (3.9)
where s 2 E is the detector energy resolution, and through the error propagation formula applied on Eq. (3.8), we obtain:

∂q ∂E e = 1 sin q m e c 2 (E 0 E e ) 2 .
(3.10)

The detector energy resolution s 2 E depends on three factors:

s 2 E = s 2 LXe + s 2 el + s 2 other . (3.11)
s 2
LXe is the LXe intrinsic energy resolution due to the statistical fluctuations in the number of electron-ion pairs produced during the ionization process. The second factor s 2 el is related to the readout electronics noise, which is as low as 100 e for the XEMIS2's experimental setup [START_REF] Oger | Développement expérimental d'un télescope Compton au xenon liquide pour l'imagerie médicale fonctionnelle[END_REF]. The last term s 2 other comes from all the other contributions to the energy resolution degradation as inefficiency of the materials and so on.

Spatial resolution

To determine the importance of the spatial resolution in the cone angular resolution, we have to remember that the cone axis is defined by the line connecting the first and the second third-gamma interaction.

If the source is positioned in the origin and r 1 and r 2 are the vectors that identify the two third-gamma interactions (Fig. 3.8), we can write: Thus the spatial resolution is proportional to the distance between the two interaction points. Applying the error propagation theory we obtain:

cos q C = (r 2 r 1 ) • (r 1 ) ||r 2 r 1 || • ||r 1 || . ( 3
s 2 s = s 2 (x 1 ,y 1 ,z 1 ) + s 2 (x 2 ,y 2 ,z 2 ) . (3.13)
where s 2 (x 1 ,y 1 ,z 1 ) and s 2 (x 2 ,y 2 ,z 2 ) represent respectively the spatial resolution of the first and second interaction of the third gamma. These quantities depend on the PMTs' light collection efficiency and on the pixelated anode's dimensions. For every charge cloud reaching the anode, the detection position is recorded in the center of the pixel taken as the centroid of the distribution.

The pixel size, 3.125 ⇥ 3.125 mm 2 in XEMIS2, is determinant in the error on the spatial resolution both for the Compton cone and for the LOR coordinates. We presented the energy and spatial resolution and now we show the angular resolution of the detector, which can be seen as a combination of the two.

Angular resolution

In Fig. 3.9 we display the angular resolution for a LXe Compton camera [START_REF] Oger | Développement expérimental d'un télescope Compton au xenon liquide pour l'imagerie médicale fonctionnelle[END_REF][START_REF] Gallego | Optimization of a single-phase liquid xenon Compton camera for 3g medical imaging[END_REF]. In (a) we show it as a function of the scatter angle: we can notice that the trend is mostly dominated by the energy resolution contribution, while the spatial resolution one is almost constant. In Fig. 3.9 (b) we have on the y-axis the energy transferred to the electron in a Compton interaction, and on the

x-axis the scatter angle going from zero to 180 . The maximum transferred energy occurs when the photon is back-scattered, and the angular resolution degrades for high scatter angles. Not all scatter angles are appropriate for a valid Compton cone determination: only for angles between 10 and 60 (between 40 and 610 keV of deposited energy) we obtain an acceptable angle resolution.

Pseudo-TOF image reconstruction

The proposed 3-g image reconstruction technique is based on the use of the isotropically emitted third photon to obtain clearer information about the annihilation localization on the LOR. This is achieved through the intersection of the LOR and the third-gamma Compton cone, which helps to determine the direction from where the third gamma is coming, and to delineate a LCI point. Once we compute and obtain the LCI coordinates, we can use this information in a similar way to the time difference in TOF-PET, to modify the way the PDF is distributed on the LOR.

In conventional PET image reconstruction, the probability distribution is constant along the LOR, while in TOF-PET it depends on the detector and event's characteristics and it is applied as a Gaussian distribution on the coincidence line. The center of the PDF is the length-equivalent of the time difference between the arrival of the two coincidence photons, computed following s = c • Dt, and the FWHM is related to the time resolution of the detector.

In 3-g image reconstruction we proceed in a similar way as for TOF-PET reconstruction, hence we refer to the proposed technique as pseudo-TOF.

Although, in this case, the Gaussian PDF used for the annihilation probability is centered on the LCI coordinates and the standard deviation depends on the uncertainty of the LCI calculation (see Fig. As already discussed in the section about 3-g imaging requirements 3.3.2, the LCI can be affected by several effects: photon undergoing Compton scattering in the phantom, positron range [START_REF] Rk Batra | Difference in ranges of positrons and electrons in rare-earth metals[END_REF], detector energy and spatial resolution [148], distance between the interactions, angle between cone and LOR, etc. Thus, we ought to introduce a reconstruction parameter called pseudo-TOF standard deviation that characterizes the LCI uncertainty and represents the PDF standard deviation, including all the possible effects.

In order to measure the impact of all the bias-carrying effects on the LCI, we studied the distance d between the computed LCI and the projection on the LOR of the third-gamma emission point (Fig. 3.11) through a similar-NEMA phantom simulation. Fig. 3.12 shows the resulting distribution: the obtained curve is wide and shows a maximum on 0, thus we decided to use as reference value the root mean square (RMS), 22 mm, that translates to 70 ps in TOF units.

[mm] FIGURE 3.12: Study on the LCI coordinates' uncertainty on a NEMA-like phantom: in the histograms the distance between the computed LCI coordinates and the projection of the thirdgamma emission pointon the LOR.

The use of a fixed pseudo-TOF standard deviation in the reconstruction underestimates the uncertainty for a part of the events, thus two different image reconstructions were carried. We used the 70-ps resolution value resulting from the NEMA-distance study and a variant pseudo-TOF resolution, which was calculated event by event as the time-equivalent of the distance d in ps. By incorporating the LCI into a PET system and accounting for the relative uncertainty, we expect improvement in the reconstructed images, especially since TOF scanners have already demonstrated SNR and CNR enhancement over conventional PET systems [START_REF] William | Recent advances and future advances in time-offlight PET[END_REF][START_REF] Moses | Time of flight in PET revisited[END_REF][START_REF] Wong | Image improvement and design optimization of the time-of-flight PET[END_REF].

Simulation

The XEMIS2 system was simulated using GATE [START_REF]GATE V6: a major enhancement of the GATE simulation platform enabling modelling of CT and radiotherapy[END_REF], which is an open-source software for Monte Carlo simulations of personalized medical systems and it is mostly used to test new designs and geometries, to optimize and experiment new methods and new data acquisition techniques.

In GATE the detector geometry is defined as a multi-layer cylinder: the active zone measures 240 mm of length and respectively 70 and 190 mm of inner and outer radius. The cylinder is filled with LXe and the FOV with air.

The active volume is enclosed in a double-layer cylindrical armor of stainless steel and aluminum. The only discontinuity in the detector active zone is due to the central copper cathode that divides the two TPCs. In Fig. 3.13 an overview of the detector geometry taken from the GATE visualization tool.

Particles and matter were simulated using GEometry ANd Tracking (Geant4), particle generation based on CLHEP libraries [START_REF] Lönnblad | CLHEP-a project for designing a C++ class library for high energy physics[END_REF], while photon interaction follows the Penelope model in GATE [START_REF] Agostinelli | GEANT4: A Simulation toolkit[END_REF]. third gamma is emitted isotropically and independently of the coincidence photons' direction. Although more particles and interactions are involved in the event, we decided to consider only the detections directly related to the 3-g. Scattering is allowed in the phantom and in the detector. For each interaction we record position and deposited energy.

Data acquisition time is 20 minutes and the initial total activity 20 kBq. The events are acquired in LM to associate an event-based pseudo-TOF standard deviation. In this work we explore the pseudo-TOF image reconstruction method and its benefits and we compare the proposed reconstruction technique to conventional two-gamma (2-g) image reconstruction. Considering this, the first simulated phantom is a NEMA-like phantom, in analogy to the imagequality studies carried out to characterize detectors, and in general in TOF versus non-TOF comparisons [START_REF] Disselhorst | Image-quality assessment for several positron emitters using the NEMA NU 4-2008 standards in the Siemens Inveon small-animal PET scanner[END_REF].

Phantoms

The phantom consists of a 12-cm-long homogeneous water cylinder of 3.5 cm radius, containing five spheres. All spheres' centers are equidistant from the cylinder central axis and positioned on the central slice of the cylinder.

The radius of the spheres measure 2, 4, 8, 10, and 12 mm, see Fig. 3.14 and Table 3.1 for further details.

The total simulated activity in the phantom is 20 kBq, with a factor 15 of contrast ratio between the sphere and the background. 30 duplicates were simulated to study the image variance.

To complete the XEMIS2 characterization, as second phantom we chose to simulate a more realistic object: a 28-g male mouse known as the Digimouse [START_REF] Stout | Creating a whole body digital mouse atlas with PET, CT and cryosection images[END_REF][START_REF] Dogdas | Digimouse: a 3D whole body mouse atlas from CT and cryosection data[END_REF]. The small animal presents a complete structure, composed of brain, muscles, eyes, glands, heart, lungs, liver, stomach, spleen, pancreas, kidneys, testes, bladder, skeleton, and skin. A 2-mm radius spherical tumor was added in the brain as the region of interest (ROI) to analyze in the image quality study. In Fig. 3.15 a transaxial view of the Digimouse central slice and in Table 3.2 the simulated uptake distribution of the mouse.

The radiopharmaceutical employed in this simulation is [ 44 Sc]Sc-J591, labeled with Sc-44, and it is an anti-prostate-specific membrane antigen monoclonal antibody, whose biological uptake is known from Carter [START_REF] Lukas | The impact of positron range on PET resolution, evaluated with phantoms and PHITS Monte Carlo simulations for conventional and non-conventional radionuclides[END_REF] and Holland [START_REF] Holland | 89Zr-DFO-J591 for immunoPET of prostatespecific membrane antigen expression in vivo[END_REF]'s studies. One of the advantages of using [ 44 Sc]Sc-J591 is that the molecule is also used for positron range studies since it is easily associated also to F-18.

The total activity simulated in the phantom was around 22 kBq and data acquisition time was 20 minutes. Due to the more detailed GATE-voxelized structure of the Digimouse phantom, longer computation time was needed and only one acquisition was simulated.

Image reconstruction software 3.5.1 Pseudo-TOF system matrix

Before getting into the details on the software implementation of 3-g image reconstruction, let us introduce the SM for pseudo-TOF. SM elements are often computed event-by-event, despite it being computationally intensive due to the matrix size. They depend on the physics and on the system, and they are characteristics of each detector and each data acquisition, thus they need to be computed for each event.

In TOF-PET the SM gives the probability p i,j,b that the annihilation of a positron emitted in the j-th voxel is detected in the i-th LOR's bin and in the time bin b, In the pseudo-TOF context, we deal with position and distance instead of time resolution, thus b cannot be considered the number of possible time bins. We translate it in the number of possible pseudo-time bins, each of them representing a different LCI position on the LOR i. Consequently, the coefficient p i,j,b becomes the probability that an annihilation that occurred in the j-th voxel was detected on the i-th LOR and that the third gamma was emitted from the pseudo-time bin b on the LOR. In absence of scattering and attenuation, this probability depends on the detector sensitivity, on the spatial and energy resolution, and on the Sc-44 positron range. To calculate the SM elements in pseudo-TOF, we can follow the same TOF method and use the pseudo-TOF resolution instead.

with i 2 [1, • • • , M] (I number of possible LORs), j 2 [1, • • • , J] (J

Pseudo-TOF algorithm in CASToR

All the images we reconstructed in this work are obtained through the LM-MLEM reconstruction algorithm in CASToR v2.3 [START_REF] Merlin | CASToR: a generic data organization and processing code framework for multi-modal and multi-dimensional tomographic reconstruction[END_REF], both for pseudo-TOF and non-TOF events, which we refer to as 3-g and 2-g respectively.

As for all the iterative optimization algorithms as MLEM, the reconstruction in CASToR is based on a projector handling both back and forward projection and there is an optimizer and a convolver, in case any kind of convolution is required (e.g. for positron range correction).

For each iteration, there are several steps to follow: first, the software calls the projector to compute one row of the SM. Then the optimizer, composed of both the objective function and the iterative optimization, performs the data update in several operations. It starts applying a forward projection on the image estimate, then it adds all the provided estimation of the background, and in the end it computes the correction term in the data space, in order to back-project them and obtain a correction image for the considered SM row.

In CASToR the computation of the system matrix elements follows the Siddon projector [START_REF] Robert L Siddon | Fast calculation of the exact radiological path for a three-dimensional CT array[END_REF][START_REF] Jacobs | A fast algorithm to calculate the exact radiological path through a pixel or voxel space[END_REF] and the so called ray-tracer method, for which the probability depends on the portion of LOR crossing each voxel. The standard optimizer used for two-gamma events is from Shepp and Vardi [START_REF] Shepp | Maximum Likelihood Reconstruction for Emission Tomography[END_REF]; while for the LM-MLEM TOF algorithm we refer to Filipovic [START_REF] Filipovic | Time-of-flight implementation for PET reconstruction in practice[END_REF].

The pseudo-TOF LM-MLEM algorithm can be expressed as: elements. In this instance, the total virtual elements' number is perfectly reasonable.

l (k+1) j = l (k) j  i,b p i,j,b  i p i,j,b  j 0 p i,j 0 ,b l (k) j 0 . ( 3
To identify the couple of elements involved in the event, we added an intermediate step between data acquisition and image reconstruction, in which we transfer the coordinates of the detected interactions on coordinates that lay on the discretized part of the detector. We intersect the event LOR with the FOV cylinder and once the couple of virtual elements are pointed out by the intersection (red elements in Fig. 3.17 (c)), we pass their ID over to CASToR to store it and subsequently use it during reconstruction. Another information needed for each event is the LOR length and the LCI position to be translated in time-units (ps).

Sensitivity image

When reconstructing in CASToR, the sensitivity image can be given as an input image file before launching the reconstruction or, if the file is missing, it is computed by the software before starting the iterations. CASToR computes the sensitivity map running a loop over all possible LORs, considering all the couples of elements in the detector geometry. In this way the sensitivity image depends only on the geometrical characteristics of the detector, attenuation and normalization are not taken into account.

In our case, due to the fact that we modified the XEMIS2 geometry discretizating the detector, and we removed the concept of depth, we provided one of our own. We obtained it simulating a homogeneous water cylinder filled with Sc-44, as large as the FOV. The hours-long simulation brought us to collect millions of events. We then reconstructed the events through an MLEM sensitivity optimizer algorithm and used the first iteration output image as our sensitivity map.

In this case, since the cylinder was filled with water and scattering was allowed in the phantom, we can consider that the scattering and attenuation effects in water are included in the image.

Figures of merit

The Gaussian smoothing filter (s = 1) was applied on all the reconstructed images before computing the metrics.

NEMA-like phantom analysis

In the NEMA phantom analysis parameters were computed throughout all the iterations for every noise replicate and the variance was calculated as follows:

s 2 c = 1 |R c | Â j2R c ✓ 1 L L Â `=1 (x j µ c j ) 2 ◆ , (4.1)
where:

µ c = 1 30|R c | L Â `=1 x c. (4.2)
c 2 {b, s} is the label for the background cylinder and for the hot spot, and

x c is the reconstructed activity from the `th simulation; a total of L = 30 repetitions were performed. The variance was used to study the trend of the image quality parameters with respect to the general noise.

The first studied parameter is the SNR, determined as the difference of activity between the sphere and the background, divided by the standard deviation in the cylinder s b :

SNR = µ s µ b s b , (4.3)
where the variables were calculated following equations (4.1) and (4.2).

For the contrast analysis we computed CNR and RC: the first as the absolute difference between the mean value in the sphere and in the cylinder, divided by the sum of the two standard deviations:

CNR = |µ s µ b | q s 2 b + s 2 s , (4.4)
while the RC was defined as the fraction of the difference between the mean value in the sphere and in the cylinder for the reconstructed images over the same quantity related to the ground truth (GT) image:

RC = µ s µ b µ GT s µ GT b . ( 4.5)
To obtain a GT image for each repetition, we mapped back to an image all the annihilation positions from the GATE simulation.

Digimouse phantom analysis

For the Digimouse phantom simulation, we worked only on one data acquisition, thus the figures of merit differ from the previous analysis.

We consider the added sphere in the brain as the volume of interest (VOI), with µ s and s s being respectively the mean value and the standard deviation in the hot spot. In this case the computation of the respective background variables is more complex due to the non-homogeneous background. We

SNR digimouse = µ s µ b s b , (4.6
)

CNR digimouse = |µ s µ b | q s 2 b + s 2 s . (4.7)
where:

µ b = 1 M M Â i=1 µ i , (4.8)
and

s b = 1 M M Â i=1 s i . (4.9)
in which M = 4 are the four background spheres in the Digimouse.

For implementation reasons, the simulation of a GATE voxelized phantom did not allow direct information about the local interactions in the mouse, and we did not reconstruct a GT image starting from the annihilation positions.

As a consequence, no RC study was carried out in this case.

NEMA-like phantom results

In this section we show results related to the NEMA-like phantom obtained through 3-g and 2-g MLEM reconstruction. For the pseudo-TOF reconstruction, we used two different pseudo-TOF approaches: in one reconstruction we used a fixed Gaussian standard deviation of 70 ps, which is the value corresponding to the average uncertainty on the LCI calculation. In the second approach the resolution was event-based, hence variable, computed for each event as the time-equivalent of the difference between the LCI and the annihilation point projection on the LOR, as shown in Fig. 3.12.

The reason for showing results with the average pseudo-time resolution and with the variable resolution is to study the potential benefits of the proposed approach in favorable and less favorable scenarios of reconstructed cone accuracy.

Digimouse phantom results

First difference between the NEMA-like and the Digimouse study is that in the second case images were obtained only through conventional 2-g reconstruction and pseudo-TOF reconstruction with several resolution values (70 ps, 100 ps, 150 ps, and 200 ps). Moreover, in this case events with and We also computed the profile along the ROI in the Digimouse's brain, which can be seen in Fig. 4.7 for all the different reconstructions.

Discussion

In the previous section 4.3 and 4.4 we presented the results for both phantoms using 3-g and 2-g image reconstruction algorithms, that we are going to discuss in this section.

We can observe the NEMA-like phantom images in Fig. 4.2, and we can notice that the smaller ROI is not easily distinguishable in the 2-g reconstruction (a), and due to this visible flaw in the image, we do not consider the results regarding the first sphere to be particularly reliable. For all the other ROIs a contrast improvement is visible in the 3-g cases, enhanced by the plot profiles We did not observe the same improvement as for the NEMA-like analysis in the Digimouse study, neither visually nor quantitatively in the figures of merit 4.6. A higher statistics is needed for such a small and detailed phantom.

As we fully complied with the regular amount of activity injected in a small animal as such [START_REF] Dogdas | Digimouse: a 3D whole body mouse atlas from CT and cryosection data[END_REF], the only way to increment the statistics was to include the 2-g events, that did not allow to evaluate the enhancement brought by the 3-g technique. In the Digimouse case we faced a prohibitive low 3-g statistics that required to include the 2-g events reconstructed following conventional 2-g MLEM. Contrarily, the NEMA case allowed for the reconstruction using only 3-g events that enhanced the 3-g reconstructed volumetric images both visually and quantitatively.

Introduction

There are several factors degrading image resolution in PET imaging: some are due to the detector's characteristics (sensitivity, scanner geometry, spatial and energy resolution), some depend on the physical interactions involved (positron range [START_REF] Blanco | Positron range effects on the spatial resolution of RPC-PET[END_REF], Compton scattering, photon attenuation). For most of these effects correction methods have been studied in PET reconstruction, some of which are very successful (e.g. for attenuation calculation or for random coincidence).

In the positron range case, the modelization is often simplified: we consider a homogeneous medium or in most of cases, if the radio-pharmaceutical presents a submillimetric positron range, the correction is not even applied, as for F-18 imaging.

Reconstructing without correcting for the positron range means obtaining the positron annihilation distribution and not directly the positron emission distribution, which corresponds to the radio-pharmaceutical uptake in the body. As already introduced in Chapter 1, the positron range increases with the positron initial kinetic energy and decreases with the material density, thus it is an important effect in the case of energetic positrons in soft tissue. In The simulation of a positron interacting in the matter is not a trivial problem: in fact, as we can see in Fig. 5.1, in the cross-section calculation we need to sum several types of interaction whose importance is energy-dependent [START_REF] Nakamura | Review of particle physics[END_REF][START_REF] Rohrlich | Positron-Electron Differences in Energy Loss and Multiple Scattering[END_REF]. Scattering for a positron is considered as ionization if the energy loss is below 0.255 MeV, while it is referred to as Bhabha scattering (electron-positron diffusion) when the energy loss is higher.

One of the difficulties in simulating positron interactions is due to the amount of computation time needed to take into account every step. A positron tends to undergo scattering until it loses all its energy and is able to annihilate and the energy loss process is long and composed of very little steps. For each of them, we have a different energy loss, due to its energy dependence, thus a thick sampling would be needed to simulate a close approximation of the real interactions.

Positron range correction

In this section, we take a look at the different approaches to deal with the positron range correction. There exist two main techniques to minimize the effect: the first is to reduce the particle range, while the second is to correct it working either on the projections, either on the reconstructed image, either on the estimated image iteration per iteration.

Positron range reduction

Positron are charged particle, subdued to the Lorentz force in presence of a magnetic field:

F = qE + qv ⇥ B. (5.1)
where F is the Lorentz force experienced by the particle, q the particle charge, v the particle velocity, E the electric field and B the magnetic field.

This force is applied perpendicularly to the magnetic field vector B and to the velocity vector v and it reaches its maximal intensity when the positron travels orthogonally to the magnetic field.

Several studies were carried out on the subject [START_REF] Wirrwar | 4.5 tesla magnetic field reduces range of highenergy positrons-potential implications for positron emission tomography[END_REF][START_REF] Raymond R Raylman | Combined MRI-PET scanner: a Monte Carlo evaluation of the improvements in PET resolution due to the effects of a static homogeneous magnetic field[END_REF][START_REF] Bertolli | PET iterative reconstruction incorporating an efficient positron range correction method[END_REF][START_REF] Burdette | A study on PET image quality using both strong magnetic fields and a ML-EM positron range correction algorithm[END_REF][START_REF] Huang | The effect of magnetic field on positron range and spatial resolution in an integrated whole-body time-of-flight PET/MRI system[END_REF],

showing that the particle range is enlarged and its direction changes to follow the Lorentz force as expected. Nevertheless, the magnetic field needs to be high to have a visible effect (> 3 T), which carries concerns on the cost, considering also that ferromagnetic materials cannot be employed in the construction of this type of PET system. In addition, the positron is affected only in one of the three directions, the one perpendicular to the magnetic field, ergo its range is reduced only in one of the three axis.

Positron range correction pre-, post-and during reconstruction

In order to apply positron range correction in distinct reconstruction steps (pre-, post-and during), we need to deal with different data types.

In the pre-reconstruction correction case, data acquisition ought to be in sinogram mode so that we can work on denoising the projections [START_REF] Stephen E Derenzo | Mathematical removal of positron range blurring in high resolution tomography[END_REF][START_REF] Andrew | EM algorithm system modeling by imagespace techniques for PET reconstruction[END_REF].

Haber and Derenzo in 1990 [START_REF] Sf Haber | Application of mathematical removal of positron range blurring in positron emission tomography[END_REF], proposed a technique for which a Fourier deconvolution was applied on the projections to remove range blurring for 68 Ga and 82 Rb acquisitions. Since the point spread function (PSF) used in the deconvolution was calculated in a homogeneous mean (water) and considered uniform in the FOV, the whole process was valid only for uniform phantoms.

Another disadvantage of this technique is that the deconvolution was applied by taking the Fourier transform of the positron range function and dividing it into the transform of the data: the division in the frequence space increases the noise level damaging the overall image quality.

A similar technique can be applied to the reconstructed image which corresponds to the annihilation distribution, when not corrected for the positron range. We can think of deconvolving the distribution with a voxel-dependent PSF in order to obtain a blurred positron emission map. This technique is not often of use due to the level of noise that we reach in the image after convolution.

The main currently used approach is positron range correction during the reconstruction process. This consists of estimating and adding the positron range effect in the system matrix, which is composed as a product of matrices [START_REF] Michel | Iterative image reconstruction for positron emission tomography based on a detector response function estimated from point source measurements[END_REF][START_REF] Vladimir Y Panin | Fully 3-D PET reconstruction with system matrix derived from point source measurements[END_REF], as we have already introduced in Section 2.6.3. In this way the positron range matrix can be used only during the projection operation and not the back-projection, reducing the computational time and the number of iterations needed for convergence [START_REF] Cal-González | Study of CT-based positron range correction in high resolution 3D PET imaging[END_REF]. In this case we talk about unmatched projector. To correct for the positron range effect we apply a convolution with a voxel-dependent kernel all along the FOV. In fact, as we have already seen, the positron range depends on the material density and on the positron energy. The kernels need to be studied ad hoc based on the phantom and on the employed radioisotope, which brings us to the main problem of this correction method, the memory storage. To simplify calculations and to avoid a huge matrix's memory allocation, it is often suggested to use an invariant and isotropic kernel [START_REF] Rahmim | Resolution modeled PET image reconstruction incorporating space-variance of positron range: Rubidium-82 cardiac PET imaging[END_REF] that is computed in an homogeneous phantom, often made of water. This might result in an over or under correction of the images if the positron emission happened in a non-water equivalent tissue (e.g. lungs [START_REF] Emond | Effect of positron range on PET quantification in diseased and normal lungs[END_REF] or bone). The best results are achieved when the imaged objects do not present high heterogeneity and the tumors are not on the boundary between different materials.

Several studies for positron range correction based on variant kernel convolution on anatomical MRI and CT images were also conducted [START_REF] Kemerink | Effect of the positron range of 18 F, 68 Ga and 124 I on PET/CT in lung-equivalent materials[END_REF][START_REF] Cal-González | Positron range effects in high resolution 3D PET imaging[END_REF][START_REF] Kraus | Simulation study of tissue-specific positron range correction for the new biograph mMR whole-body PET/MR system[END_REF][START_REF] Alessio | Spatially variant positron range modeling derived from CT for PET image reconstruction[END_REF][START_REF] Jacobo | Positron range and prompt gamma modeling in PET imaging[END_REF]. Other techniques were studied for heterogeneous objects [START_REF] Wan | Probability based positron range modeling in inhomogeneous medium for PET[END_REF][START_REF] Bai | Positron range modeling for statistical PET image reconstruction[END_REF]; Bai in 2005 [START_REF] Bai | Evaluation of MAP image reconstruction with positron range modeling for 3D PET[END_REF] proposed the truncation of the PSF to reduce artifacts due to the boundary effect. Results showed that at matched noise levels this positron range model gives better contrast, nevertheless the best trade-off between noise and resolution is achieved without truncation.

The approach of Alessio and MacDonald in 2008 [START_REF] Alessio | Spatially variant positron range modeling derived from CT for PET image reconstruction[END_REF] consisted of estimating variant kernels for all materials in a homogeneous mean and build a voxel-dependent kernel during reconstruction, interpolating the previously estimated kernels.

Monte Carlo simulations are the alternative to analytical models, although they highly extend the reconstruction time to propagate the positron interactions iteration per iteration on forward and back-projection.

Proposed correction

For the XEMIS2 project the interest in positron range correction concerns the use of Sc-44 as 3-g radioisotope, whose positron range is not negligible as seen in Table 1.1 in Chapter 1.

In this section we discuss the two implemented approaches in CASToR to correct for the effect: first through an invariant kernel approach and second through a material-based variant kernel. The materials map was obtained from the GATE phantom's file.

Simulation and phantoms

The phantom chosen for this study is composed of a 12-cm-long cylinder with 3.5-cm radius positioned in the center of the FOV, within which we find five hot spheres, see Fig. 5.3 and Table 5.2. The materials were chosen to be water, lung and bone, to explore three mediums of very different densities (1.00 g/cm 3 water, 0.26 g/cm 3 lung and 1.42 g/cm 3 bone) [START_REF] Emond | Effect of positron range on PET quantification in diseased and normal lungs[END_REF][START_REF] Cal-González | Tissue-dependent and spatially-variant positron range correction in 3D PET[END_REF].

Data are simulated in GATE v6 using the Penelope model. The total simulated activity is 20 kBq for a data acquisition time of 20 minutes, with a factor 15 of contrast between the hot spots and the background.

Reconstruction and positron range correction estimation

The reconstruction was carried on in CASToR using a LM-MLEM algorithm in which we implemented the positron range correction through convolution.

We show results obtained with conventional 2-g and 3-g reconstruction to be able to compare the two approaches, with and without positron range correction.

CASToR is designed as a flexible and generic software in which we can easily develop an image convolution in any algorithm, at any point during the program execution. In particular, for a stationary kernel, the convolution is already implemented in the software, with multiple choices as of when to apply it: on the image to be forward projected, on the back-projection's correction terms, on the reconstructed image as a post-processing step, or on the current estimated image. Moreover, CASToR's convolution operation is fast since the image is saved in a padded buffer, in which zeros are added along each dimension following the maximum kernel size on the current Material FWHM [mm] FOV position [mm] Water 1.9 0, 0, 0 Bone 1.9 0, 0, -60 Lung 2.2 0, 0, 60 dimension.

In our implemented correction techniques with a variant kernel, we chose to apply an isotropic convolution on the image to be forward projected at each iteration. The developed algorithm accepts 4 different kernels, for the 4 materials in the phantom, and for each of them, we can specify the number of Gaussian sigmas and the FWHM (axial and transaxial). In order to associate the right kernel to every voxel, the convolution algorithm takes in input a labeled image with numbers from 0 to 3, each of them corresponding to a different kernel to be used during convolution. In our case the labeled image was directly obtained from the GATE phantom and based on the simulated materials.

Results

We begin this section showing the metrics used in the image analysis, then we present the results for both invariant and variant kernel convolution.

Figures of Merit

The

Discussion

The results in this section are obtained applying the positron range correction to the image before each iteration forward projection. We present images and plots concerning 2-g and 3-g reconstruction with variant and invariant kernel convolution. The variant kernel correction was based on a material study, carried out through three PSFs' simulation related to the three materials found in the phantom.

A first effect of positron range correction can be seen in Fig. 5.8, where the reconstructed images using conventional PET reconstruction (first column) and our proposed 3-g approach (second column) are compared to the ones obtained through the same image reconstruction techniques with the positron range convolution correction.

Apart from the contrast improvement due to the use of 3-g instead of classical PET reconstruction, the general contrast of all the five spheres is improved in the corrected images, especially for 3-g in (d) and (f), respectively variant and invariant kernel convolution. The improvement is properly quantified in the plot in Fig. 5.9, where we show the profile along the yellow lines for all the six reconstructions. Between the profiles of the 2-g reconstructions and the ones related to the 3-g we have a gap of 20% on the grayscale, while there is no remarkable difference between convolved images and not corrected ones for the same image reconstruction technique. A proper comparison cannot be carried on between two methods so different from each other, since a very important feature in image reconstruction is the convergence rate. Indeed, the noise in the image increases with the iteration number, as well as the time needed to obtain results. Those two characteristics are of the highest importance in a medical context and cannot be ignored.

In the shown results we see that the SNR and CNR peak values are much higher for the 3-g technique than for the conventional 2-g reconstruction, nevertheless, the values at the 30-th iteration tend to keep up only for notcorrected 3-g images, due to the blurring introduced by the convolution.

In this particular case, the best way to acknowledge an improvement once the correction is applied lies in the visual results in Fig. 5.8.

three clusters of interaction in LXe : the two clusters corresponding to the annihilation photons and the one corresponding to the third gamma.

The main challenges related to the detection of three gamma events in real data concern the clustering of the LOR and 3-g groups, as well as the determination of the Compton angle and of the first two points of interactions for the prompt gamma used to determine the cone.

For the first issue, several clustering methods have been studied through GATE simulation. An initial approach for XEMIS2 proposed by N. Beaupère There are, therefore, N! possible combinations of interactions to consider while looking for the optimal path. The selected sequence of interactions is chosen using a c 2 test on the difference between the geometrical angle and the one obtained through Compton kinematics.

The main issues with this method are the low number of events registered as 3-g due to the many sources of error in the LCI computation process, and the long computational time associated with the determination of the right path-sequence for each event. A novel approach to improve the determination of the LCI coordinate using machine learning methods is exposed in the next section.

3-g detection using Neural Networks

Several works have proposed to use neural networks for the determination of the interaction point in monolithic detectors [START_REF] Iborra | Ensemble of neural networks for 3D position estimation in monolithic PET detectors[END_REF][START_REF] Wang | 3D position estimation using an artificial neural network for a continuous scintillator PET detector[END_REF][START_REF] Bruyndonckx | Evaluation of machine learning algorithms for localization of photons in undivided scintillator blocks for PET detectors[END_REF][START_REF] Bruyndonckx | Neural network-based position estimators for PET detectors using monolithic LSO blocks[END_REF][START_REF] Zoglauer | Application of neural networks to the identification of the compton interaction sequence in compton imagers[END_REF]and some studies have also focused on establishing Compton kinematics using machine learning methods [START_REF] Peng | Compton PET: A simulation study for a PET module with novel geometry and machine learning for position decoding[END_REF][START_REF] Michaud | Sensitivity in PET: Neural networks as an alternative to compton photons LOR analysis[END_REF]. These techniques are mainly based on 511-keV photons scattering as they were proposed for PET but can nevertheless be considered in the present context of 3-g imaging.

Recent studies carried on in our group have demonstrated that an ensemble of deep multi-layer neural networks trained on extensive Monte Carlo simulations is able to robustly predict the interaction point in monolithic detectors even in the presence of imperfect geometries [START_REF] Iborra | Ensemble of neural networks for 3D position estimation in monolithic PET detectors[END_REF].

We suggest to leverage the predictive power of similar neural network architectures to better determine the LCI point in future XEMIS2 experiments.

The LCI determination could be divided into two sub-problems related to:

1. The clustering of the three groups of interaction corresponding to a 3-g event.

2. The reconstruction of the Compton cone from the third gamma group.

These two problems can be tackled using different feed-forward multilayer perceptrons (MLPs) architectures.

For the clustering we can consider a 4-class MLP classifier taking as inputs the positions and energies of each interaction within an event and predicting as output the classification in four groups: two annihilation groups, one third-gamma group and one group corresponding to spurious interactions.

For the prediction of the Compton cone, two strategies may be adopted:

• Similarly to the clustering case, a MLP classifier can be trained to identify the first and second interaction within the third gamma group, followed by conventional Compton kinematics as described in this thesis.

• Alternatively, a regressor MLP could be directly trained to predict the LCI using the third gamma without explicit Compton kinematics. The criterion would be the minimization of the distance to the annihilation point through e.g. least squares.

While both classification and regression strategies can be considered, classification MLPs seem more appealing as they could directly serve as a replacement for the current clustering approach. The impact of such a procedure could be evaluated straightforwardly using the same reconstruction method developed in this thesis. In the near future, we will therefore study the potential improvements achieved by this strategy on image reconstruction quality.

Contributions

In this manuscript we presented a new image reconstruction technique for 3-g imaging, that we called pseudo-TOF.

This thesis' work consisted in:

-the implementation of the proposed technique within an existing image reconstruction framework (CASToR);

-a GATE simulation study using Sc-44 and XEMIS2 geometry using two phantoms: a NEMA-like and a more realistic Digimouse;

-a quantitative analysis with respect to the pseudo-TOF resolution as reconstruction parameter;

-the implementation of positron range correction through variant and invariant kernel convolution in CASToR;

I 2 . 6 4 . 7

 2647 photon scatters with a small angle before reaching the liquid xenon (LXe), causing a change in the LOR/cone intersection (LCI) coordinates. (c) Event with an 1157-keV gamma scatter in the phantom: depending on the angle width the cone can or not intersect the LOR, anyways the LCI determination is highly affected. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Scheme of the vectors involved in the spatial resolution calculation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5. 8

 8 Transaxial view of the NEMA-like phantom, reconstructed with four different techniques: (a) conventional 2-g PET and (b) 3-g, with no positron range correction; (c) 2-g and (d) 3-g reconstruction convolved with a variant kernel, (e) 2-g and (f) 3-g reconstruction convolved with a stationary kernel. . . . . 5.9 Plot profiles of the six different reconstructions through the yellow line seen in Fig. 5.8. . . . . . . . . . . . . . . . . . . . . . 5.10 SNR and CNR plots of the three smallest spheres, water-made:

5. 11 44 CASToRF- 18 Fluor- 18 FBPXEMIS2 xenon medical imaging system 2 A

 114418182 photons dans les détecteurs. Pour maintenir des normes de qualité d'image requises par la médecine moderne tout en diminuant l'activité totale, un certain nombre de méthodes software de traitement d'image ont été mises en oeuvre en parallèe de l'utilisation de nouveaux systèmes d'imagerie. De nouveaux systèmes, appelés TOF-PET (imagerie temps de vol), ont d'abord été conçus dans les années 1970. Leur objectif est de mesurer la différence de temps entre les deux arrivées de photons par coïncidence dans des détecteurs opposés dans l'espace pour estimer la position le long de la ligne de réponse (LOR) où l'annihilation a eu lieu. L'incertitude de localisation est déterminée par la résolution temporelle du détecteur, qui peut être caractérisée par la largeur à mi-hauteur de la distribution centrée sur x = c • Dt, où c est la vitesse de la lumière et Dt la différence de temps. En utilisant un détecteur idéal avec une résolution temporelle infiniment précise, on pourrait localiser exactement la position d'annihilation et donc avoir directement accès à la biodistribution des radiotraceurs. Cependant, en raison des limites de techniques d'instrumentation actuelles, les techniques tomographiques sont toujours nécessaires. Les algorithmes utilisés pour reconstruire les images tomographiques peuvent être classés en deux groupes principaux: les méthodes analytiques et itératives. Dans le premier cas, nous traitons un problème d'inversion entre une fonction et son intégrale curviligne. Dans dans le second cas, nous modélisons tous les aspects possibles de l'acquisition de données et toutes les interactions physiques qui peuvent se produire et nous cherchons l'image qui correspond le mieux aux données de manière itérative. Les données sur lesquelles nous basons notre la reconstruction sont acquises en détectant des photons de coïncidence de 511 keV en TEP et photons uniques en SPECT. Dans les deux cas, le nombre limité d'événements et le comportement stochastique des particules entraîne l'apparition de

FIGURE 1 :

 1 FIGURE 1: Geometrie du XEMIS2.

FIGURE 2 :

 2 FIGURE 2: Vision transversale du scanner.

FIGURE 3 :FIGURE 4 :FIGURE 5 :

 345 FIGURE 3: Schéma de la désintégration b du Sc-44 [1].

FIGURE 8 :

 8 FIGURE 8: Coupe axiale du fantôme type NEMA. Radius [mm] Position (X,Y,Z) [mm] 2 (0.0, 17.5, 0.0) 4 (-8.8, 15.1, 0.0) 8 (-17.5, 0.0, 0.0) 10 (0.0, -17.5, 0.0) 12 (17.2, 3.7, 0.0)

Fig. 7 .

 7 Fig. 7. La raison de montrer les résultats avec la résolution pseudo-temporelle moyenne et avec la résolution variable est d'étudier les bénéfices potentiels de l'approche proposée dans des scénarios favorables et moins favorables de précision du cône reconstruit.

FIGURE 10 :FIGURE 11 :

 1011 FIGURE 10: Coupe axiale du fantôme type NEMA après avoir applique un filtre gaussien (s = 1 mm). Images reconstruite utilisant (a) reconstruction TEP classique, et la technique proposée dans ce travail (3-g) avec (b) 70-ps et (c) un écart type variable événement par événement. Les images sont les résultats a la 30 th itération. (d), (e) et (f) montrent les profiles sur la ligne jaune.

FIGURE 12 :

 12 FIGURE 12: Profile sur la ligne qui croise la ROI dans le cerveau du Digimouse.

FIGURE 13 :

 13 FIGURE 13: Coupe axiale du fantôme pour l'étude sur la correction de range du positon.

ResumeFIGURE 14 :

 14 FIGURE 14: Coupe transversale du fantome. Images reconstruites avec differents techniques: (a) TEP classique et (b) 3-g,pas corrigee; (c) 2-g et (d) 3-g convolues avec des noyaux variantes, (e) 2-g et (f) 3-g convolues avec un noyau stationaire.

FIGURE 15 :

 15 FIGURE 15: Profile sur la ligne jaune.

  Nous avons montré une étude basée sur la simulation d'une acquisition de données de 20 minutes avec 20 kBq de radioactivité initiale dans le détecteur XEMIS2. Nous avons simulé un fantôme de type NEMA et Digimouse et montré une amélioration de la qualité d'image pour des objets de différentes densités et dimensions. Nous avons présenté le calcul des coordonnées LCI, sur lequel nous basons la méthode de reconstruction d'image 3-g. Nous avons discuté des caractéristiques obligatoires du détecteur en termes de résolution spatiale, énergétique et angulaire pour réduire l'incertitude sur le point LCI. Le calcul des coordonnées LCI est une étape cruciale vers l'application de l'approche proposée dans un contexte clinique réel. Pour cela, il faut être capable de déterminer avec suffisamment de précision la localisation de trois clusters d'interaction dans LXe: les deux clusters correspondant aux photons d'annihilation et celui correspondant au troisième gamma. Les principaux défis liés à la détection de trois événements gamma dans des données réelles concernent le clustering des groupes d'interactions, ainsi que la détermination de l'angle de Compton et des deux premiers points d'interactions pour le troisième gamma utilisé pour déterminer le cône.

  coincidence photons' arrival in opposite detectors to estimate the position along the line of response (LOR) where the annihilation took place. The localization uncertainty is determined by the time resolution of the detector, which can be characterized by the full width half maximum (FWHM) of the distribution centered in x = c • Dt, where c is the speed of light and Dt the time difference. Using an ideal detector with infinitely accurate time resolution, one could exactly locate the annihilation position and therefore have direct access to the radiotracer biodistribution. However, due to the limitations of current instrumentation techniques, tomographic techniques are still required.

Contents 1 . 1 24 1. 2

 11242 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . Physical principles . . . . . . . . . . . . . . . . . . . . . . . . 24 1.2.1 Radioactive decay . . . . . . . . . . . . . . . . . . . . 24 1.2.2 Charged particle's energy loss . . . . . . . . . . . . . 26 1.2.3 Positron range . . . . . . . . . . . . . . . . . . . . . . . 27 1.2.4 Photon interactions with matter . . . . . . . . . . . . 28 1.2.5 Photon attenuation . . . . . . . . . . . . . . . . . . . . 32 1.3 PET imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 1.3.1 Detection . . . . . . . . . . . . . . . . . . . . . . . . . . 34 1.3.2 Coincidences . . . . . . . . . . . . . . . . . . . . . . . 36 1.3.3 Sensitivity and Depth Of Interaction . . . . . . . . . . 37 1.3.4 Time of Flight PET . . . . . . . . . . . . . . . . . . . . 38

  Fig. 1.3 we can see how the distance d influences the LOR positioning. This quantity is not directly proportional to the three-dimensional (3D) range of the positron, r. The distance d is computed as the projection of the positron range distribution on the LOR [42].

2 :

 2 FIGURE 1.2: Normalized energy spectra distributions for Sc-44 (green) and F-18 (blue).

FIGURE 1. 3 :

 3 FIGURE 1.3: Scheme of positron decay: the range r is the distance from the annihilation point. The real effect is the LOR displacement, measured as the distance d between the annihilation point and the LOR.

FIGURE 1. 4 :

 4 FIGURE 1.4: The Klein-Nishina distribution of Compton scattering angles over a various range of energies.

FIGURE 1 . 5 :

 15 FIGURE 1.5: Photon interactions with matter: (a) photoelectric effect, (b) Compton scattering and (c) pair production [43].

Fig. 1 .

 1 Fig.1.[START_REF] Herrmann L Blumgart | STUDIES ON THE VELOC-ITY OF BLOOD FLOW: IV. The Velocity of Blood Flow and Its Relation to Other Aspects of the Circulation in Patients with Arteriosclerosis and in Patients with Arterial Hypertension 1[END_REF] shows the schematic functioning of a PET exam: first a radiopharmaceutical labeled with a b+ emitter is injected in the patient body. After a certain amount of time determined by the marker half life, the isotope decays emitting a positron. The particle travels in the body losing its energy until it annihilates with one of the surrounding electrons. Positron-electron annihilation results in two nearly co-linear photons, each carrying 511 keV of kinetic energy. Data acquisition starts when the patient is positioned in the scanner.

FIGURE 1. 7 :

 7 FIGURE 1.7: Complete process of a PET scan exam [53].

FIGURE 1 . 8 :

 18 photodetector used in PET imaging. The functioning is the following: a 511-keV photon is converted into a visible photon inside the scintillator, and the light particle strikes the PMT surface and, due to photoelectric absorption, some photoelectrons are emitted. The electrons are directed, through a focusing electrode, to the first of a series of dynodes. On each dynode is applied a voltage higher than on the previous dynode, in order to accelerate the electrons and multiply them through secondary emissions. The signal is then collected at the end of the PMT on the last dynode, called anode.

FIGURE 1 .

 1 FIGURE 1.11: TOF-PET scan model: A and B represent the detectors and X the position of the annihilation respect to the scanner center. The time registered in A and B is proportional to the distance TOF A and TOF B [60].

FIGURE 1 .

 1 FIGURE 1.12: (a) classical PET reconstruction, (B) TOF probability, in which the time resolution dt limits the number of voxels contributing to the event [60].

  analytical and model-based or iterative. Analytical techniques are based on back-projection, an inversion of the data, collected as projections. Statistical or iterative image reconstruction uses complex computational models, which are based on a statistical pattern of the acquisition, including all prior information about the radiotracer and the biodistribution.

FIGURE 2 . 1 :

 21 FIGURE 2.1: From projections to sinograms.[START_REF] Bailey | Physics and Instrumentation in PET[END_REF]

  2.1. EveryLOR is collected in a sinogram, a bi-dimensional array in which each voxel corresponds to a single LOR and the voxel value is incremented every time a coincidence matches the same LOR. Events along the same row in the sinogram represent events associated with parallel LORs. If the acquisition occurs in LM the events are collected separately, one by one.

FIGURE 2 . 2 :

 22 FIGURE 2.2: Projection of an image and its Radon transform [68].

transform F 2 {

 2 (x, y)} of the function f (x, y) to the 1D Fourier transform F 1 {(x, y)} of the Radon transform p(s, f) = R f (s, f), an illustration of the theorem is shown in Fig. 2.3.

FIGURE 2 . 3 :

 23 FIGURE 2.3: Central section theorem in 2D [72].

FIGURE 2 . 4 :

 24 FIGURE 2.4: Integral of a 2D plane used for the Radon transform of a 3D object [75].

FIGURE 2 . 5 :

 25 FIGURE 2.5: From left to right: reconstructions from fully 3D PET data on the left, 2D rebinned non-TOF, 3D rebinned non-TOF and fully 3D PET-TOF [82].

FIGURE 2 . 6 :

 26 FIGURE 2.6: Convergence properties of MLEM algorithm from simulated SPECT brain data. First row noise-free images, second row noisy, third row noisy filtered reconstruction [92].

FIGURE 2 . 7 :

 27 FIGURE 2.7: The images show the iterated image estimates of OSEM for SPECT brain data. For every image, we indicate the iteration number while the subsets number is reported on the left of each row [92].

3. 4 80 3. 5

 4805 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 3.4.1 Phantoms . . . . . . . . . . . . . . . . . . . . . . . . . Image reconstruction software . . . . . . . . . . . . . . . . . 82 3.5.1 Pseudo-TOF system matrix . . . . . . . . . . . . . . . 82 3.5.2 Pseudo-TOF algorithm in CASToR . . . . . . . . . . . 83 3.5.3 XEMIS2 geometry in CASToR . . . . . . . . . . . . . . 84

FIGURE 3 .

 3 FIGURE 3.1: Transversal view of the scanner to show where the small animal is positioned, how the three gammas are emitted and where they are detected.

FIGURE 3 . 2 :

 32 FIGURE 3.2: Scheme of the XEMIS2: layer and dimensions for one quarter of the scanner, the rest being similar through symmetry.

FIGURE 3 . 3 :

 33 FIGURE 3.3: XEMIS2 experimental setup at the Subatech laboratory: on the left the cryostat and on the right the camera.

FIGURE 3 . 4 :

 34 FIGURE 3.4: Decay scheme of Sc-44 [1].

FIGURE 3 . 6 :

 36 FIGURE 3.6: Scheme of the influence of the scattering angle's error on the LOR.

FIGURE 3 . 7 :

 37 FIGURE 3.7: (a) Event with no scattering. (b) Event in which one coincidence photon scatters with a small angle before reaching the LXe, causing a change in the LCI coordinates. (c) Event with an 1157-keV gamma scatter in the phantom: depending on the angle width the cone can or not intersect the LOR, anyways the LCI determination is highly affected.

FIGURE 3 . 8 :

 38 FIGURE 3.8: Scheme of the vectors involved in the spatial resolution calculation.

FIGURE 3 . 9 :

 39 FIGURE 3.9: (a) Angular resolution as a function of the scatter angle, in red the energy resolution, in green the spatial resolution, in black the total resolution. (b) The recoil energy of the electron coming from Compton scattering as a function of the scatter angle. The red lines determine the energy interval in which we have an acceptable angular resolution [4].

  3.10).

FIGURE 3 . 10 :

 310 FIGURE 3.10: Probability on the LOR for different image reconstruction techniques: (a) classical PET with a uniform probability, (b) TOF-PET with a Gaussian distribution in which the FWHM is due to the time resolution and (c) pseudo-TOF with a Gaussian distribution centered on the LCI coordinates.

FIGURE 3 . 11 :

 311 FIGURE 3.11: LCI coordinates' uncertainty scheme: d is the distance between the projection of the third-gamma emission point on the LOR and the LCI coordinates computed through the geometrical intersection.

FIGURE 3 . 13 :

 313 FIGURE 3.13: Overview of the XEMIS2 simulated as a multilayer cylinder in GATE. We can see the phantom in pink, the LXe in green and the PMTs around the active volume in red.

FIGURE 3 . 14 :

 314 FIGURE 3.14: Axial view of the central slice of the NEMA-like phantom showing the five hot spheres of 2, 4, 8, 10, and 12-mm radius.

FIGURE 3 . 15 :

 315 FIGURE 3.15: Axial view of the central slice of the Digmouse phantom showing some of the internal structures. Organ activity/voxel Organ activity/voxel [Bq/0.5 mm 3 ] [Bq/0.5 mm 3 ] Skin 0.1 Testes 0.1 Spinebone 1.7 Stomach 0.1 Eye 0.1 Spleen 3.2 Brain 0.1 Pancreas 0.1 Muscle 0.1 Liver 1.5 Bladder 2.25 Kidney 0.7 Lung 0.1 Tumor (brain) 2.25

  number of voxels), and b 2 [1, • • • , B] (B number of possible time bins on the LOR). In the TOF case, the SM elements make use of a Gaussian function to weight the LOR length component, whose FWHM depends on the time resolution.

. 14 )

 14 where i refers to the LOR bin and b to the pseudo-time bin corresponding to the considered event; l (k) j is the image value in the voxel j for the k-th estimation. We denote as pseudo-time bin the time bin corresponding to the LCI-position equivalent in ps.

FIGURE 3 .

 3 FIGURE 3.17: (a) Discretization of the XEMIS2 considering the total active volume, we obtain 234337008 virtual elements, a number too high to be treated by any software in a reasonable amount of time. (b) Discretization of the XEMIS2 considering only the first layer around the FOV. In (c) we show how the events are registered: the LOR length considered in CASToR is the length between the two red elements, whose element ID is used during reconstruction to get the event LOR position.

  image quality metrics studied in this work mainly relate to signal and contrast measurements (SNR, CNR and RC). In the following we use label with b and s subscript to indicate variables related respectively to the background and to the ROI. Due to the low statistics-induced noise in the images, a

delineated 4

 4 spheres positioned in the low-activity part of the phantom that we can consider uniform, as shown in Fig 4.1. We computed SNR and CNR as for the NEMA-like phantom, where the background variables are averaged on the 4 spherical volumes. The formulas become:

FIGURE 4 .

 4 FIGURE 4.1: Transaxial section of the Digimouse phantom mask to show the VOI (red sphere) and the 4 background (green, purple, blue and yellow).

Fig. 4 .

 4 Fig. 4.2 shows the axial view of the phantom central slice on the last iteration for the three reconstruction methods: (a) 2-g, (b) 70-ps pseudo-TOF and (c) pseudo-TOF variant resolution reconstruction, and the relative plot profiles along the yellow line in (d) and (e) vs (f).

Fig. 4 .FIGURE 4 . 2 :

 442 Fig. 4.3 shows the same view of Fig. 4.2 after applying a Gaussian filter of scale s = 1 to reduce noise induced by low statistics. The metrics in Fig. 4.4 were calculated following the equations in section 4.2, normalizing the values inside the cylinder. We show the SNR, CNR

FIGURE 4 . 3 :

 43 FIGURE 4.3: Axial view of the NEMA phantom after Gaussian filtering (s = 1 mm) using (a) classical PET reconstruction and the proposed 3-g technique with (b) 70-ps and (c) variant standard deviation. Images are the results of the 30 th iteration. In (d), (e) and (f) the corresponding plot profiles along the yellow line. The high contrast ring around the FOV in the first image on the left is due to a CASToR misinterpretation of the activity value on the FOV edge and it does not affect image reconstruction.

FIGURE 4 . 4 :

 44 FIGURE 4.4: Sphere-by-sphere metrics results, from top to bottom the (1) 2-mm, (2) 4-mm, (3) 8-mm, (4) 10-mm, and (5) 12-mm radius sphere with the respective values of SNR, CNR and RC.

FIGURE 4 . 5 :FIGURE 4 . 6 :

 4546 FIGURE 4.5: Reconstructed images of the Digimouse, from left to right: 2-g reconstruction, pseudo-TOF with 200-ps, 150-ps, 100-ps, and 70-ps resolution.

FIGURE 4 . 7 :

 47 FIGURE 4.7: Line profile plot of the mouse's head's hot spot for all the different reconstructions on the last iteration.

  (e) and (f) vs (d), where the lower values indicate less reconstructed activity in the hot spots. Looking at the reconstructed images after filtering (Fig.4.[START_REF] Holland | 89Zr-DFO-J591 for immunoPET of prostatespecific membrane antigen expression in vivo[END_REF] we can appreciate even more the contrast enhancement in the 3-g case (b) and (c), quantitatively visible in the profiles (e) and (f) versus (d) 2-g.One of the main amelioration brought by the XEMIS2 geometry can be noticed observing the axial view of the NEMA phantom. All five spheres appear perfectly round and in the exact position they are expected to be, with no distortion. Despite the depth of the detector (12 cm), no DOI correction was applied, since the XEMIS2 technology is able to detect the three-coordinate position of interaction, which is the main advantage of our monolithic detector.Another quantitative analysis of the amelioration brought by the techniqueproposed in this work, can be found in the metrics plots, Fig.4.4. The improvement is visible in the SNR, CNR and RC higher values, and it is mostly due to the technique's ability to narrow the decay's position on the LOR and to focus on the right ROIs[START_REF] Laval | Design and performance of a new positron computed tomograph (PCT) using the time-of-flight information[END_REF][START_REF] Joel | Benefit of time-of-flight in PET: experimental and clinical results[END_REF]. While we know that for conventional PET reconstruction the decay's probability is uniform on the LOR. The higher CNR and SNR, together with the fast convergence, are the main advantages of the proposed pseudo-TOF algorithm. The acceleration in the convergence rate can be seen observing the SNR, CNR and RC peak in the plots in Fig.4.4. For all five spheres, we reach the maximum value earlier in the 3-g reconstruction case than in conventional image reconstruction.In Tables 4.1 and 4.2 we present some values related to the SNR and CNR percentage gain of 3-g over 2-g, both for variant and for 70-ps pseudo-TOFresolution. The values related to the smallest sphere are surprisingly high, while in a TOF-like reconstruction we do not expect to gain in SNR or CNR for very little objects. It is most likely due to the fact that in the 2-g case the sphere is confounded in the background, and again we do not consider these result well grounded. For the other spheres, we can notice how the proposed method shows an evident percentage gain, which spans between the 28% and 34% and that is dependent on the object's size. The best values refer to the 70-ps resolution reconstruction, supporting the choice of the RMS value in the distance distribution as pseudo-TOF resolution.The results of the Digimouse phantom differ from those presented using the previous studied objects. The Digimouse phantom presents a complex internal geometry, made of small and detailed substructures, which are challenging to reconstruct in low count conditions. While the cylindrical NEMA phantom is simulated as a superposition of simple GATE objects, the small animal is defined voxel by voxel in GATE. Due to the Digimouse body and organs' dimensions, the voxel size was reduced to 0.8 mm 3 to avoid the reconstruction of artifacts in the images. Since the calculation of the pseudo-TOF resolution requires a geometrical intersection between the Compton cone and the LOR, the error on the LCI points is proportional to the distance between the cone center and the line. Any uncertainty on the cone determination, for example concerning the Compton angle miscalculation, is directly reported on the LCI coordinates as an error increasing with the distance between the LOR and the cone vertex. The Digimouse is a more complex phantom than the NEMA, with 16 different materials and small organs, which increase the uncertainty in the LCI, due to Compton scattering in the phantom and to the positron range, which is increased in low-density materials. Furthermore, the lack of information about the interactions happening inside the Digimouse phantom, due to GATE simulation output for voxelized phantoms, led us to use pseudo-TOF resolution values that were not adapted to the simulated object, since the study on the distance between the third-gamma emission point and the LCI position on the LOR could not be carried out. Fig. 4.5 shows results for all the reconstructions tested: the pseudo-TOF resolution values chosen for the Digimouse reconstructions were taken from the NEMA phantom study going from the average to the maximum distance. They spanned from 70 ps (e) to 200 ps (b), which corresponds to a range from 22 to 60 mm of uncertainty on the LOR, typical state of the art TOF values for preclinical scanners. In order to gain image quality on small animal's dimensions, we ought to work with very precise resolution. The 22-mm resolution on the LOR calculated for the NEMA phantom represents a sufficiently high value to show improvement over conventional PET reconstruction for a larger object. Lowering the resolution leads to an approximation of the 2-g case, which can be seen in the Digimouse 3-g image reconstruction with 200-ps of pseudo-TOF resolution, comparing Fig. 4.5 (a) and (b).

5. 5 109 5. 5 . 1 109 5. 5 . 2 110 5. 6

 510951109521106 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figures of Merit . . . . . . . . . . . . . . . . . . . . . . Invariant and variant kernel convolution . . . . . . . Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

FIGURE 5 . 1 :

 51 FIGURE 5.1: Fractional energy loss for positron and electrons in lead as a function of the particle energy [189].

FIGURE 5 . 2 :

 52 FIGURE 5.2: Reconstructed images using different positron range correction, b indicates the noise level [203].

42 TABLE 5 . 2 :

 4252 Sphere nb Position [mm] Radius [mm] Material Density [g/cm 3 Description of the spheres in the cylindrical phantom: reference number, position in mm considering the center of the FOV as the origin, and radius size.

FIGURE 5 . 3 :

 53 FIGURE 5.3: Axial view of the central slice of the phantom showing the three spheres of 2, 4, and 8-mm radius made of water and the two with a 10 and 12-mm radius made of lung and bone.

FIGURE 5 . 4 :

 54 FIGURE 5.4: Transaxial view of the central slice of the PSF study phantoms showing the three point-source spheres made of bone, water and lung, from left to right.

FIGURE 5 . 5 :

 55 FIGURE 5.5: PSF plot and Gaussian fit for the bone sphere in (0, 0, -60).

  metrics employed in the positron range image analysis are the SNR and the CNR, that were calculated as follows: the label b refers to the background and the label s to the hot spot.

FIGURE 5 . 8 :

 58 FIGURE 5.8: Transaxial view of the NEMA-like phantom, reconstructed with four different techniques: (a) conventional 2-g PET and (b) 3-g, with no positron range correction; (c) 2-g and (d) 3-g reconstruction convolved with a variant kernel, (e) 2-g and (f) 3-g reconstruction convolved with a stationary kernel.

FIGURE 5 . 9 :

 59 FIGURE 5.9: Plot profiles of the six different reconstructions through the yellow line seen in Fig. 5.8.

FIGURE 5 . 10 :

 510 FIGURE 5.10: SNR and CNR plots of the three smallest spheres, water-made: in numerical order sphere of 2-mm, 4-mm and 8-mm radius. With "conv" we refer to variant kernel convolution, while "stat conv" refers to the invariant kernel convolution correction.

FIGURE 5 . 11 :

 511 FIGURE 5.11: SNR and CNR plots of the two biggest spheres: sphere number 4, 10-mm radius with lung density and sphere number 5, 12-mm radius with bone density. With "conv" we refer to variant kernel convolution, while "stat conv" refers to the invariant kernel convolution correction.

(

  Subatech, 2017) consists in accepting only events with at least 4 interactions, one for each coincidence photon and two for the third gamma to build the Compton cone. The division in three groups is based on a distance-study to determine a distance threshold. Through a long realistic GATE simulation he studied the distribution of distances between interactions belonging to the same photon or to different photons, and established that these cases follow two distributions partly overlapping (Fig.5.12). The distance of D cut = 74 mm is taken as the distance-cut to determine if the interaction belongs or not to the same group. The distinction between the third-gamma group and the 511-keV-photons groups is based on the group released energy amount: since the prompt gamma is more energetic (1.157 keV) the group with more than 511-keV energy is the one used to compute the Compton cone.

FIGURE 5 . 12 :

 512 FIGURE 5.12: Distance study to determine the D cut to use in the interactions clustering.

8 et Table 1 pour plus de détails. L'activité totale simulée dans le fantôme est de 20 kBq, avec un facteur 15 de rapport de contraste entre la sphère et l'arrière-plan. 30 doublons ont été

	third gamma
	emission point
	LCI coordinates
	d
	FIGURE 6: Étude sur d, la distance entre la projection du point
	d'émission du troisième gamma sur la LOR et les cordonnées
	LCI.
	simulés pour étudier la variance de l'image.
	Pour compléter la caractérisation, en tant que deuxième fantôme, nous
	avons choisi de simuler un objet plus réaliste: une souris mâle de 28 g (fantôme
	Digimouse). Le petit animal présente une structure complète, composée du
	cerveau, des muscles, des yeux, des glandes, du coeur, des poumons, du foie,
	de l'estomac, de la rate, du pancréas, des reins, des testicules, de la vessie, du
	squelette et de la peau. Une tumeur sphérique de rayon de 2 mm a été ajoutée
	Le fantôme est constitué d'un cylindre d'eau homogène de 12 cm de long dans le cerveau en tant que ROI à analyser pour l'étude de qualité d'image.
	et d'un rayon de 3,5 cm, contenant cinq sphères. Les centres de toutes les Le figure 9 montre une vue transaxiale de la tranche centrale de Digimouse et
	sphères sont équidistants de l'axe central du cylindre et positionnés sur la le tableau 2 la distribution simulée d'absorption de la souris.
	tranche centrale du cylindre. Le rayon des sphères mesure 2, 4, 8, 10 et 12 mm, Le radiopharmaceutique utilisé dans cette simulation est le [ 44 Sc]Sc-J591

[mm]

FIGURE 7: Incertitude sur la position du LCI.

voir Fig.

TABLE 1

 1

	: Dimensions et positions des cinq sphères dans le fan-
	tôme.
	marqué au Sc-44, un anticorps recombinant monoclonal marqueur de l'antigène
	PSMA spécifique de la membrane prostatique l'absorption biologique est
	connue de Carter [2] et des études de Holland [3]. Un des avantages de
	l'utilisation du [ 44 Sc]Sc-J591 est que la molécule est également utilisée pour
	les études de parcours du positon puisqu'elle est facilement associée aussi au
	F-18.
	L'activité totale simulée dans le fantôme était d'environ 22 kBq et le temps
	d'acquisition des données était de 20 minutes. En raison de la structure plus
	détaillée du fantôme Digimouse, un temps de calcul plus long a été nécessaire
	et une seule acquisition a été simulée.
	Toutes les images que nous avons reconstruites dans ce travail sont obtenues
	grâce à l'algorithme de reconstruction LM-MLEM dans CASTOR v2.3, pour le
	pseudo-TOF et le non-TOF, que nous appelons respectivement 3-g et 2-g .
	Comme pour tous les algorithmes d'optimisation itérative comme MLEM,
	la reconstruction dans CASTOR est basée sur un projecteur gérant à la fois la

TABLE 2

 2

	: Biodistribution du [ 44 Sc]Sc-J591 dans le corp du Digi-
	mouse.
	projection arrière et avant et il y a un optimiseur et un convolveur, au cas où
	tout type de convolution serait nécessaire (par exemple pour la correction du
	parcours du positon).
	Dans CASTOR, le calcul des éléments de la matrice système suit le pro-
	jecteur Siddon et la méthode dite ray-tracer, pour laquelle la probabilité dépend
	de la partie de LOR traversant chaque voxel. L'optimiseur standard utilisé
	pour les événements à deux gamma est de Shepp et Vardi; tandis que pour
	l'algorithme LM-MLEM/TOF nous nous référons à Filipovic.

TABLE 3 :

 3 Description des sphères a l'intérieur du fantôme.

TABLE 1

 1

.1: Mean positron range in water for isotopes of interest in PET

[START_REF] Bai | Positron range modeling for statistical PET image reconstruction[END_REF][START_REF] Michael | Effect of positron range on spatial resolution[END_REF]

.

TABLE 1 .

 1

	or plastic as well as liquid or gaseous (e.g. LXe, argon) [54]. Table 1.3 shows
	some conventional scintillators and their density, atomic number, and decay
	constant.

3: Physical properties of some common detectors in PET.

1.3.1 Detection

Conventional PET cameras are composed of a first layer of detector mean (e.g. scintillators) and a second layer of photo-detectors (e.g. photomultiplier tubes). The detection material needs to convert the radiation signal into visible light and to fulfill some properties in order to be suitable for PET imaging: good stopping power for 511-keV photons, short signal decay time, high light output, and good energy resolution. Some of the most common PET scintillators are organic, such as bismuth germanate (Bi 4 Ge 3 O 12 ), lutetium oxyorthosilicate (Lu 2 SiO 5 ), gadolinium oxyorthosilicate (Gd 2 SiO 5 ),

  human head) we have a gain around 2.1, with a scanner time resolution of 300 ps. No high improvements in terms of SNR is expected in TOF small animal systems, however for clinical scanners the SNR gain is exploited to reduce the acquisition time and the injected dose [63]. Detector geometry . . . . . . . . . . . . . . . . . . . . 59 2.6.2 Attenuation and normalization . . . . . . . . . . . . . 59 2.6.3 Positron range . . . . . . . . . . . . . . . . . . . . . . . 60

2.5 System modelization and corrections . . . . . . . . . . . . . 58 2.6 System matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 2.6.1 2.7 Scattering and random coincidences . . . . . . . . . . . . . 60

TABLE 3 .

 3 2: Radiopharmaceutical uptake distribution for the mouse phantom with 0.8 ⇥ 0.8 ⇥ 0.8-mm 3 voxels: [ 44 Sc]Sc-J591 uptake values in the Digimouse computed from the data in Holland, 2010[START_REF] Holland | 89Zr-DFO-J591 for immunoPET of prostatespecific membrane antigen expression in vivo[END_REF].

TABLE 4 .

 4 1: Percentage gain in SNR for pseudo-TOF over non-TOF for all spheres, computed on the last iteration.

	Radius CNR % gain	CNR % gain
	mm	70 ps	variant resolution
	2	43.1	46.0
	4	34.1	33.3
	8	30.1	28.4
	10	29.3	28.0
	12	28.8	27.6

TABLE 4 .

 4 2: Percentage gain in CNR for pseudo-TOF over non-TOF for all spheres, computed on the last iteration.

Table 5

 5

	.1 some common radionuclides with the corresponding estimation of
	positron range.
	Isotope Half-life E mean E max	R mean in	R max in	Ref.
		[h]	[keV] [keV] water [mm] water [mm]
	18 F	<2	250	635	0.5	2.3	[165]
	64 Cu	12.7	216	653	0.6	2.9	[166]
	89 Zr	78.4	403	902	1.3	4.2	[167]
	45 Ti	3.08	442	1040	1.5	5.2	[168]
	44 Sc	3.97	632	1474	2.4	6.9	[37, 169]
	86 Y	14.7	640	3141	2.5	11.1	[170]

TABLE 5 .

 5

1: Some positron range properties of common radioisotopes used in PET imaging.

TABLE 5 .

 5 3: Transaxial and axial FWHM values of the kernel used to correct for Sc-44 positron range in water, bone and lung.

-a qualitative and quantitative analysis comparing conventional and pseudo-TOF reconstruction; -realization of the first reconstructed images for XEMIS2, through conventional and 3-g reconstruction.

Chapter 2

Image Reconstruction Summary

This chapter is an introduction to image reconstruction as historically developed, from analytical to iterative algorithms. We begin introducing the concept of sinogram and projection in data acquisition and we present the image reconstruction theory starting with analytical methods, in two-dimensional (2D) and in 3D. Afterward we pass over to iterative image reconstruction algorithms mentioning many of them (algebraic reconstruction technique (ART), simultaneous ART (SART), simultaneous iterative reconstruction technique (SIRT), multiplicative ART (MART)) and going in the specific of MLEM and OSEM. At last, we introduce the system modelization through the system matrix (SM), with a look into the correction factors for geometrical effects, attenuation and positron range.

3-g Image Reconstruction

Summary

In this chapter, we get into the details regarding the geometry and the functioning of the XEMIS2 and we present the proposed 3-g image reconstruction algorithm. We expose the calculation of the LCI coordinates, using Compton kinematics and LOR/Compton cone intersection. We then present the 3-g or pseudo-TOF algorithm, whose name is due to the similarity between the proposed technique and TOF image reconstruction.

We also present the Monte Carlo GATE simulation used for data acquisition, and the simulated phantoms. At the end of the chapter, we deal with the details of the image reconstruction software, CASToR. We explain the theory behind the 3-g image reconstruction algorithm, defining the SM elements and presenting the implementation of the reconstruction algorithm into the CASToR framework.

XEMIS2 geometry in CASToR

The first information required by CASToR to start the image reconstruction process is the detector geometry. Normally, through the use of a tool provided with the package, we are able to directly transfer the GATE geometry into a CASToR geometry file. The software needs the number of modules, sectors and crystal elements to rebuild the coincidences and identify them with its own notation (CASToR-ID), in order to recognize the events in its system of reference.

This functioning cannot be applied in a monolithic detector case, due to the absence of detector elements. In this case, we need to produce a geometry file in which the detector is discretized in smaller elements, with dimensions comparable to the spatial resolution of the system to not lose information (Fig. 3.16).

In our case, we chose to discretize with 1-mm 3 virtual elements, which brought us to a total of more than 23 millions elements due to the large size of the XEMIS2 active volume, that was not acceptable in terms of computational time to process. The solution to this matter was the discretization of only a thin layer of the active zone around the FOV, as shown in Fig. 3.17 (b). We considered a 0.5-mm-deep layer of elements all along the FOV cylinder for a total amount of Chapter 4

Results

Summary

In this chapter we present the results obtained in this work regarding image reconstruction with XEMIS2. First, we introduce the figures of merit employed in the image analysis and we then show the results for the two simulated phantoms. We compare 3-g images to images obtained through conventional PET reconstruction for both phantoms in XEMIS2. We show the reconstructed images, ROIs's plot profiles, CNR, SNR and recovery coefficient (RC) trends, and we briefly discuss the obtained results.

Introduction

In this chapter we discuss the results obtained for image reconstruction with the XEMIS2. We hold a comparison between 2-g and 3-g reconstructed images in order to assess the benefits of the pseudo-TOF technique over conventional PET reconstruction.

Positron range estimation

The positron range estimation is the basis for the effect's correction. Two different approaches can be followed: an analytical modelization in which the range and other parameters are obtained through a Monte Carlo simulation, and a more experimental way based on actual measurements.

In particular positron range estimation's experiments were conducted starting in 1975 by Phelps [START_REF] Michael | Effect of positron range on spatial resolution[END_REF], Cho [START_REF] Cho | Positron ranges obtained from biomedically important positron-emitting radionuclides[END_REF], Derenzo (1986) [START_REF] Stephen E Derenzo | Mathematical removal of positron range blurring in high resolution tomography[END_REF] and Levin (1999) [START_REF] Craig | Calculation of positron range and its effect on the fundamental limit of positron emission tomography system spatial resolution[END_REF], while the interest in PET imaging was growing.

On this occasion, the words range and path length were defined: the range as the penetration depth and the path as the integral of the reciprocal of the stopping power, an approximation of the range.

In fact, since the positron undergoes several elastic and inelastic interactions while propagating in the medium, its path appears larger than its range, and due to this difference a correct estimation of the range was needed.

Early experiments were not completely successful since the resolution of the detector employed in the measurements was comparable to the positron range [START_REF] Cal-González | Positron range estimations with PeneloPET[END_REF]. Derenzo experiments in 1979 [START_REF] Stephen E Derenzo | Precision measurement of annihilation point spread distributions for medically important positron emitters[END_REF] were centered on the measurement of the positron range in polyurethane foam and in aluminum. The idea is to remove the blurring caused by this effect through a deconvolution.

The main problem with this approach is the loss of information in the deconvolution and the need for scaling the range to a water-equivalent material other than foam, since the range does not exclusively depend on the material density.

Later Palmer and Brownell [START_REF] Matthew | Annihilation density distribution calculations for medically important positron emitters[END_REF] [START_REF] Lehnert | Analytical positron range modelling in heterogeneous media for PET Monte Carlo simulation[END_REF]. The simulation of the positron range can be carried out using different frameworks and software: GATE [START_REF] Sébastien | GATE: a simulation toolkit for PET and SPECT[END_REF] with its option empenelope allows accurate positron interactions, as well as PeneloPET [START_REF] España | PeneloPET, a Monte Carlo PET simulation tool based on PENELOPE: features and validation[END_REF] based on PENELOPE [START_REF] Salvat | PENELOPE-2006: A code system for Monte Carlo simulation of electron and photon transport[END_REF] and Geant4. Results are plotted in function of the iteration number, considering each sphere as a separate VOI and the rest of the cylinder as the background.

Invariant and variant kernel convolution

We reconstructed images via conventional PET reconstruction and via pseudo-TOF reconstruction, using a 70-ps resolution. In both algorithms, we added the positron range correction through variant and invariant kernel convolution.

For the invariant kernel case we convolved with a FWHM resulting from the average of the values in Table 5.3, i.e. 2.0 mm. For the variant case, the 10-mm-radius sphere and the 12-mm-radius sphere were associated with the corresponding values in Table 5.3 and the rest of the image was convolved with a FWHM equal to 1.9 mm, considering that we are in presence of water and air.

In In Fig. 5.9 we show the yellow line profiles: the line passes through the two spheres with the bone (12-mm radius) and lung (10-mm radius) density.

The plots concerning the SNR and CNR values for the three water-dense spheres are showed in Fig. 5.10 while in Fig. 5.11 we show the same results

for the other two spheres.

Conclusion and perspectives

The work presented in this thesis focused on the development of a novel reconstruction technique for 3-g imaging, which was evaluated on the XEMIS2 imaging system geometry. The XEMIS2 study is the preclinical part of the more general XEMIS project, which already developed a first prototype, the XEMIS1. The characterization of the XEMIS1 Compton camera gave the experimental demonstration of the feasibility of 3-g imaging, we find the previous studies in LXe , Oger [START_REF] Oger | Développement expérimental d'un télescope Compton au xenon liquide pour l'imagerie médicale fonctionnelle[END_REF], and Grignon [START_REF] Grignon | Étude et développement d'un télescope Compton au xénon liquide dédié à l'imagerie médicale fonctionnelle[END_REF].

In this work we use the third-gamma interactions in LXe to determine a Compton cone, whose point of intersection with the LOR, obtained from the coincidence photons, allows to narrow the localization of the LOR section involved in the combined annihilation + emission event. The event is then reconstructed using a pseudo-TOF MLEM algorithm where the line-cone intersection is taken as the center of a Gaussian PDF used in the reconstruction, in a similar fashion as in TOF-PET, whence the denomination of the proposed method as pseudo-TOF. A main advantage associated with the proposed approach is to allow for seamless embedding of the third gamma information through existing reconstruction frameworks with TOF capability, such as CASToR. Doing so, we leverage state of the art reconstruction techniques, which is an important step towards feasibility in a real clinical setting.

We showed a simulation-based study of a realistic 20-minute data acquisition with 20 kBq of initial radioactivity in the XEMIS2 detector. We simulated a NEMA-like and Digimouse phantom and showed image quality improvement for objects of various densities and dimensions. In Chapter 3 we presented the LCI coordinate calculation, on which we base the 3-g image reconstruction method. In particular in section 3. Abstract: In this thesis we present three-gamma imaging, where the acquisition system relies on a beta+ and gamma emitter. The rationale of 3-gamma imaging is that the third gamma detection information may help to provide better localization of the annihilation point, thus enabling higher image quality and fewer dose delivered to the patient. We present the 3-gamma system, the XEMIS2, developed at Subatech, Nantes, that is a LiquidXenon detector suitable for 3-gamma imaging due to its stopping power, its scintillation characteristics and its continuous geometry. The principle of 3-gamma image reconstruction is based on the intersection of a LOR, obtained from the coincidence photons, with a Compton cone, determined by the third gamma. The idea is to find the LOR\cone intersection and use it to locate the most probable annihilation position on the line, as for the time difference in TOF-PET. We present a complete GATE simulation study of two phantoms (similar-NEMA and Digimouse), to assess the improvements of 3-gamma image reconstruction over conventional PET and we study the positron range correction, which is important for our beta+gamma emitter, Sc44.