Je Remercie

Uwe Bergman

Guillaume Bal

Mohammed El Rhabi

Loic Bertrand

Mathieu Thoury

Sebastien Schoeder

Regina Oprandi

Sophie David

Pierre Gueriau

Giulia, Céline, pardon Cécile, Josiane, Susy, Rafaella, Tulin, Selwin Dominique Une Mention Aussi Pour Jiayi

Jeremy Luis

Nadège, Liliane, Tahar, Catherine, Vahagn, Christophe, Mohamed Christine Laure

PE Sybille Salim

Hélène Thomas

Maxime Jonas Hugo

Colin Felix

Arsen, Tamara, Florien, Sébastian, Mamadou... Benjamin

François

Géraldine

Antoine Florien

Dr Guil- Laume Omar

Martina Fernando

Oscar

Esteban Diego

Mabe, Martín, Eli,... y sus Daniela Yandira

Amarunas

Mis Otros Colegas Latinoamericanos Max

Pavo

Luis

Yadira Lucas

J'ai

Diego

Adrian Santiago

Edmundo Nicky Christian

Carlos Nestor

Mi Papá

Mi Mamá

M K Guerrero

L Nguyen

S X Dumas

Cohen

X-Ray

la excellente ambiance de travail qu'on y trouve

qui m'ont apporté.

Introduction

The present work deals with the mathematical aspects of a radiation imaging modality to characterize heterogeneous ancient materials objects. Such characterization is performed in terms of the electronic density of objects. Indeed, the knowledge of the internal electronic density distribution of an object provides us informations about the morphology and chemical composition of the object and then a study focused on other research fields, such as material sciences or art history, can be realized [2]. Unfortunately, the flattened geometry of the objects prohibits some well studied imaging techniques like photo-absorption or phase contrast tomography. The radiation imaging modality used here, is based on Compton Scattering Tomography, developed in recent years mostly for applications in emission imaging techniques [START_REF] Morvidone | On the V-line radon transform and its imaging applications[END_REF][START_REF] Nguyen | Radon transforms on a class of cones with fixed axis direction[END_REF][START_REF] Truong | On new V-line radon transforms in R 2 and their inversion[END_REF][START_REF] Rigaud | Modeling and simulation results on a new compton scattering tomography modality[END_REF]. More precisely, it consists in a external scanning imaging system allowing a three-dimensional mapping of heritage objects.

Modelling of imaging systems requires to study both the direct and inverse problems. The direct problem is called the image formation process, which aim is to provide a model of the image generated from a given input object. The inverse problem is called object reconstruction, or inversion, and is focused in providing an approximation of the distribution of an unknown specific physical characteristic of the object. In Compton Scattering Tomography, the output image is an spectral image resolved both in energy and in space whereas the reconstructed object is a mapping of the electronic density distribution of the input object.

After a theoretical study of the direct and inverse problem, numerical simulations have been performed on discrete input objects, representative of flattened ancient material objects. The output spectral image is obtained through a numerical implementation of the direct problem. A numerical reconstruction of the object from the simulated image is then performed by solving the inverse problem. The quality of the reconstructed object is finally measured by calculating some distance between the initial object and its reconstruction.

Two major objectives can be identified in this PhD thesis. The first one is to develop a mathematical model of the direct problem of image formation. This model must take into account all the parameters encountered in a realistic scenario and then to be able to numerically estimate for a given input object the output spectral image. The image formation process will be compared with Monte Carlo simulations of photon transport as defined in [START_REF] Salvat | PENELOPE, a code system for Monte Carlo simulation of electron and photon transport[END_REF]. The second main objective of this thesis is the definition of a reconstruction procedure for the inverse problem. An analytical inverse formula will be obtained and then numerically discretized to lead to a numerical reconstruction of the input object. A real synchrotron experience collecting Compton scattering data under the configuration proposed is the remaining step to fully validate the presented approach.

A more detailed description of the particular context of ancient material objects and the associated experimental setup is proposed below. This introduction ends with a brief outline of the next chapters.

Specificities of ancient materiel 3D studies

Ancient material objects are remnants of the past that have historical value. Such objects are encountered in archaeology, palaeoenvironments, palaeontology, or cultural heritage. The study of their chemical composition and their morphology not only provides informations about their historical context (origin, use...) but also allows to understand their alteration processes and suggest the most suitable strategies for their restoration and conservation. Details about the specific characteristics of heterogeneous ancient materials in synchrotron examinations can be found in [2].

Imaging or spectroscopic studies, such as absorption, emission or phase contrast tomography among other X-ray based imaging techniques, have been developed to access to this information. When one deals with heritage objects, the non-invasiveness and nondestructiveness properties of inspections are a requirement that X-ray imaging methodology provides, enabling two-dimensional imaging of the sample. Nevertheless, one can easily be faced with samples presenting a flattened geometry, that is, samples presenting a large ratio between their front area and thickness. The challenge is then to perform a three-dimensional probing without using a relative rotation between the sample and the imaging system as it would be done in conventional tomography, using either absorption or phase contrast modality, since probing would suffer from the high differential light path in distinct directions.

Samples presenting such characteristics are encountered in particular in studies of conservation and restoration of easel painting. Indeed, they require the characterization of the stratigraphical assemblage of pigments often over a very dense background layer such as one made of lead white. An example of this kind of studies is shown in Figure 1.1 which was performed by means of an invasive method. Figure 1.2 presents another example of objects possessing this morphology, namely in palaeontology with the Lagerstätten fossils [START_REF] Gueriau | Trace elemental imaging of rare earth elements discriminates tissues at microscale in flat fossils[END_REF] which are mechanically flattened during the fossilisation process and stand on one side of a thick sedimental slab which can not be thinned for the study. In both cases, the volume of interest forms a layer on top of a material support which is opaque to X-ray either due to its density or its thickness.

Two alternatives are available to work out this issue : either to perform a stratigraphical section of the sample which is an invasive method, or to limit the study to a bi-dimensional analysis of the front surface of the sample, for example with synchrotron X-ray fluorescence spectral raster-scanning as performed in [START_REF] Gueriau | Trace elemental imaging of rare earth elements discriminates tissues at microscale in flat fossils[END_REF], remaining however a three-dimensional study still to be done.

To overcome the limitations of the mentioned methodologies for such samples, an imaging modality based on Compton Scattering Tomography is proposed in this work.

Main principles of Compton Scattering Tomography

Classical X-ray tomography considers mainly four basic photon-matter interaction phenomena, namely, the photoelectric absorption, electron-positron pair production, Rayleigh scattering which is both elastic and coherent, and Compton scattering which conversely is both inelastic and incoherent. In X-ray imaging and tomography, Compton scattered signal is considered as noise added to photoelectric absorption and coherent scattered data in the image formation process. This is because X-ray transmission signal is dominated by the photoelectric absorption whilst coherent scattering may produce significant amplitude variations at low scattering angles thanks to constructive and destructive interference effects due to the coherent nature of this scattering.

In Compton scattering tomography (CST), an incident photon of energy E 0 is absorbed by a target electron, who re-emits a secondary photon scattered by an angle ω relative to the direction of the original photon (see Figure 1.3). The scattered photon has an energy E ω which is related to the scattering angle ω by the identity (1.1) called the Compton equation

E ω = E 0 1 + E 0 m e c 2 (1 -cos ω) , (1.1)
where m e c 2 = 511 keV is the rest mass energy of the electron. Depending on the material, if the incident radiation has a larger energy of about 4×10 4 electron-volts, Compton scattering becomes the dominant phenomenon in the process, even more when detection is performed outside the direct transmission area (ω = 0). Under these conditions, by using CST, we are able to avoid relative rotations between the sample and the imaging system. Furthermore, because of a potentially opaque supporting material, transmission and forward Compton scattering data, that is, with a scattering angle inferior of π 2 , is impossible to collect. This motivates the proposal of a modality tapping on back-scattered data, that is, data collected with a scattering angle ω comprised between π 2 and π (see figure 1.5). The idea of exploiting scattered radiation by Compton effect in imaging techniques has been introduced and studied simultaneously. It has given birth to Compton scattering tomography [6,[START_REF] Morvidone | On the V-line radon transform and its imaging applications[END_REF][START_REF] Nguyen | Radon transforms on a class of cones with fixed axis direction[END_REF][START_REF] Norton | Compton scattering tomography[END_REF][START_REF] Truong | On new V-line radon transforms in R 2 and their inversion[END_REF], which focuses on reconstructing the electron density map of the object. The first CST scanner was proposed in 1994 by Norton [START_REF] Norton | Compton scattering tomography[END_REF] through a Radon transform over circle arcs starting at a γ-ray source and ending at a moving detecting site. Radon transforms over conical surfaces having fixed axis directions and variable opening angle were proposed and studied in [START_REF] Nguyen | Radon transforms on a class of cones with fixed axis direction[END_REF] where a first analytic inversion formula is proposed using circular component analysis, with applications to emission imaging based on Compton scattered radiation. Generalisations of this kind of transforms to higher dimensions spaces with related inverse formulae in a filtered back-projection type are presented in [START_REF] Haltmeier | Exact reconstruction formulas for a radon transform over cones[END_REF]. A back-projection inversion algorithm for a conical Radon transform in R 3 was developed recently in [6].

The imaging configuration device

Compton Scattering Tomography aims to reconstruct the electron density map of the studied object. In this work, the electron density is represented mathematically by a nonnegative function defined on the Schwartz space (both smooth and rapidly decreasing) on R 2 f : R 2 × R + → R + .

A synchrotron radiation setup with a parallel monochromatic X-ray beam (about 50 keV) and a space-energy resolved detector are considered. As presented in Fig. 1.4, an incident photon of energy E 0 , in the Oz-axis direction and with a square section, is Compton scattered by an electron situated inside the object at M subtending an angle ω varying in π 2 < ω < π with the direction of incidence. The scattered photon, of energy E ω approximated by (1.1), reaches a detecting site D = (ζ, ξ, 0) on the 2D detector that is located over the xOy-plane.

As represented in Fig. 1.4, the detector will be placed between the source and the object to capture back-scattered photons. It will have a hole in the middle of area 4ζ 0 ξ 0 for some two positive real numbers (ζ 0 , ξ 0 ) to allow the beam to go through. Therefore, we will have data for values of (ζ, ξ) on the xOy-plane verifying |ζ| > ζ 0 or |ξ| > ξ 0 . Horizontal and vertical translations of the sample will be needed to allow the imaging of the full object. When all electrons of the material are both free and at rest, the Compton equation establishes in (1.1) a diffeomorphism between the scattering angle and the scattered energy, provided that the incident beam is monochromatic.

In such circumstances the Compton scattering image formation process mathematically corresponds to a Radon transform over the surface of a cone, called the conical Radon transform.

For sake of clarity, a 2D configuration is going to be studied first. This is due to fact that the azimuthal scattering angle is uniformly distributed in a Compton event, that is to say, the angular distribution of scattered photons has axial symmetry around the direction of incidence. This framework is represented in Figure 1.5. Here, the detector is located over the Oy-axis of the Cartesian reference plane and noted D ζ = (0, ζ). This 1D detector is placed between the source and the object to capture back-scattered photons as in the 3D case and it will have a hole in the middle of length 2ζ 0 to allow the beam to go through. We will thus have data for values of ζ verifying |ζ| > ζ 0 for some ζ 0 > 0. Vertical translations of the sample will also be needed. In this case, the integration is performed over two lines in a half-space, leading to the so-called V-line Radon transform.

Brief outline of chapters 2 to 6

The Radon transform models the image formation process in classical photo-absorption tomography. In the simplest two-dimensional case, a continuous function f belonging to the object input space X of functions defined on R 2 is given. This integral transform returns a function Rf defined in the image output space Y that models the output of the imaging system. This output function performs integrations over lines in R 2 of the input function following a given radiation direction y that characterizes the line. Each line integral is called a projection of the object, and can be written as

Rf (y) = y f (x) dx.
(1.2)

Solving this last equation for f is the main issue of the inverse problem, together to a workable numerical implementation for large output datasets.

Chapter 2 reviews the direct and inverse problem for the originally introduced Radon transform, covering the continuous and discrete framework. The Radon transform in cartesian coordinates is studied as a preliminary step to the half space-Radon transform, which turns to be itself a preliminary step to the V-line Radon transform exposed in chapter 3. For the cartesian transform, functions in the output image space are defined on the cartesian space Y = R 2 . The half-space Radon transform has the same output image space Y but the input object space is restricted to functions defined on R + × R. Continuous and discrete formulations of the direct and inverse problem are treated and illustrated by numerical simulations for each Radon transform.

Various generalizations of the Radon transform has been studied since the seminal Radon work. They are mostly applied in imaging sciences where integration is performed over different manifolds according to the imaging configuration. The proposed Compton scattering modality studied in this thesis is first restricted to a two-dimensional framework for a first understanding of the model in Chapter 3. In that case, integrations of the input function must be performed over two half lines with common vertex that forms a V-line in R 2 . The V-line Radon transform V is applied on functions defined on R + × R whereas the output image space Y is defined by some parametrization of these V-lines:

Vf (y) = V-line(y) f (x) dx.
(1.

3)

The inverse problem associated to this Radon transform is then defined by means of a filtered back-projection procedure.

Chapter 4 details the numerical implementation of both direct and inverse problems, the sampling details and simulations results for two test input objects representing some flat heritage sample.

A complete three-dimensional framework is then studied in Chapter 5 by means of the conical Radon transform C. This integral transform acts on functions defined on some domain of R 3 and integrate theses functions over conical surfaces with a fixed axis direction. The output image space is thus made of functions on a space defined by some parametrization of these cones: Cf (y) = cone surface(y) f (x) dx.

(1.4)

These manifolds are considered because the opening angle of the cone stands for a scattering angle in a Compton event. The model will take into account a realistic scattering scenario with the correspondent physical parameters, for both the direct and inverse problem.

Finally, Chapter 6 presents the associated discrete model for both the direct and inverse problem and numerical simulations results to validate the workability of the approach and to propose the development of the device.

Publications

Journal paper

Introduction

Le présent travail porte sur les aspects mathématiques d'une modalité d'imagerie par rayonnement dans le but de caractériser des matériaux anciens hétérogènes. Cette caractérisation est effectuée du point de vue de la densité électronique des objets. En effet, la connaissance de la distribution de la densité électronique interne d'un objet nous fournit des informations sur la morphologie et la composition chimique de l'objet qui va nous permettre d'accomplir une étude focalisée sur d'autres domaines de recherche, tels que les sciences des matériaux ou l'histoire de l'art [2]. Malheureusement, la géométrie aplatie des objets interdit certaines techniques d'imagerie bien étudiées comme la photo-absorption ou la tomographie par contraste de phase. La modalité d'imagerie par rayonnement utilisée ici est basée sur la tomographie de rayonnement Compton diffusé, développée ces dernières années principalement pour des applications dans les techniques d'imagerie par émission [START_REF] Morvidone | On the V-line radon transform and its imaging applications[END_REF][START_REF] Nguyen | Radon transforms on a class of cones with fixed axis direction[END_REF][START_REF] Truong | On new V-line radon transforms in R 2 and their inversion[END_REF][START_REF] Rigaud | Modeling and simulation results on a new compton scattering tomography modality[END_REF]. Plus précisément, il s'agit d'un système d'imagerie avec une source parallèle externe sans rotation relative permettant une cartographie tridimensionnelle des objets du patrimoine.

La modélisation des systèmes d'imagerie nécessite d'étudier à la fois le problème direct et inverse. Le problème direct est appelé le processus de formation d'image, qui vise à fournir un modèle de l'image générée à partir d'un objet d'entrée donné. Le problème inverse est appelé reconstruction d'objet, ou inversion, et fournit une approximation de la distribution d'une certain caractéristique physique inconnue de l'objet. Dans la tomographie de rayonnement Compton diffusé , l'image de sortie est une image spectrale résolue à la fois en énergie et en espace alors que l'objet reconstruite est une cartographie de la distribution de la densité électronique de l'objet d'entrée.

Après une étude théorique du problème direct et inverse, des simulations numériques ont été effectuées a partir d'objets d'entrée discrets représententatifs d'objets aplatis du patrimoine. L'image spectrale de sortie est obtenue grâce à une implémentation numérique du problème direct. Une reconstruction numérique de l'objet à partir de l'image simulée est ensuite effectuée en résolvant le problème inverse. La qualité de l'objet reconstruit est 11 finalement mesurée en calculant une certaine distance entre l'objet initial et sa reconstruction.

Deux objectifs majeurs peuvent être identifiés dans ce travail doctoral. Le premier est de développer un modèle mathématique du problème direct de la formation d'image. Ce modèle doit prendre en compte tous les paramètres rencontrés dans un scénario réaliste et ensuite être capable d'estimer numériquement pour un objet d'entrée donné l'image spectrale de sortie. Le processus de formation d'image sera comparé aux simulations Monte Carlo du transport de photons telles que définies dans [START_REF] Salvat | PENELOPE, a code system for Monte Carlo simulation of electron and photon transport[END_REF]. Le deuxième objectif principal de cette thèse est la définition d'une procédure de reconstruction pour le problème inverse. Une formule inverse analytique sera obtenue et discrétisée numériquement pour ainsi obtenir une reconstruction numérique de l'objet d'entrée. Une expérience synchrotron réelle où des données associées à la diffusion Compton seraient récupérées avec la configuration d'imagerie proposée est l'étape restante pour une validation de l'approche présentée. Ce travail a été publié partiellement dans [START_REF] Guerrero | Three-dimensional imaging of flat natural and cultural heritage objects by a compton scattering modality[END_REF][START_REF] Guerrero | Modelling of a New Compton Imaging Modality for an In-depth Characterisation of Flat Heritage Objects[END_REF].

Une description plus détaillée du contexte particulier des objets aux matériaux anciens et de la configuration expérimentale associée est proposée dans la section suivante. Cette introduction se termine par un bref aperçu des chapitres suivants.

Spécificités des études tridimensionnelles des matériaux anciens

Les objets aux matériaux anciens sont des vestiges du passé qui ont une valeur historique. De tels objets sont rencontrés en archéologie ou en paléontologie et appartiennent au patrimoine culturel. L'étude de leur composition chimique et de leur morphologie fournit non seulement des informations sur leur contexte historique (origine, utilisation ...) mais permet aussi de comprendre leurs processus d'altération et de proposer les stratégies les plus appropriées pour leur restauration et leur conservation. Une étude complète des caractéristiques spécifiques des matériaux anciens hétérogènes dans le contexte des expériences synchrotron est présentée dans [2]. Des modalités d'imagerie ou de spectroscopie, telles que la tomographie d'absorption, d'émission ou de contraste de phase parmi d'autres techniques d'imagerie à base de rayons X ont été développées pour accéder à cette information. Quand on travaille avec des objets du patrimoine, les propriétés non invasives et non destructives des inspections sont une exigence que la méthodologie d'imagerie par rayons X fournit, permettant l'imagerie bidimensionnelle de l'échantillon. Néanmoins, on peut facilement faire face à des échantillons présentant une géométrie aplatie, c'est-à-dire des échantillons présentant un rapport important entre leur surface frontale et leur épaisseur. Le défi consiste alors à effectuer une exploration tridimensionnelle sans utiliser de rotation relative entre l'échantillon et le système d'imagerie comme cela serait fait en tomographie conventionnelle, en utilisant la modalité d'absorption ou de contraste de phase.

Des échantillons présentant de telles caractéristiques se rencontrent notamment dans des études de conservation et de restauration de la peinture de chevalet. En effet, ils demandent une caractérisation de l'assemblage stratigraphique de pigments souvent sur une matrice de fond très dense telle que celle en plomb blanc. Un exemple de ce genre d'études est montré dans la figure 1.1 qui a été réalisée par une méthode invasive. La figure 1.2 présente un autre exemple d'objets possédant cette morphologie, à savoir en paléontologie les fossiles Lagerstätten [START_REF] Gueriau | Trace elemental imaging of rare earth elements discriminates tissues at microscale in flat fossils[END_REF] qui sont mécaniquement aplatis au cours du processus de fossilisation et se tiennent d'un côté d'une matrice sédimentaire épaisse qui ne peut pas être éclaircie pour l'étude de l'objet. Dans les deux cas, le volume d'intérêt forme une couche au-dessus d'un support matériel opaque aux rayons X, soit par sa densité, soit par son épaisseur.

Deux alternatives sont disponibles pour résoudre ce problème: soit effectuer une section stratigraphique de l'échantillon qui est une méthode invasive, soit limiter l'étude à une analyse bidimensionnelle de la surface frontale de l'échantillon, par exemple avec une technique basée sur la flourescence de rayons X synchrotron réalisée dans [START_REF] Gueriau | Trace elemental imaging of rare earth elements discriminates tissues at microscale in flat fossils[END_REF], ce qui laisse cependant une étude tridimensionnelle à réaliser.

Pour surmonter les limites des méthodologies mentionnées pour de tels échantillons, une modalité d'imagerie basée sur la tomographie de rayonnement Compton diffusé est proposée dans ce travail.

Principes généraux de la tomographie de rayonnement Compton diffusé

La tomographie classique aux rayons X considère principalement quatre phénomènes fondamentaux d'interaction photon-matière, à savoir l'absorption photoélectrique, la création de paires électron-positon, la diffusion de Rayleigh à la fois élastique et cohérente et la diffusion Compton qui est à la fois inélastique et incohérente. Dans la tomographie par rayons X, le signal Compton diffusé est considéré comme du bruit ajouté à l'absorption photoélectrique et aux données de diffusion cohérente dans le processus de formation d'image. En effet, le signal de transmission des rayons X est dominé par l'absorption photoélectrique tandis que la diffusion cohérente peut produire des variations d'amplitude importantes aux faibles angles de diffusion grâce aux effets d'interférence constructifs et destructifs dus à la nature cohérente de cette diffusion.

Dans la tomographie de rayonnement Compton diffusé (TCD), un photon incident d'énergie E 0 est absorbé par un électron cible, qui réémet un photon secondaire diffusé d'un angle ω par rapport à la direction du photon original (voir Figure 1.3). Le photon diffusé a une énergie E ω qui est liée à l'angle de diffusion ω par l'identité (1.1) appelée Introduction l'équation de Compton qui est donnée par

E ω = E 0 1 + E 0 m e c 2 (1 -cos ω) , (1.5)
où m e c 2 = 511 keV est l'énergie d'un électron au repos. En fonction du matériau, si le rayonnement incident a une énergie supérieure à environ 4 × 10 4 électrons-volts, la diffusion Compton devient le phénomène dominant dans le processus, d'autant plus lorsque la détection est effectuée en dehors de la zone de transmission directe (ω = 0). Dans ces conditions, en utilisant la TCD, nous pouvons éviter les rotations relatives entre l'échantillon et le système d'imagerie.

D'autre part, en raison de la potentielle présence d'un matériau de support opaque, les données de transmission et de diffusion Compton vers l'avant, c'est-à-dire avec un angle de diffusion inférieur à π 2 , sont impossibles à collecter. Ceci motive la proposition d'une modalité basée sur des données rétro-diffusées, c'est-à-dire des données collectées avec un angle de diffusion ω compris entre π 2 et π (voir figure 1.5). L'idée d'exploiter le rayonnement diffusé par l'effet Compton dans les techniques d'imagerie a été introduite et étudiée simultanément. Elle a donné naissance à la TCD [6,[START_REF] Morvidone | On the V-line radon transform and its imaging applications[END_REF][START_REF] Nguyen | Radon transforms on a class of cones with fixed axis direction[END_REF][START_REF] Norton | Compton scattering tomography[END_REF][START_REF] Truong | On new V-line radon transforms in R 2 and their inversion[END_REF], dont le but est la reconstruction de la carte de densité électronique de l'objet. Le premier TCD scanner a été proposé en 1994 par Norton [START_REF] Norton | Compton scattering tomography[END_REF] à travers une transformation de Radon sur des arcs de cercle partant d'une source γ vers un site de détection en mouvement. La transformation de Radon sur des surfaces coniques ayant des directions axiales fixes et des angles d'ouverture variables a été proposé et étudié dans [START_REF] Nguyen | Radon transforms on a class of cones with fixed axis direction[END_REF] où une première formule d'inversion analytique est proposée en utilisant l'analyse en composants circulaires, avec des applications en imagerie d'émission. Des généralisations de ce type de transformations vers des espaces de dimensions supérieures avec des formules inverses associés du type rétroprojection filtrée sont présentées dans [START_REF] Haltmeier | Exact reconstruction formulas for a radon transform over cones[END_REF]. Un algorithme d'inversion par rétroprojection filtrée d'une transformée de Radon conique dans R 3 a été développé récemment dans [6].

La configuration du système d'imagerie

La TCD vise à reconstruire la carte de la densité électronique de l'objet étudié. Dans ce travail, la densité électronique est représentée mathématiquement par une fonction non négative définie sur l'espace de Schwartz des fonctions lisses et à décroissance rapide sur R 2 .

Une configuration synchrotron avec un faisceau parallèle de rayons X monochromatique (environ 50 keV) et un détecteur résolu spectralement et spatialement est considéré. Comme présenté sur la figure 1

Bref résumé des chapitres 2 à 6

La transformée de Radon modélise le processus de formation d'image en tomographie de photo-absorption classique. Dans le cas bidimensionnel le plus simple, une fonction continue f appartenant à l'espace des objets d'entrée X des fonctions définies sur R 2 est donnée. La réponse de cette transformation intégrale est une fonction Rf définie dans l'espace des images de sortie Y qui modélise la sortie du système d'imagerie. Cette fonction de sortie effectue des intégrations sur des lignes dans R 2 de la fonction d'entrée suivant une direction de rayonnement donnée y qui caractérise la ligne. Chaque intégrale curviligne Introduction est appelée une projection de l'objet, et peut s'écrire comme

Rf (y) = y f (x) dx.
(1.6)

La résolution de cette dernière équation pour f est le sujet principal du problème inverse, conjointement à une implémentation numérique réalisable à partir d'un ensemble de données de taille importante.

Le chapitre 2 rappelle le problème direct et inverse de la transformation de Radon originale, couvrant le cadre continu et discret. La transformée du Radon en coordonnées cartésiennes est étudiée comme une étape préliminaire à la transformation en demi-espace de Radon, qui devient elle-même une étape préliminaire à la transformée de Radon en V exposée au chapitre 3. Pour la transformée cartésienne, les fonctions de m'image de sortie sont définies dans l'espace cartésien Y = R 2 . La transformée en demi-espace de Radon a le même espace d'image de sortie Y mais l'espace d'objet d'entrée est limité aux fonctions définies sur R + × R. Les formulations continues et discrètes du problème direct et inverse sont traitées et illustrées par des simulations numériques pour chaque transformée de Radon.

Diverses généralisations de la transformée de Radon ont été étudiées à partir du travail initial de Radon. Elles sont principalement appliquées dans les sciences de l'imagerie où les intégrations sont réalisées sur des différentes variétés selon la configuration d'imagerie. La modalité de diffusion Compton proposée dans cette thèse est d'abord limitée à un cadre bidimensionnel pour une première compréhension du modèle dans le chapitre 3. Dans ce cas, les intégrations de la fonction d'entrée doivent être effectuées sur deux demi-droites de sommet commun qui forment une ligne brisée d'allure V dans R 2 . La transformée de Radon en V,V, est appliquée sur les fonctions définies sur R + ×R alors que l'espace d'image de sortie Y est défini par une certaine paramétrisation de ces lignes en V:

Vf (y) = V-ligne(y) f (x) dx.
(1.7)

Le problème inverse associé à cette transformation de Radon est ensuite défini avec une procédure de rétroprojection filtrée.

Le chapitre 4 détaille l'implémentation numérique des problèmes directs et inverses, les détails de l'échantillonnage et les résultats des simulations pour deux objets d'entrée représentant un échantillon plat du patrimoine.

Un cadre tridimensionnel complet est ensuite étudié dans le chapitre 5 basé sur la transformée de Radon conique C. Cette transformée intégrale agit sur des fonctions définies sur un domaine de R 3 et intègre ces fonctions sur des surfaces coniques avec une direction axiale fixe. L'espace d'image de sortie est donc constitué de fonctions sur un espace défini par une certaine paramétrisation de ces surfaces coniques:

Cf (y) = surface conique(y) f (x) dX.
(1.8)

Ces variétés sont considérés sachant que l'angle d'ouverture du cône représente un angle de diffusion dans un événement Compton. Le modèle prendra en compte un scénario de diffusion réaliste avec les paramètres physiques correspondants, aussi bien pour le problème direct que pour le problème inverse.

Finalement, le chapitre 6 présente le modèle discret associé ainsi que des résultats de simulations numériques de problèmes directs et inverses afin de valider la faisabilité de l'approche.

2

Recall of the forward and inverse standard Radon transforms

Introduction

The Radon transform is an integral operator that maps a function defined on R n into the set of its integrals over the hyperplanes of R n [START_REF] Natterer | Mathematical Methods in Image Reconstruction[END_REF]. It was first studied by H.A. Lorentz who gave a solution, i.e., an inverse formula, for the case n = 3. Although it was not published and it is not known either the motivation of his work or which techniques he used, it was Later on, the Radon transform was formally introduced by the Austrian mathematician Johann Radon in his 1917 paper [START_REF] Radon | Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten[END_REF] where he investigates the problem of recovering a function of two dimensions from its line integrals. The article includes in particular a solution for the inverse problem when n = 2. An english translation of this paper is found in [START_REF] Deans | The Radon Transform and Some of its Applications[END_REF]. Although the transform was studied only from a theoretical point of view, it became afterwards a key tool in applied mathematics, namely in tomography. Since Radon's paper, there have been various rediscoveries of the transform. A. Cormack gave a remarkable solution to the bi-dimensional case and applied it on tomography in 1963 [START_REF] Cormack | Representation of a function by its line integrals, with some radiological applications[END_REF]. He used Fourier techniques to invert the transform and proposed to apply the solution to reconstruct the absorption coefficient of an object using a X-ray scanner.

In 1972 G. Hounsfield developed an imaging system based on parallel X-rays. This gave birth to the so-called X-ray computed tomography. It uses the Radon transform and Cormack's inverse procedure to produce medical diagnostics at first. The principle is to recover the attenuation coefficient of a studied object. To achieve this, a parallel X-ray source radiates the object and the non-absorbed radiation that cross the object is collected by X-ray detectors, as illustrated in Figure 2.1. The procedure is repeated by rotating the imaging system or the object to obtain the data or projections of the object, corresponding to the Radon transform. The target attenuation coefficient is then recovered by an inversion procedure of the transform.

A complete review of the history of the Radon transform is exposed in [START_REF] Cormack | Computed tomography: some history and recent developments[END_REF]. For deeper reviews about the Radon transform see references [START_REF] Deans | The Radon Transform and Some of its Applications[END_REF][START_REF] Helgason | The Radon Transform[END_REF][START_REF] Natterer | Mathematical Methods in Image Reconstruction[END_REF] whereas a good summary can be found in [START_REF] Rigaud | Principe de la tomographie 2D par transmission et transformation de Radon[END_REF] and recent developments or generalisations for instance in [START_REF] Rigaud | Study of Generalized Radon Transforms and Applications in Compton Scattering Tomography[END_REF][START_REF] Truong | New properties of the V-line radon transform and their imaging applications[END_REF]. Numerical implementation details are presented in [START_REF] Toft | The Radon Transform: Theory and Implementation[END_REF].

In this chapter, the original Radon transform is first recalled. Then, an equivalent definition is derived, called the Radon transform in cartesian coordinates. Finally, the halfspace Radon transform is introduced as a preliminary step to the V-line Radon transform that will be applied in Compton scattering tomography in section 3.1.1.

The main inversion techniques for these three Radon transforms are also presented, all based on Fourier transforms. The chapter ends with simulations for each cases, illustrating the X-ray computed tomography applications.

The original Radon transform

Radon's original definition of the transform introduced in [START_REF] Radon | Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten[END_REF] can be formulated as a transform R mapping functions from the Schwartz space S(R 2 ) to S(R × S 1 ) where S 1 represents here the one-dimensional unit sphere. A general definition related to functions defined on R n is first given below.

Definition 1. Let f : R n → R + be a S(R n )-function, (p, n) ∈ R × S n-1 and H(p, n) = {x ∈ R n : p = x • n}. (2.1)
the affine hyperplane perpendicular to the vector pn.

The Radon transform of f is the function Rf defined for (p, n) as the integral of f over the hyperplane H(p, n), i.e.,

Rf (p, n) = H(p,n) f (x) dx. (2.2)
Or it is equivalently defined by means of the Dirac δ function as

Rf (p, n) = R n f (x)δ(p -x • n) dx. (2.
3)

It is convenient to let the function f belong to a nice class of functions such as the Schwartz space S(R n ) of rapidly decreasing C ∞ functions. Among the advantages of working with this space, there is the fact that the function Rf can be differentiated as often as desired and the possible application of the Fubini's theorem [START_REF] Deans | The Radon Transform and Some of its Applications[END_REF]. A more general definition is possible, for instance in the space L 2 loc (R n ).

Remark 1. The Radon transform is a linear mapping

R : S(R n ) → S(R × S n-1 ) f → Rf. (2.4)
For a proof of the last assertion, see [START_REF] Helgason | The Radon Transform[END_REF].

Remark 2. Since H(p, n) = H(-p, -n), the mapping (p, n) → H(p, n) is a double cover- ing of R × S n-1
onto the space of hyperplanes in R n . Hence, we can consider p ∈ R + or n in a half of the unit sphere S n-1 .

For each couple (p, n), Rf (p, n) is called a projection of f over the hyperplane H(p, n), and for all (p, n) ∈ R × S n-1 , Rf (p, n) is called the image formed from the function or object f . This process is known as image formation or the direct problem of the Radon transform.

The obvious question is how to recover f (x), for all x ∈ R n , given the projections Rf (p, n). Such inverse process is known as object reconstruction. This section deals with both problems for the two-dimensional (2D) case.

The 2D direct problem

This section gives an outline of the direct formulation, i.e., the study of the resulting function Rf for a given function f in the two dimensional case. The 2D framework is particularly relevant due to its applications in computed tomography. Following previous notations when n = 2, x belongs here to R 2 and n to the unit sphere S 1 . Let us introduce the notations, for ϕ ∈ [0, 2π[ x = (x, y), n = (cos ϕ, sin ϕ).

(2.5)

The affine hyperplane perpendicular to pn is in this case the line of parameters (p, ϕ) defined by L(p, ϕ) = {(x, y) ∈ R 2 : p = x cos ϕ + y sin ϕ}.

(2.6)

Definition 2. The two-dimensional Radon transform, or the (p, ϕ)-Radon transform of a function f ∈ S(R 2 ), noted f , is given by

f (p, ϕ) = L(p,ϕ) f (x, y) dl, (2.7)
where dl is the measure along the line L(p, ϕ). The geometry of the (p, ϕ)-Radon transform is represented in Figure 2.2.

As in R n , it can also be defined through the Dirac δ function as

f (p, ϕ) = R 2 f (x, y) δ(p -x cos ϕ -y sin ϕ) dxdy. (2.8)
This (p, ϕ)-Radon transform can also be written in a more explicit form, by means of a change of variables or an axis rotation by ϕ via x = p cos ϕ -r sin ϕ, y = p sin ϕ + r cos ϕ, dr = dl.

(2.9)

Then, (2.8) can be rewritten as

f (p, ϕ) = R f (p cos ϕ -r sin ϕ, p sin ϕ + r cos ϕ) dr. (2.10)
This last integral is more adapted to perform numerical simulations of the Radon transform.

Remark 3. To let L(p, ϕ) cover the space of affine hyperplanes of R 2 , as pointed out in Remark 2, ϕ only needs to vary in [0, π[ supposing we still have p ∈ R.

The 2D inverse problem

The goal of this section is to solve (2.7) for f (x, y), i.e., to reconstruct an object represented by the function f from the measured projections f (p, ϕ). Note that to cover the topics treated in this work, we only need to focus on the inversion of the two-dimensional Radon transform (see [START_REF] Deans | The Radon Transform and Some of its Applications[END_REF][START_REF] Helgason | The Radon Transform[END_REF] for a general case).

2.2.2.a A first reconstruction: the Fourier slice theorem

The crucial connection between the Radon transform and the Fourier transform has been known with no doubt shortly after the discovery of the first one. But the most contributory connection was brought to light by Bracewell in 1956 [START_REF] Bracewell | Strip integration in radio astronomy[END_REF]. It is the Fourier slice theorem, that states the relationship between the 1D Fourier transform of a projection and the 2D Fourier transform of the original object.

In the following, let F n be the n-dimensional Fourier transform operator and for a given angle ϕ, let fϕ be the 1D function fϕ : p → f (p, ϕ).

(2.11)

The Fourier slice theorem relates the 1D Fourier transform of a projection for a given angle F fϕ (q) with a line slice of the 2D Fourier transform of the object F 2 f (u, v) subtending an angle ϕ with the u-axis, i.e., a line of equation v cos ϕ = u sin ϕ.

Theorem 1 (Fourier slice theorem). Let f ∈ S(R 2 ) and ϕ ∈ [0, π[, then we have

F fϕ (q) = F 2 f (q cos ϕ, q sin ϕ), (2.12
)

or reciprocally, for (u, v) ∈ R 2 F 2 f (u, v) = F fϕ(u,v) (||(u, v)||), (2.13)
where ϕ(u, v) and ||(u, v)|| are respectively the angle and the Euclidean norm of (u, v).

Proof. The 1D Fourier transform of fϕ reads

F fϕ (q) = ∞ -∞ fϕ (p)e -2πipq dp = ∞ -∞
f (p, ϕ)e -2πipq dp.

(2.14)

Using definition (2.8) and applying Fubini's theorem, we get from the last integral

F fϕ (q) = R 2 f (x, y) ∞ -∞ e -2πipq δ(p -x cos ϕ -y sin ϕ) dp dxdy, (2.15)
where the last inner p-integration can be easily solved and gives

F fϕ (q) = R 2
f (x, y)e -2πi(xq cos ϕ+yq sin ϕ) dxdy.

(2.16)

The result claims by definition of the 2D Fourier transform.

This result gives us a first inversion of the Radon transform using Fourier transforms. We actually have the 2D Fourier transform of f on radial lines v cos ϕ = u sin ϕ. The inverse formula is announced in the following corollary.

Corollary 1 (inverse formula). The 2D Fourier transform of the object F 2 f (u, v) is available through the Fourier slice theorem for all (u, v) ∈ R 2 . Therefore, the object f is recovered by the 2D inverse Fourier transform, that is by the inverse formula

f (x, y) = R 2 F fϕ(u,v) (||(u, v)||)e 2πi(ux+vy) dudv.
(2.17)

This reconstruction method is often referred as the direct Fourier method. Numerically, it requires a 2D interpolation in the Fourier domain and computations of 2D Fourier transforms.

2.2.2.b A second reconstruction: the filtered back-projection

The filtered back-projection reconstruction method, detailed in [START_REF] Deans | The Radon Transform and Some of its Applications[END_REF][START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF][START_REF] Natterer | Mathematical Methods in Image Reconstruction[END_REF] was introduced due to its numerical advantages compared to the Fourier slice theorem. It only involves 1D Fourier transforms and 1D (linear) interpolations. It is the most important reconstruction algorithm in tomography [START_REF] Natterer | Mathematical Methods in Image Reconstruction[END_REF].

Basically, it consists in modifying, or filtering the projections f of a target object in the Fourier domain, then computing the resulting filtered projections f * and finally to back-project them by means of the adjoint operator R † defined in this section.

We first need to introduce the adjoint Radon transform, concept related to backprojecions.

The adjoint transform. In the following, let us consider the notation

X = S(R 2 ), Y = S(R × [0, π[), (2.18)
with respective inner products

(f 1 , f 2 ) X = R 2 f 1 (x, y)f 2 (x, y) dxdy, (g 1 , g 2 ) Y = π 0 R g 1 (p, ϕ)g 2 (p, ϕ) dpdϕ.
(2.19)

Definition 3. The adjoint to the Radon transform R † is the unique operator

R † : Y → X g → R † g, (2.20)
that verifies, ∀f ∈ X, ∀g ∈ Y, (Rf, g) Y = (f, R † g) X .

(2.21)

In order to derive an expression for R † , we start from the left hand side of the previous identity

(Rf, g) Y = π 0 R Rf (p, ϕ)g(p, ϕ) dpdϕ, (2.22)
and insert definition (2.8) together with Fubini's theorem to have

(Rf, g) Y = R 2 f (x, y) π 0 R g(p, ϕ)δ(p -x cos ϕ -y sin ϕ) dpdϕ dxdy. (2.23)
Finally, inner p-integration is performed and we have

(Rf, g) Y = R 2 f (x, y) π 0 g(x cos ϕ + y sin ϕ, ϕ) dϕ dxdy, (2.24)
from where one is able to extract the adjoint transform R † of R that verifies (2.21). R † is therefore defined by

R † g(x, y) = π 0 g(x cos ϕ + y sin ϕ, ϕ) dϕ. (2.25)
In the way it is defined, the adjoint transform can be interpreted as a back-projection procedure of projections at the same angle ϕ. It will be clarified subsequently.

Theorem 2 (Filtered back-projection). Under same hypotheses of Theorem 1, define the filtered projections of f by

f * (p, ϕ) = F -1 1 |q|F fϕ (q) (p), (2.26
)
where fϕ is defined in (2.11). Then, f is recovered by means of identity the f = R † f * , i.e., by

f (x, y) = π 0 f * (x cos ϕ + y sin ϕ, ϕ) dϕ. (2.27)
Proof. The Fourier slice theorem plays a key role in this proof. Let us express the function f through its 2D Fourier transform:

f (x, y) = R 2
F 2 f (u, v)e 2πi(ux+vy) dudv.

(2.28)

The 2D Fourier transform of f can be replaced by the 1D Fourier transform of the projections thanks to (2.13) and then

f (x, y) = R 2 F fϕ(u,v) (|(u, v)|)e 2πi(ux+vy) dudv. (2.29)
Applying a change of variables to polar coordinates:

u = q cos ϕ, v = q sin ϕ, dudv = q dqdϕ, (2.30) equation (2.29) becomes f (x, y) = 2π 0 ∞ 0
F fϕ (q)e 2πiq(x cos ϕ+y sin ϕ) q dqdϕ. (2.31)

In order to make appear an inverse Fourier transform and the adjoint operator in the last double integral, the limits of integration can be adapted as

ϕ ∈ [0, π[, q ∈ R. (2.32)
due to the fact that

F fq (ϕ) = F f-q (ϕ + π), (2.33)
Then (2.31) reads

f (x, y) = π 0 R
|q|F fϕ (q)e 2πiq(x cos ϕ+y sin ϕ) dqdϕ.

(2.34)

The q-integration corresponds to the inverse Fourier transform of the function q → |q|F fϕ (q) (2.35) evaluated at x cos ϕ + y sin ϕ, i.e., (2.41)

f (x, y) = π 0 F -1 1 |q|F fϕ (q) (x cos ϕ + y sin ϕ) dϕ. ( 2
This last equation give rise to the inversion formula found by Radon in [START_REF] Radon | Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten[END_REF], by different means, that reads

f = - 1 2π R † H ∂ ∂p f .
(2.42)

The Radon transform in cartesian coordinates

By means of a change of variables, the Radon transform in Cartesian coordinates can be formulated, that is, a transform R mapping from the Schwartz space S(R 2 ) to S(R 2 ). This version of the transform is applied in fields like seismology or slant stacking. It is numerically interesting because it allows the replacement of the computationally expensive back-projection by a series of fast Fourier transforms [START_REF] Edholm | The linogram algorithm and direct Fourier method with linograms[END_REF][START_REF] Truong | New properties of the V-line radon transform and their imaging applications[END_REF]. The associated direct and inverse problem of the transform are studied in this section.

The direct problem

(u, v) ∈ R 2 via u = p cos ϕ , v = tan ϕ which implies 1 + v 2 = 1 cos ϕ . (2.43)
The geometrical representation of the new variables (u, v) of the transform is depicted in Figure 2.3.

With these new variables, (2.8) can be rewritten as

f (u, v) = 1 + v 2 R 2 f (x, y) δ(u -x -vy) dxdy, (2.44)
and after performing x-integration:

f (u, v) = 1 + v 2 ∞ -∞
f (u -vy, y) dy.

(2.45)

Note the use of the same notation f for sake of simplicity.

Remark 6. The cartesian Radon transform is now a linear mapping

S(R 2 ) → S(R 2 ) f → f . (2.46)

The inverse problem

As in the (p, ϕ)-Radon transform, we write in this section the inversion procedure for the cartesian Radon transform starting with an inverse formula derived from the Fourier slice theorem. Then, a filtered back-projection expression type is derived.

2.3.2.a A first inversion formula

Fourier techniques are used here to write an inversion formula for (2.45). The Fourier slice theorem corresponding to the cartesian Radon transform is the following.

Theorem 3 (Fourier slice theorem). Let f ∈ S(R 2 ) , µ ∈ R and denote Let F a f be the 1D Fourier transform of f according to the variable a. Then, the following identity is verified

F u f (µ, v) = 1 + v 2 F 2 f (µ, vµ).
(2.47)

Proof. Let f (x, y) be expressed using its 1D x-Fourier transform,

F x f : µ → f (µ, •), as f (x, y) = ∞ -∞ F x f (µ,
(u, v) = 1 + v 2 R 2
F x f (µ, y)e 2πiµ(u-vy) dµ dy, (2.49) and applying Fubini's theorem,

f (u, v) = 1 + v 2 ∞ -∞ e 2πiuµ ∞ -∞ F x f (µ, y)e -2πiyvµ dy dµ. (2.50)
The last inner integral is the 1D y-Fourier transform of the function y → F x f (•, y) evaluated at vµ, i.e., it is F y (F x f )(µ, vµ). Thus, (2.50) becomes

f (u, v) = 1 + v 2 ∞ -∞ F y (F x f )(µ, vµ)e 2πiuµ dµ = 1 + v 2 ∞ -∞ F 2 f (µ, vµ)e 2πiuµ dµ.
(2.51) From where, by definition of the 1D u-Fourier transform of f , the result of the theorem claims.

From (2.47), we can apply the inverse and scaling property of the Fourier transform to get the inverse formula of (2.45). This result is announced in the following corollary. x+vy) dvdµ.

Corollary 2. If f is the cartesian Radon transform of a function f ∈ S(R 2 ), f is then recovered via the formula f (x, y) = ∞ -∞ |µ| ∞ -∞ 1 √ 1 + v 2 F u f (µ, v)e 2πiµ(
(2.52)

2.3.2.b Inversion via filtered back-projection

Back-projections reconstructions can be adapted to Radon transforms defined over different manifolds. In all cases, the back-projection process is obtained via the associated adjoint transform.

In the present case of the cartesian Radon transform, the adjoint transform is introduced and then a back-projection inversion is presented below. Definition 4. The adjoint transform of the cartesian Radon transform is defined in the same way as in (2.21). It is then easy to verify that it is given by

R † g(x, y) = R 1 + v 2 g(x + yv, v) dv.
(2.53)

From the inverse formula (2.52), we can apply apply Fubini's theorem to have a filtered back-projection inverse with the ramp filter |µ| in the Fourier domain, namely x+vy) dµ.

f (x, y) = ∞ -∞ dv √ 1 + v 2 ∞ -∞ |µ|F u f (µ, v)e 2πiµ(
(2.54)

The previous inverse formula also writes

f (x, y) = ∞ -∞ dv √ 1 + v 2 F -1 u (|µ|F u f )(x + vy, v). (2.55)
This last formula can also be expressed in terms of the adjoint transform after introducing the filtered projections

f * (u, v) = F -1 u (|µ|F u f )(u, v). (2.56)
Theorem 4 (filtered back-projection for the Radon transfom in Cartesian coordinates).

A function f ∈ S(R 2 ) can be recovered via

f = R † ( 1 1 + v 2 f * (u, v)), or equivalently f (x, y) = R 1 √ 1 + v 2 f * (x + yv, v) dv.
(2.57)

The half-space Radon transform

The half-space Radon transform and its analytic inverse formula were introduced recently in [START_REF] Truong | New properties of the V-line radon transform and their imaging applications[END_REF]. It is a Radon transform defined on half-lines, starting on points laying over a fixed line (here the line x = 0) towards the direction of where an object is placed.

The direct problem

The half space Radon transform acts on functions having support in R + × R and has the following definition.

Hf (ω, ζ) = ∞ 0 f (r cos ω, ζ -r sin ω) dr, (2.58)
where dr is the measure on the half-line.

This transform acts as the integral of a function over all points laying on the half-line starting at (0, ζ) towards the half-plane x > 0, and having a slope equal to tan(π -ω). This half-line is noted L ω,ζ and is represented in Figure 2.4.

Note that if this transform is applied on a function defined in R 2 , the reconstruction that will be defined below will only recover the part of the function located in the half space x > 0.

Through the following change of variables

x = r cos ω and t = tan ω, (2.59)
the transform can also be defined in terms of the Dirac δ function as the transform

S(R + × R) → S(R 2 ) f → Hf, (2.60)
by

Hf (t, ζ) = 1 + t 2 ∞ 0 f (x, ζ -tx) dx = 1 + t 2 R×R + f (x, y) δ(y -ζ + tx) dxdy.
(2.61)

Remark 7.
The half space transform is an intermediary step before defining the V-line Radon transform to be done on the next chapter. It is defined for functions on R + × R but can also be written through the Radon transform on R 2 (cartesian or the (p, ϕ)-transform) via the Heaviside step function H, restricting the definition domain to R + × R. In other words, the half-space Radon transform also writes as

f (p, ϕ) = R 2 f (x, y)H(x) δ(p -x cos ϕ -y sin ϕ) dxdy.
(2.62)

The inverse problem

In order to obtain an inverse formula for (2.58), we still use Fourier techniques. An analytical inverse formula and then an inversion via a filtered back-projection procedure are presented below.

2.4.2.a A first inversion formula

The Fourier slice theorem relates the 1D ζ-Fourier transform of H with a slice of the 2D Fourier transform of f given by the equation x = ty.

Theorem 5 (Fourier slice theorem). Let f ∈ S(R + × R), q ∈ R, then

F ζ Hf (t, q) = 1 + t 2 F 2 f (qt, q), (2.63)
or equivalently, with the angular variable w:

cos ωF ζ Hf (ω, q) = F 2 f (q tan ω, q). ( 2

.64)

We are then able to recover the function f via the inverse formula

f (x, y) = F -1 2 |q| √ 1 + t 2 F ζ Hf (t, q) (qx,
F ζ Hf (t, q) = R Hf (t, ζ)e -2πiqζ dζ = 1 + t 2 R×R + f (x, y) R δ(y -ζ + tx)e -2πiqζ dζ dxdy, (2.66)
from where the result claims by definition of the δ function as

F ζ Hf (t, q) = 1 + t 2 R×R + f (x, y)e -2πiq(tx+y) dxdy = 1 + t 2 F 2 f (qt, q).
(2.67)

2.4.2.b Inversion via filtered back-projection

As the classical Radon transform on straight lines, the filtered back-projection inversion on a half-space combines a back-projection operation and a filtering operation. It says that the back-projection b ω (x, y) at (x, y) for an angle ω assigns to (x, y) the value of the projection recorded at ζ = y + x tan ω where (x, y) was projected. That is to say, the back-projection of (x, y) at ω is given by

b ω (x, y) = Hf (ω, y + x tan ω), (2.68)
and equivalently, with t = tan ω as b t (x, y) = Hf (t, y + tx).

(2.69)

The filtering operation is applied to the ζ-Fourier transform of Hf (t, ζ). To obtain the correct filter, the inversion formula of (2.61) can be derived from (2.65) and Fubini's theorem as

f (x, y) = ∞ -∞ ∞ -∞ |q|F ζ Hf (t, q)e 2πiq(tx+y) dq dt √ 1 + t 2 , ( 2.70)
The inner integral can be identified as the inverse Fourier transform of the function q → |q|F ζ Hf (•, q) evaluated at y

+ xt f (x, y) = ∞ -∞ F -1 |q|F ζ Hf (t, q) (t, y + xt) dt √ 1 + t 2 . (2.71)
Equivalently, in terms of the angle ω, the last inversion formula writes:

f (x, y) = π/2 -π/2 F -1 (|q|F ζ Hf )(ω, y + x tan ω) dω cos ω . (2.72)
Using as before the ramp filter q → |q|, the filtered projections are noted H * f and defined by

H * f (t, ζ) = F -1 [|q|F ζ Hf ( t, q)](t, ζ). (2.73)
Then, one is able to write the back-projection of (x, y) related to filtered projections as bt (x, y) = H * f (t, y + xt).

(2.74)

A filtered back-projection reconstruction of f (x, y) for all (x, y) is also expressed in the following compact form: Theorem 6. For all f ∈ S(R + × R), one has:

f (x, y) = ∞ -∞ H * f (t, y + xt) dt √ 1 + t 2 .
(2.75) Remark 8. The last reconstruction can also be written with the angular parametrization as

f (x, y) = π/2 -π/2 H * f (ω, y + x tan ω) dω cos ω .
(2.76)

Point Spread Function

The Point Spread Function (PSF) or the impulse response characterizes the output of an integral transform. More precisely, for (x 0 , y 0 ) ∈ R 2 , the PSF is defined as the response of the Radon transform applied to the Dirac function

f 0 (x, y) = δ(x -x 0 , y -y 0 ). (2.77)
Due to the linearity of the transform, the PSF gives us an idea of the general response for a given input function.

Remark 9. Until now, the Radon transform has been defined for functions on S(R 2 ) and the Dirac delta function doesn't belong to the mentioned space as requested. Nevertheless, the Radon transform of distributions can also be defined in a week sense [START_REF] Helgason | The Radon Transform[END_REF] by using identity (2.21). In particular, the integral (2.8) when f is a Dirac function does exist in terms of distributions.

As mentioned before, the PSF characterizes a linear transform and then its support suggests us the shape of a corresponding output image. In this section, PSF are computed for (x 0 , y 0 ) ∈ R 2 (with x 0 > 0) for the three Radon transforms described above. The corresponding PSF supports are plotted in Figure 2.5, correspond to a case where (x 0 , y 0 ) = (1, 1). Regarding the (p, ϕ)-Radon transform, the PSF support has a sinusoidal shape. Indeed, from (2.8), the PSF at (x 0 , y 0 ), denoted f0 , is given by f0 (p, ϕ) = δ(p -x 0 cos ϕ -y 0 sin ϕ), (2.78) which has a sinusoidal shape support in the Radon space (p, ϕ), given by the equation

p = x 0 cos ϕ + y 0 sin ϕ. (2.79)
This output image is called a sinogram, and is represented in Figure 2.5a. Note that the ϕ axis of sinogram is located vertically following most of the literature.

In the case of the cartesian transform, by using (2.44), the PSF related to the point object

f 0 writes as f0 (u, v) = 1 + v 2 δ(u -x 0 -vy 0 ). ( 2

.80)

This PSF has a rectilinear shape support of equation u = x 0 + vy 0 and is called a linogram (see the case (x 0 , y 0 ) in Figure 2.5b).

Concerning the half-space transform, the PSF support related to this transform is obtained from expression (2.61) either in terms of t or ω:

Hf 0 (t, ζ) = 1 + t 2 δ(y 0 -ζ + tx 0 ), (2.81) Hf 0 (ω, ζ) = 1 + tan 2 ω δ(y 0 -ζ + x 0 tan ω).
(2.82)

If x 0 > 0, the first support is a line (as expected) having a slope t. For the second expression, the support is given by the equation

ζ = y 0 + x 0 tan ω, (2.83)
and is represented in Figure 2.5c.

Some simulations

The Shepp-Logan phantom [START_REF] Shepp | The fourier reconstruction of a head section[END_REF], that represents a human head, is commonly used to test reconstruction algorithms. It is illustrated in figure 2.6 and it will be used to expose reconstructions by means of the filtered back-projection method. It has dimensions of 240 × 240 pixels where each pixel value varies between 0 and 1 representing a target feature of the object. Simulations are carried out with the Radon transform in (p, ϕ)-coordinates, in cartesian coordinates and finally by using the half-space Radon transform. Reconstructions are performed with the filtered back-projection formulas given in (2.27), (2.57), (2.76) respectively. Correlations between the original Shepp-Logan phantom and its different reconstructions will be computed to give a numerical value of their quality. Simulations are implemented in the R-language [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF]. As n increases, the effect of the filter does the same and may over-smooth reconstructions.

Window functions

A so-called window function or apodization function q → w(q) is applied to the ramp filter q → |q| in numerical reconstructions to control high frequencies due to the divergence of the ramp filter and to attenuate noise in reconstructions.

In such case, filtered projections written in (2.26), (2.56) and (2.73) for the (p, ϕ)-Radon, cartesian and half space Radon transforms are rewritten to be numerically implemented respectively in the form

f * (p, ϕ) = F -1 1 |q|w(q)F fϕ (q) (p), (2.84) f * (u, v) = F -1 u |µ| w(µ)F u f (u, v), (2.85)
H * f (ω, ζ) = F -1 |q|w(q) Hf (ω, ζ). (2.86)
Cosine functions are widely used as window functions:

w(q) = cos n (πq).

(2.87)

The behaviour of this function according to n is shown in Figure 2.7. As the value of n increases, reconstructions are smoother due to the apodization reducing for example the presence of noise. Important details in object reconstructions may also disappear if n is too large so the choice of n depends on the imaging configuration and expected results. In other words, it is a compromise between regularity and resolution. For the numerical simulations presented in this section, in the absence of noise, the value n = 1 will be chosen. Other window functions are presented in literature such as Hann, Hamming, Shepp Logan (see [START_REF] Easton | Fourier Methods in Imaging[END_REF] for more details about window functions).

The discrete formulation

Up to now, the direct and inverse problem have been worked out on a continuous framework, that is, the object function is a continuous function defined on R 2 and the related projections are also continuous functions. However, imaging systems collects only a finite number of discrete projections that will consist on the data. Therefore, reconstruction algorithms recover the input function as a discrete and finite object defined on some grid. Discrete versions of the direct and inverse problem must be derived.

2.6.2.a Sampling details

In simulation works, test objects as the Shepp-Logan phantom are known on a finite and discrete array of size N × M represented by the set

{f (x i , y j ), 0 ≤ i < M, 0 ≤ j < N }.
(2.88)

Each (x i , y j ) represents a pixel position in the object domain having dimensions ∆x and ∆y. In the case of the Shepp-Logan phantom, we have N = M = 240 and ∆x = ∆y = 1.

Image reconstructions will also be a discrete array of same dimensions as the input object.

The image projection, discrete as well, will be an array which size depends on the number of detectors positions and the number of rotations made. Consequently, the direct problem modelled by equations (2.10) for the (p, ϕ)-Radon transform, (2.45) for the cartesian case and (2.58) for the half-line transform are no longer integrals but finite sums performing numerical integrations. In presents simulations, the trapezoidal rule is used to implement those integrals with an uniform sampling step ∆r. Discretization parameters are given in what follows.

Integrals (2.10), (2.45) and (2.58) for the presented transforms request sampled values that may not coincide with values on the grid (2.88). Interpolation of the discrete object f is then needed. The most common methods to face this problem are, from the less to the most computationally expensive, the nearest neighbour, linear and sinc interpolation. Details and comparisons of each approach applied to the Radon transform can be found in [START_REF] Toft | The Radon Transform: Theory and Implementation[END_REF]. In present simulations, linear interpolation is chosen due to the fact that the Shepp Logan phantom does not present strong sharp edges.

The inverse problem is solved with the filtered back-projection procedure. The method also involves discrete Fourier transforms, implemented with the fast Fourier transform algorithm. Back-projections need to be performed from the discrete filtered data by means of the adjoint transform appearing in inverse formulas (2.27), (2.57), (2.76) respectively. As in the direct formulation, these integrals are numerically computed through the trapezoidal rule. They also need sampled values that are obtained by performing linear interpolation on discrete filtered projections.

2.6.2.b Setup configurations

In the context of computed tomography, particularly in absorption tomography, as explained in section 2.1, a target object is radiated from an external source of some kind of radiation, e.g., X-rays. The radiation interacts with the object and it absorbs some of the radiation. The remaining radiation is collected on detectors placed on the other side of the source with respect to the object.

The geometrical parameters are adapted for each Radon transform. Two parameters should be detailed, namely one describing the 1D detector sampling and the second describing angular variations of radiations. They are chosen for each transform as follows.

Regarding the (p, ϕ)-Radon transform, relative rotations between the object and the system can be made. Thus, every pixel in the input object is radiated equally for all sampled values of the ϕ angle with an angular step ∆ϕ = 0.01 radians. Detectors can be placed perpendicular to radiations and then p is sampled with a step ∆p = 1. 340 samples of p are sufficient to detect radiation crossing the maximal cross section of the object when ϕ = π 4 . The object is centred at the origin and then -120 < x i , y i ≤ 120 while 170 < p ≤ 170. The integration step ∆r is set equal to 1.

In the case of the cartesian Radon transform, the detectors are always parallel to the x-axis (no rotations are made). Radiation directions are described by the tangent of the angle or the line slope and therefore, to compensate the rotationless process, detectors are sampled with -512 < u ≤ 512 and ∆u = 1. The object is also centered at the origin and we have v = tan ϕ sampled with 314 values within [-5, 5]. ∆r is also set to 1.

Finally, as the half-space transform only integrates the function over the plane x > 0, the object is placed obeying 0 < x i ≤ 240 and -120 < y i ≤ 120. As in the cartesian transform, detectors are not rotating. They are parallel to the y-axis and sampled by ζ suche that -512 < ζ ≤ 512 with ∆ζ = 1. Radiation directions are given by ω sampled in ] -π 2 , π 2 [ with ∆ω = 0.01 radians.

Numerical results

Graphical reconstructions of the Shepp-Logan phantom as well as projections are showed in figures 2.8, 2.9 and 2.10. Correlations between input and output objects together with a summary of sampling parameters are presented in table 2.1.

Clearly, better results in terms of correlation are obtained through the (p, ϕ)-Radon transform, then through the cartesian version and lastly through the half-line transform. This is clearly due to the rotation by ϕ of detectors for the first one. Rotations allow to irradiate in all sampled directions the object and then back-project the filtered data. It explains the homogeneous quality in reconstruction for the first transform.

For the second one, some artifacts can be observed on the top and bottom of the reconstruction. It is due to the the lack of data in those sections related to large angle values and thus limited by the length of detectors placed horizontally.

ω ∈] -π 2 , π 2 [ ζ ∈] -512, 512]
The half-space Radon transform suffers from the same limitations than the cartesian one. Moreover, as it only allows radiation towards one side, in this case from right to left to the y-axis, a strong artifact on the right side of the reconstruction is visible, as we have less data on that region related to large angles.

Note that computing time is also more important in the cartesian and half-space versions due to the size of detectors, leading to larger projection arrays to compute, to filter and to back-project.

In this chapter, a simplified two-dimensional preliminary version of the proposed Compton Scattering Tomography configuration is mathematically described. It is based on the Vline Radon transform which models the scattered photon flux density. The related inverse problem is solved by using a filtered back-projection procedure. Numerical simulations on two phantoms representative of flattened heritage objects will be presented in chapter 4.

In the presented 2D framework, note that electrons are considered to be free and at rest. In other words, objects are made of a single material with different electronic densities. In this case, the Doppler broadening and binding effect are not considered in the direct problem and then leading to a bijective relationship between the scattering angle and the related energy-loss. Doppler and binding phenomena will be taken into account in the 3D framework presented in chapters 5 and 6.

The work presented in this chapter has been partly published in [START_REF] Guerrero | Modelling of a New Compton Imaging Modality for an In-depth Characterisation of Flat Heritage Objects[END_REF].

The direct problem

The V-line Radon transform is defined in this first section. Although it has been already exploited in emission configurations in [START_REF] Morvidone | On the V-line radon transform and its imaging applications[END_REF][START_REF] Truong | On new V-line radon transforms in R 2 and their inversion[END_REF], it is adapted and applied here to an imaging configuration based on transmission tomography, that is an external scanning system with a parallel monochromatic source and limited energy resolved detectors.

In the direct problem, the scattered photon flux density is studied, that is the number of photons recorded by a detector per unit time and at a given energy range. Physical factors playing a role in a real scenario, such as the Klein-Nishina cross-section will have to be discussed. 43

The V-line Radon transform

The V-line Radon transform and its inverse formula were first introduced in [START_REF] Morvidone | On the V-line radon transform and its imaging applications[END_REF]. Extensions and applications of the transform in emission imaging are discussed in [START_REF] Truong | On new V-line radon transforms in R 2 and their inversion[END_REF]. More properties of the transform written in cartesian coordinates are developed in [START_REF] Truong | New properties of the V-line radon transform and their imaging applications[END_REF].

Briefly speaking, it is a Radon transform defined on a pair of half-lines of the 2D plane having a V shape. This V-line has a fixed axis direction, parallel to the horizontal axis and its vertex is located on the vertical y-axis. Therefore, it has two geometrical parameters, namely, the opening angle ω, varying in ]0,

V : S(R + × R) → S(]0, π 2 [×R) f → Vf, (3.1)
such that

Vf (ω, ζ) = ∞ 0 f (r cos ω, ζ -r sin ω) + f (r cos ω, ζ + r sin ω) dr r , (3.2)
where dr is the measure on the V-line.

The V-line Radon transform is well adapted in the context of Compton scattering tomography. Indeed, it can be interpreted as the integration of a function, representing an electron density distribution, over all locations of a given object where photons, scattered by the same angle, have the same energy.

Note that the factor 1 r in (3.2) is added because of the photometric law or the difference in solid angle from the scattering site to the detecting site in CST (see below).

Remark 10. This transform can be seen as the sum of two half-space Radon transforms, already discussed in chapter 2, with a modified measure on the half-line and with respective parameters (ω, ζ) and (-ω, ζ).

An equivalent definition for the V-line Radon transform will be useful. It comes from the change of variables t = tan ω, x = r cos ω.

(3.

3)

The V-line transform then becomes a mapping

V : S(R + × R) → S(R + × R), (3.4)
such that

Vf (t, ζ) = ∞ f (x, ζ -xt) + f (x, ζ + xt) dx x . (3.5)

Point Spread Function and kernel function

As pointed out in chapter 2, the Point Spread Function, that is the image of an impulse function, characterizes a transform due to its linearity property. It can easily be obtained after computing the kernel of the integral transform, that is the function k such that

Vf (ω, ζ) = R + ×R f (x, y)k(ω, ζ|x, y) dxdy, (3.6)
Following [START_REF] Morvidone | On the V-line radon transform and its imaging applications[END_REF], we can express the kernel function k with Dirac functions and Heaviside step function H as

k(ω, ζ|x, y) = ∞ 0 δ(ζ + r sin ω -y)H(y -ζ) + δ(ζ -r sin ω -y)H(ζ -y) δ(r cos ω -x) dr r . (3.7)
By making the change of variables u = r sin ω, the kernel k writes as

k(ω, ζ|x, y) = ∞ 0 1 u δ(-ζ + y -u)H(y -ζ) + δ(ζ -y -u)H(-y + ζ) × δ(u cot ω -x) du, (3.8)
and then thanks to the relationship

x tan ω = |ζ -y|, (3.9)
we have an expression of the kernel after performing u-integration

k(ω, ζ|x, y) = 1 ζ -y δ (ζ -y) cot ω -x H(y -ζ) + 1 -ζ + y δ (-ζ + y) cot ω -x H(-y + ζ).
(3.10)

Due to the relationship (3.9), to 0 < ω < π 2 and to the scaling property of the delta function, the kernel is finally expressed as

k(ω, ζ|x, y) = sin ω |ζ -y| δ |ζ -y| cos ω -x sin ω = cos ω x δ |ζ -y| cos ω -x sin ω . (3.11)
The Point Spread Function at P 0 = (x 0 , y 0 ) related to the V-line transform is then derived by writing

Vf 0 (ω, ζ) = R + ×R k(ω, ζ|x, y)δ(x -x 0 )δ(y -y 0 ) dxdy = k(ω, ζ|x 0 , y 0 ), (3.12)
where, from (3.11) we have the expression of the PSF as

Vf 0 (ω, ζ) = cos ω x 0 δ |ζ -y 0 | cos ω -x 0 sin ω . (3.13)
This PSF has an arctangent support in the Radon space (ω, ζ) given by the equation

ω = arctan |ζ -y 0 | x 0 . ( 3

An image formation model in CST

The V-line Radon transform defined above was addressed in a continuous framework. However, the image formation processes has to considerer a realistic imaging configuration with radiation detectors having a finite spectral and spatial resolution. In this section, we give a discrete model related to the 2D Compton scattering tomography. The experimental configuration is a backward scattering configuration and is visible on figure 3.3. Thus, the scattering angle lies in ] π 2 , π]. For sake of simplicity, we call from now on ω = π -ω the supplementary angle of ω, therefore, we will have 0 < ω < π 2 . Also on this figure, it can be seen that the detecting sites D ζ are placed over the Oy-axis, centred at a position (0, ζ) and with a length equal to τ . Incident photons are coming parallel to the 0x-axis towards the positive direction having a constant incident energy E 0 .

3.1.3.a Klein-Nishina differential cross-section

Klein and Nishina in [START_REF] Klein | The scattering of light by free electrons according to dirac's new relativistic dynamics[END_REF] introduced an expression for the differential cross-section for Compton scattering of photons by a free electron at rest. It describes the angular distribution of scattered photons as a function of the scattering angle ω returning the probability of a photon to be scattered into this specific direction. It is equivalently defined as the normalisation of the scattered photon flux density by both the solid angle of the detector and the incident photon flux density [START_REF] Als-Nielsen | Elements of modern X-ray physics[END_REF].

Assuming that the azimuthal angle is uniformly distributed in the interval [0, 2π[ (no polarisation), the Klein-Nishina differential cross-section (or DCS) in a three dimensional medium is given by

σ(ω) = 1 2 r 2 e P (ω), (3.15)
where r e is the classical electron radius and P (ω) is the Klein-Nishina Compton scattering probability given by

P (ω) = E ω E 0 2 E ω E 0 + E 0 E ω -sin 2 ω . (3.16)
The factor 1 2 r2 e comes from the diametrical transversal area of an electron πr 2 e and the uniformity of the azimuthal angle 1 2π

. We also recall that E ω the scattered energy is given by the Compton equation (1.1). The Klein-Nishina differential cross-section in a two dimensional medium is a function giving the probability of a photon to be scattered by ω limited to a bi-dimensional scatter plane. It will not incorporate the uniformity of the azimuthal angle but the factor 1 a two dimensional cone, i.e., a V-line [START_REF] Driol | Modeling and simulation results on high sensitivity scattered gamma-ray emission imaging[END_REF]. Therefore, it is given by

σ 2D (ω) = 1 2
πr 2 e P (ω).

(3.17 -I 0 : the incident photon flux density before the scattering event.

-σ 2D (ω): the 2D Klein-Nishina differential cross-section at an angle ω.

-f (M): the electron density at M.

-Ω(M, D ζ ): the solid angle from M to D ζ .

-dω: the elementary variation of ω.

-dM: the elementary length around M over the V-line.

3.1.3.c Scattered photon flux energy density

It is convenient to express the last angular density as a energy density. This is due to the fact that when dealing with limited energy resolved detectors, we need to integrate the output density into a given energy channel. A density in terms of the scattered energy will be more practical to handle when performing this integration numerically. Th angular density in (3.21) can also be expressed through a change of variables from ω to E ω given by the Compton equation in (1.1) as

dI(E ω , ζ) = τ I 0 cos ω σ 2D (ω)Vf (ω, ζ) dω dE ω dE ω , (3.22)
where

dω dE ω (E ω ) = m e c 2 E ω 1 -cos 2 ω -1 2 = m e c 2 E ω 1 -(1 -m e c 2 ( 1 E ω - 1 E 0 )) 2 -1 2 .
(3.23)

3.1.3.d Limited spectral resolved images

A realistic simulation of a Compton scattering modality requires to consider a discrete energy-resolved detector. Let ∆E be the energy resolution of the detector. Hence, for a fixed detector, we will have a limited number of energy channels recording photons laying on a specific energy interval. Let C i be the energy channel defined for an energy interval

C i = [E i , E i + ∆E].
We are thus interested in I(C i , ζ), the photon counting for an energy channel C i at a detecting site D ζ per unit time, expressed as the integral

I(C i , ζ) = Eω∈C i dI(E ω , ζ). (3.24)
Note that an integration over an energy channel is equivalent to an integration over an angular channel with the correct change of variables and the corresponding boundaries. The boundaries of the energy channel can be equivalently defined in terms of the scattering angle by the Compton equation. Actually, the V-line transform returns a function of the scattering angle from where we get the angular density by (3.21) and then the energy density by (3.22). The boundaries of the channel are set first in terms of energy and then in terms of angle as needed to confine the angular projections.

The direct problem can be then finally modeled by the angular density as

I(C i , ζ) = τ I 0 ω∈C i cos ω σ 2D (ω)Vf (ω, ζ) dω, (3.25)
or equivalently by the energy density as

I(C i , ζ) = τ I 0 Eω∈C i cos ω σ 2D (ω)Vf (ω, ζ) dω dE ω dE ω , (3.26)
where C i is an energy or angular interval according to the case.

The inverse problem

We now need an inverse formula for the V-line Radon transform defined in (3.2) to write down a reconstruction process. A filtered back-projection procedure has been built to this aim and published in [START_REF] Morvidone | On the V-line radon transform and its imaging applications[END_REF]. It is presented in this section.

A first inverse formula

In order to derive a first inverse formula, it is more convenient to work with the (t, ζ) definition of the V-line transform, given in (3.5). The corresponding Fourier Slice Theorem is given below, linking the ζ-Fourier transform of the V-transform with 2D Fourier transforms of the input function f or equivalently with its 2D coupled cosine-Fourier transform.

Theorem 7 (Fourier Slice Theorem). Let f ∈ S(R + × R), q ∈ R and let ḟ be the function

ḟ : (x, y) → 1 x f (x, y), (3.27
)

then F ζ Vf (t, q) = F 2 ḟ (qt, q) + F 2 ḟ (-qt, q), (3.28)
or equivalently

F ζ Vf (t, q) = 2F (cos,y) ḟ (qt, q), (3.29)
where F (cos,y) represents, for the function ḟ , a 2D coupled cosine and Fourier transform related to variables x and y respectively

F (cos,y) ḟ (t, q) = F cos F y (x, q) (t, q) = R 2
ḟ (x, y) cos(2πtx)e -2πiqy dxdy.

(3.30)

Proof. As for the half-space transform, we start by writing the ζ-Fourier transform of Vf and by means of Fubuni's theorem we get

F ζ Vf (t, q) = R×R + ḟ (x, y) R δ(y -ζ -tx)e -2πiqζ dζ dxdy + R×R + ḟ (x, y) R δ(y -ζ + tx)e -2πiqζ dζ dxdy. (3.31)
By definition of the delta function, the 2D Fourier transform appears as

F ζ Vf (t, q) = R×R +
ḟ (x, y)e -2πiq(tx+y) dxdy + R×R + ḟ (x, y)e -2πiq(-tx+y) dxdy = F 2 ḟ (qt, q) + F 2 ḟ (-qt, q).

(3.32)

The proof of (3.29) is derived directly from (3.28). Actually, for a given function g and a given frequency ξ ∈ R, we have the identity Fg(ξ) + Fg(-ξ) = 2F cos (ξ), by developing (3.28) as

F ζ Vf (t, q) = ∞ 0 F y ḟ (x, q) e -2πiqxt + e 2πiqxt dx = 2 ∞ 0 F y ḟ (x, q) cos(2πqxt) dx, (3.33)
where the result stands by definition of the cosine transform.

The integral equation (3.33), that is the cosine transform of the function

x → 1 x f (x, •)
evaluated at qt, can be solved with the related inverse transform and its scaling property. Hence, we obtain:

F y ḟ (x, q) = 1 2 |q| ∞ 0 F ζ Vf (t, q) cos(2πqxt) dt. (3.34)
Finally, f (x, y) can be reconstructed from (3.34) for all points (x, y) ∈ R + × R through the y-inverse Fourier transform giving the following corollary.

Corollary 3. An inverse formula for the V-line Radon transform is given for all

f ∈ S(R + × R) and for all (x, y) ∈ R + × R by f (x, y) = 1 2 x ∞ 0 F -1 y |q|F ζ Vf (t, q) cos(2πqxt) dt. (3.35)

A filtered back-projection inversion

The ramp filter in the Fourier domain appearing in the previous formula suggests to build as before a filtered back-projection reconstruction. This filtered back-projection inversion on V-lines was obtained in [START_REF] Morvidone | On the V-line radon transform and its imaging applications[END_REF]. It combines a back-projection operation and a filtering operation. Lets first define X = S(R + × R) and Y = S(]0, π 2 [×R) with their respective L 2 inner products (•, •) X and (•, •) Y and the relative adjoint transform.

Definition 7. The adjoint transform or the back-projection operator, is then defined as the unique operator

V † : Y → X g → V † g, (3.36) verifying ∀f ∈ X, ∀g ∈ Y, (Vf, g) Y = (f, V † g) X .
(3.37)

This operator performs back-projections for all ω ∈]0, π 2 [, that is, it acts on a function g ∈ Y as

V † g(x, y) = 1 x π/2 0 g(ω, y + x tan ω) + g(ω, y -x tan ω) dω. ( 3

.38)

A back-projection inversion formula for (3.2) can then be written with the adequate filters derived from the inverse formula (3.35). This results is expressed below. Theorem 8. Let the filtered projections by the ramp filter be defined as the function

V * f (t, ζ) = F -1 [|q|F ζ Vf (t, q)] (t, ζ), (3.39)
where t can be replaced by ω if the angular variable is considered. Then, a filtered back-projection inversion for the V-line Radon transform can be performed via the expression

f (x, y) = 1 4 x π/2 0 [V * f (ω, y + x tan ω) + V * f (ω, y -x tan ω)] dω cos 2 ω , (3.40) for all (x, y) ∈ R + × R.
Proof. Starting from the inversion formula (3.35) through Fubini's theorem we have

f (x, y) = 1 2 x ∞ 0 ∞ -∞ |q|F ζ Vf (t, q) cos(2πqxt) e 2πiqy dq dt, (3.41)
and by developing the cosine function

f (x, y) = 1 4 x ∞ 0 dt ∞ -∞ |q|F ζ Vf (t, q)e 2πiq(y+xt) dq + ∞ -∞ |q|F ζ Vf (t, q)e 2πiq(y-xt) dq , (3.42)
where inner integrals are inverse Fourier transforms of the function q → |q|F ζ Vf (•, q), and then

f (x, y) = 1 4 x ∞ 0 dt F -1 (|q|F ζ Vf )(t, y + xt) + F -1 (|q|F ζ Vf )(t, y -xt) . (3.43)
This last expression can be rewritten in terms of filtered projections defined in (3.39) and using the angular variable ω as

f (x, y) = 1 4 x π/2 0 [V * f (ω, y + x tan ω) + V * f (ω, y -x tan ω)] dω cos 2 ω . (3.44)
Note that this formula can also be expressed with the adjoint transform and filtered projections by

f (x, y) = 1 4 x 2 V † 1 cos 2 ω V * f (ω, ζ) (x, y).
(3.45)

Reconstruction from an infinitely energy resolved image

The previous reconstruction by a back-projection procedure requests to perform an angular summation of the filtered projections with respect to the opening angle ω of the V-line. If the output of the system is infinitely resolved in energy, the resulting spectral image is a density in the form dI dE (E, ζ). (3.46) This density in terms of energy needs then to be adapted to a density in terms of angle before the filtering and back-projection process. Moreover, the factors related to the 2D Klein-Nishina differential cross section and the solid angle not included in the V-line transform must also be corrected to this output density.

To convert this energy density to an angular one, we apply a change of variables given by the Compton formula to have an angular density in the form

dI dE (E, ζ) dE dω dω, (3.47)
where

dE dω (ω) = - E 2 0 m e c 2 sin ω 1 + E 0 m e c 2 (1 -cos ω) -2 = - 1 m e c 2 sin ω E 2 ω .
(3.48)

Finally, a correcting factor must be applied to (3.47) by doing

1 τ I 0 cos ω σ 2D (ω) dI dE (E, ζ) dE dω dω. (3.49)
The expression in (3.49) is now an angular density, we can apply to it the filtering and back-projections process detailed previously in this section to obtain an object reconstruction.

Remark 11. In the case of dealing with a realistic limited energy resolved detector, the output dataset {I(C i , ζ)} represent actually a summation of photons arriving in the channel defined by C i , that is, the discrete spectral image. However, our reconstruction procedure is valid from a spectral density dI dE that will be converted to an angular one as explained in the previous section 3.2.3. It remains that we need then to estimate an energy density from the data {I(C i , ζ)} numerically. This will be elaborated on in the next chapter 4 that treats the problem numerically.

Conclusion

A mathematical model for a backward Compton Scattering Tomography configuration has been presented in this chapter. The direct problem, or the image formation process, is modeled by the V-line Radon transform which allows to compute the scattered photon flux density for given imaging parameters, namely the incident energy and the spatial and spectral resolution of detectors. The inverse problem is then solved by a filtered backprojection algorithm related to the V-line Radon transform.

The image formation process should also take into account the chemical composition of the input object. Electrons are considered here to be free and at rest without taking into account the Doppler broadening and binding effects in the direct problem. These effects depend on the different materials that can be present in a heterogeneous object and then the relationship between scattered energy and scattering angle is no longer dominated by the Compton equation (1.1). This additional factors will be treated in the 3D framework, in chapters 5 and 6. Before that, chapter 4 is devoted to some numerical simulations in this 2D context.

Numerical simulations for the 2D framework

This chapter presents numerical simulations for the image formation process as well as for the inverse problem or object reconstruction related to the 2D backward Compton scattering tomography. As explained in the last chapter, simulations will be performed under a back-scattering setup, i.e., data collected with a scattering angle comprised between π 2 and π and following the setup detailed in Chapter 3. Spatial and spectral sampling details, as well as the numerical implementation of the transforms are given.

The discrete version of the V-line Radon transform and the corresponding photon flux density computation is presented for the direct problem, as well as the discrete object reconstruction for two test stratigraphical samples. Different setup parameters are taken into account to evaluate their influence in reconstructions.

The discrete formulation, sampling details

The V-line transform (3.2) integrates over the V-line a regular input function f . For simulation purposes, we deal with discrete objects where we know their value on a finite discrete grid, i.e., an array of M × N values as

{f (x i , y j ), 0 ≤ i < M, 0 ≤ j < N }, (4.1)
where each (x i , y j ) represents the coordinates of a sampled value or a pixel of the object of area ∆x∆y. Here ∆x = ∆y = 2 µm. The V-line transform represents the projections of the input object in the (ω, ζ) domain, these projections are available only from some sampling of both the opening angle and the detector site. If uniform sampling is considered, the projections we dispose belong to the set

{Vf (ω k , ζ l ), 0 ≤ k < K, 0 ≤ l < L}, (4.2
)

57
where K is the number of angular samples and L the number of detecting sites. In presented simulations, sampling rates are chosen to be ∆ζ = 2 µm for the spatial resolution and ∆ω = 0.03 radians regarding the angular resolved image formation simulations. The simulation of spectral images correspondong to finite energy resolution must be obtained from non-uniform angular projections corresponding to the uniform energy sampling. Obtained dataset is then

{I(C i , ζ l ), 0 ≤ i < K ′ , 0 ≤ l < L}, (4.3)
where K ′ is the number of energy channels.

Discrete image formation

To compute the data from an input object, integral (3.2) must be performed numerically to obtain a projection Vf (ω k , ζ l ). The trapezoidal rule is employed and then the V-line transform (3.2) on a point (ω k , ζ l ) is expressed as the sum

Vf (ω k , ζ l ) = ∆r 2 R-2 i=0 (f ikl + f (i+1)kl ), (4.4)
where {r i , 0 ≤ r < R} is a discretization of the V-line onto R points with an uniform step ∆r, and

f ikl = 1 r i f (r i cos ω k , ζ l + r i sin ω k ) + f (r i cos ω k , ζ l -r i sin ω k ) (4.5)
is a projection for each (r i , ω k , ζ l ). The values f (r i cos ω k , ζ l ± r i sin ω k ) requested in sum (4.4) must be obtained from the cartesian sampling of the input object in (4.1). Bilinear interpolation of the discrete object f is used in that purpose as mostly done for the Radon transform due to its good compromise between computational complexity and quality. A nearest neighbour interpolation can be used as well being less computationally expensive while sinc interpolation allows better reconstructions being more expensive related the original Radon transform, this as proven in [START_REF] Toft | The Radon Transform: Theory and Implementation[END_REF].

Photon flux density. By means of the discrete V-line transform in (4.4), we will write down the discrete version of the scattered photon flux density given in (3.26) for detectors resolved both in energy and space. In order to do that, we need to considerer the data we have {Vf (ω k , ζ l )} as in (4.2) recalling that ωk is the supplementary angle of ω k . These data can be uniformly sampled in space on ζ and randomly in angles on ω to produce an angular resolved image. Detector pixels are resolved in energy and not in angle, then for a given detector pixel, an integration over each energy channel with respect to E ω , as explained in section 3.1.3.c, must be done. As the angular spread corresponding to each energy channel has a large dynamic we chose to have an adaptative sampling of ω, setting ∆ω in each energy channel to have the same number of ω samples in each energy channel.

For a single detector site, i.e., at given ζ l and a given channel C i , the discrete version of (3.25) reads by the trapezoidal rule as

I(C i , ζ l ) = τ 2 I 0 ∆ω ω: Eω k ∈C i Vf kl , (4.6)
where

Vf kl = cos ωk σ 2D (ω k )Vf (ω k , ζ l ) + cos ωk+1 σ 2D (ω k+1 )Vf (ω k+1 , ζ l ). (4.7)
The boundaries values of the interval C i may not coincide with a value ω k of the angular sampling due to the non-linearity of the Compton equation and then an interpolation method would be needed to obtain these values. To avoid interpolations, for each energy channel we performed a fixed number of angular projections then each boundary value of C i will coincide with a sampled value ω k . The different energy channels have indeed a different angular width therefore for each energy channel, a different ∆ω sampling is performed.

Discrete filtered back-projection

The input object will be reconstructed over a rectangular grid of the same dimensions, i.e., the reconstruction will be the set {f (x i , y j ), 0 ≤ i < M, 0 ≤ j < N }.

(4.8)

Using the trapezoidal rule to perform the back-projection integral or the adjoint transform (3.38), the back-projection operator, applied to some data {g(ω k , ζ l )} sampled angularly with a constant step, written on its discrete version reads

V † g(x i , y j ) = ∆ω M -1 k=0 g(ω k , y j + x i tan ω k ) + g(ω k , y j -x i tan ω k ) , (4.9)
where ∆ω = π 2K is the angular sampling rate.

This discrete adjoint transform acts on a discrete angular sampling {g(ω k , ζ l )}. To recover the object, we need to apply first the ramp filter q → |q| in the Fourier space of projections Vf , x 2 in the reconstruction space and 1 cos 2 ω in the image space as explained in section 3.2.2.

Linear interpolation from the dataset {g(ω k , ζ l )} is applied to compute the values g(ω k , y j ± x i tan ω k ) requested in the sum (4.9), nearest neighbour interpolation can be used as well. Comparisons of different interpolation methods are studied in [START_REF] Toft | The Radon Transform: Theory and Implementation[END_REF] for the original Radon transform.

A discrete version of reconstruction formula (3.40) is then given, for all (x i , y j ) on the considered rectangular grid by

f (x i , y j ) = 1 4 x i V † 1 cos 2 ω k V * f (ω k , ζ l ) (x i , y j ), (4.10)
through the discrete version of V † in (4.9). V * f (ω k , ζ l ) stands for the discrete filtered projections evaluated by means of (3.39). Filtered of projections in the Fourier domain is done by applying the fast Fourier transform implemented in the R [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF] language.

From the spectral image to angular density. Last reconstruction procedure works acting on an output photon flux angular density, uniformly sampled in the angular variable. However, the simulated data expressed as the spectral image correspond to the summation over finite energy channels of the photon flux density.

To obtain an estimation of the angular density from the spectral image, we consider this spectral image as a piecewise constant function over each energy channel. Then we sample the scattering angle ω from ] π 2 , π[ to have {ω k } with an uniform angular step ∆ω small enough to have more than one sample of ω on each energy channel. From the sampling {ω k } a (non-uniform) energy sampling {E k } is obtained through the Compton equation. Then we assign to each energy sample E k the corresponding value given by the mentioned piecewise constant function. In other words, a discrete energy density is obtained from the spectral image by nearest neighbour interpolation of this image over a energy sampling {E k } that corresponds to a uniform angular sampling.

Finally the uniform sampling of the angular density is estimated from this energy sampling {E k } via the change of variable given the Compton equation

{E dE dω } k . (4.11)
This density can now be filtered and back-projected to obtain the reconstructed object.

Window functions

As pointed out and detailed in section 2.6.1, window functions q → w(q) can be applied to the ramp filter in reconstructions to control high frequencies and then the filtered projections can be written as

V * f (t, ζ) = F -1 [|q|w(q)F ζ Vf ] (t, ζ). (4.12)
From the family of cosine windows having the form w(q) = cos n (πq), (4.13) where the case n = 2 corresponds to the well characterised Hann Window function, we tested various values of n during our simulation work and found out that n = 2 was the one maximising the correlation between the object and its reconstruction given the numerical noise produced by the simulation. Still different values of n gives satisfactory reconstructions as well.

y ζ m -ζ m x l τ x 0 O detectors object incident beam (E 0 ) E ω ζ 0 -ζ 0 D ζ M ω dΩ

Imaging configuration details

Simulations will considerer the following physical parameters for the imaging system: Incident energy of the X-ray source and its cross section length, spectral and spatial resolution of the detector, length of the detector, position relative between sample, source and detector with the corresponding maxima and minima detection angle limitations.

These parameters are illustrated in Figure 4.1. As mentioned there, vertical translations are needed to overcome the limitations of the beam size l which induces systematic missing data around ω = 0. We are considering in these simulations parameters detailed in table 4.1 for an ideal detector resolved angularly and in table 4.3 for a limited energy resolved detector.

In present simulations, attenuation and multiple scattering are neglected. The Doppler broadening and binding effects will be considered in the 3D case in chapters 5 and 6.

Test objects

Two test phantoms are going to be considered in the simulations for the 2D case. Namely, a numerical stratigraphic phantom created by means of random layers and interior grains and a real sample composed mostly by copper and related corrosion products reconstructed by classical absorption X-ray tomography and filtered back projections. Each value of a pixel represents the electron density of the object.

The first one presents three layers of about 150 µm of thickness each and 128 µm of height. They have electron densities of 0.9, 1.1 and 1.0 respectively. Grains of random diameter and position were inserted inside the layers, having densities of 1.3 ± .1, 6.0 ± 0.5 and 2.0±0.3 respectively. We consider 1 pixel to represent 2×2 µm 2 and then the phantom is a matrix of 284 × 64 pixels.

The second one is a real sample composed by different copper corrosion products. Reconstructed from a real imaging experience and segmented for sake of simplicity into 6 classes of electron densities by taking the mean values of each class. It is a matrix of 380 × 64 pixels where each pixel represents again 2 × 2 µm 2 .

The first phantom presents a more defined morphology where each element has a regular boundary and big discontinuities between them. Reconstructions on it will expose the capacity of the model to identify these discontinuities. The second phantom has irregular boundaries for the different layers and elements but the contrasting between different elements are less important. Reconstructions on it will characterise the capacity to recover irregular boundaries. Both phantoms are represented in Figure 4.2. The first performed simulations regard a simple application of the V-line Radon transform and its inversion without considering physical constraints. The phantoms considered are the simulated analytical stratigraphical and the real heritage sample presented in Fig- Sampling parameters, namely the number of projections and the number of detectors are given as follows. Concerning the angular sampling rate, 100 and 200 projections uniformly distributed in angle are considered. Concerning detectors spacial sampling, 2 10 and 2 11 detector sites are considered, in every case they are separated by 1 pixel, i.e, ∆ζ = 2µm.

Simulation results

Uniform angular sampling

First, projections are obtained from a direct application of equation (4.2) to the whole object and applying the physical factor in the image formation cos(ω)σ 2D (ω), they are showed in Figure 4.3 for 200 uniform sampled projections and 2 11 detecting positions. Due mostly to the Klein-Nishina DCS, most of the scattered energy is detected at low angles that is used to recover the global depth structure of the object at a correct scale. The output angular image exposed at logarithmic scale shows the shape of the image for large angular values, data coming for these angles records the contrasting and morphological details of the object.

These projections are then filtered by the ramp filter and the Hann window function in the form of (4.13) when n = 2 to obtain the discrete filtered projections written in (4.12). The choice n = 1 could be used as well as this image formation is done without introducing form uniformly sampled on the opening angle given in (4.9) for every pixel in the reconstruction grid.

Graphical results are show in Figures 4.5 and 4.6 for both phantoms and for the different sampling parameters summarised in table 4.2. Correlations coefficients between the two input objects and their reconstruction are presented in the mentioned table. The correlation coefficient measure the degree of linear dependence between two vectors in R n , we applied it here to obtain a scale invariant measure of comparison between the phantom and its reconstruction. Geometrically, the correlation coefficients is also the cosinus of the angle between the two centred vectors compared and is then bounded in the interval [-1, 1]. 1 being the correlation of two linearly dependent vectors. The arccos of the correlation coefficient is also presented in table 4.2 that stands for the angle between those two vectors; the original and the reconstructed objects, providing a finer contrast to compare correlation values close to 1.

Conversely, to be able to characterise the variation in scale between the original phantom and its reconstructions, we are using the ratio of the variances, that is the measure Var(reconstruction) Var(phantom) ,

where a value of 1 means that the reconstruction is on the same absolute scale as the phantom.

More projections implies less (almost) horizontal artifacts in reconstruction due to a better interpolation in back-projections for small angles. Larger detectors in absolute length means that we have more data for big angles and so less (almost) verticals artifacts presented in reconstructions generated from the lack of back-projections coming from those missing angles. Simulations taking into account a spectral resolved detector are going to be presented in the next section.

Uniform energy sampling

We present in this section the image formation and object reconstruction for a 2D backward Compton scattering tomography considering a more realistic scenario, i.e., finite energy resolved detectors, a beam of a limited cross-section length and incomplete data due to the hole in the detector. Setup parameters are summarized in table 4.3. The beam cross-section width is set to be 32 µm and the phantoms have a width of 128 µm therefore we performed 8 vertical translations of 16 µm to allow the imaging of the full object. As mentioned in the angular case, regions of 16µm of width are partially overlapping between two consecutive translations and only central part of each is conserved to avoid edge artifacts in reconstructions, these artifacts are due to the different scale between reconstructed bordering and central pixels. These boundaries artifacts are caused both by interpolation problems of bordering pixels in the direct problem creating aliasing effects and also for the rectangular support of the region of interest.

4.4.2.a Spectral image formations

The image formation process must be performed for each slice. We selected first a heterogeneous slide of each phantom showed both in Figure 4.7. We recall that the direct problem consist in obtaining the scattered photon flux density per both detecting site and energy channel, i.e., we aim to compute I(C i , D ζ ) defined in (4.6) related to the selected slices to then estimate the angular density dI(C i , D ζ ).

The three stages of the problem are performed. Namely, from the sample slice we compute the projections or the angular density, i.e., the V-line projections times the Klein-Nishina and solid angle factors. This density is presented in the left columns of Figure 4.8. Then, summations over each energy channel of projections are performed to obtain the finite resolved spectral image in the second column. Finally, the angular density performed initially is estimated from the spectral image by nearest neighbour interpolations in the third column. This latter step will be clarified in the next section.

The stratigraphical phantom slice presents two grains of about 6 times the electronic density of the rest of the phantom. This two grains can be identified in Figures 4.8 having a support in the form of the PSF of the V-transforms studied in section 3.1.2.

As the unit length considered is 2 µm then the spatial resolution in reconstructions is 2 × 2 µm 2 per pixel. In this backscattering setup, a single experience simulate a X-ray beam crossing a hole of 20 µm in a 1D detector, hence we have ζ 0 = 10 µm. The 1D detector is located vertically over the line x = 0 and the phantom is placed at different distances x 0 to the detector. The detector has a vertical length of L = 2 10 pixels including the hole. Thereby, we have L -24 detecting sites and a resolution of 2 µm. (24 sites not considered due to the hole.)

Regarding the trapezoidal rule in calculating integrals. In the direct problem the integration spatial step over V-lines is set to dr = 2 µm and dω = 0.005 radians the angular step to compute the limited spectral image from angular sampled projections.

4.4.2.b Reconstructions

Phantoms are reconstructed by means of the filtered back-projection procedure detailed in section (4.1) for each slice. First, as explained in section 4.1.2, from the spectral image, the angular density is estimated by nearest neighbours of each energy channel. The resulting angular density is presented in Figure 4.8 where we can appreciate a good estimation for all channels but for the first one, the channel covering the larger angular interval.

Figure 4.9 shows the piecewise constant function representing the spectral image that is constant on each channel for two detecting sites, namely ζ ∈ {-12, 88}. Different spectral resolutions are represented where the case ∆E = 1 eV suggests the true energy density. We can appreciate the errors in interpolation mostly on channels where the derivative of the density is seen to be larger than other channels. The first channel is a particular case because the derivative is bigger on the centric zone of the detector where most of the signal is detected. There, errors in interpolations implies huge blurring artifacts in reconstructions.

This issue could be addressed by performing a more regular interpolation of the spectral image than nearest neighbours to estimate the energy density, a C 1 interpolation would avoid the seen step effect between channels and also a better approximation of the derivative of the density inside a channel.

We present the reconstruction procedure and results detailed in section 3.2 for the slice of the copper sample. The data to filter and back-project is the estimated angular density i.e., data in the right column of density. The filtering process can be performed, via (3.39) or also via (4.12) if window apodization functions are used in the Fourier domain, applied to the corrected angular density.

Figure 4.10 illustrates the various steps of the filtering process for the stratigraphical phantom and a 50 eV resolved detector. The 1D Fourier transform of the angular density per each angle are performed. The large errors introduced by the interpolation can be appreciated mostly in the first channel, these are reduced by assigning to these channel the values of the second channel in the next figure. Then we apply the ramp filter where we can appreciate numerical noise due to the divergence of the filter. In the next image, a cosine apodization window in the form of (4.13) with n = 2 is applied to the ramp filter to control the numerical noise, the value of n is chosen due to the fact that we have not introduced Poisson noise and so there is not the need to over-smooth the reconstructions with higher values of n, still n = 1 is not enough due to interpolations numerical noise.

Finally, the 1D inverse Fourier transform is performed to obtain filtered projections. We can appreciate the acting of filters on the filtered projections comparing them with the original projections. Contrasting features as those produced by the frontiers of different layers or elements are better appreciated. The consequences of systematically missing data due to the presence of the hole can be seen in figures, even if there are only few missing pixels. The intensity of the effect is due to fact that a large part of the signal would have been detected on these missing sites. Reconstruction will suffer mostly a error in scale due to this missing data.

From filtered projections as those showed in Figure 4.10f, back-projections and then a reconstruction of the slice is obtained through (3.40).

The incomplete data due to the presence of a hole in the detector affects the reconstruction at low opening angle projections. Placing the object further to the detector allows us to reduce the missing data related to low angles but we loose data corresponding to large angles. That is, for a given half-detector determined by ζ 0 and ζ m in Figure 4.1, an input pixel in the object located at the middle of the incident beam of coordinates (x, 0) generates diffracted radiation at respective minimum and maximum opening angles ω m (x), ω M (x) given by

ω m (x) = arctan( ζ 0 x ) ω M (x) = arctan( ζ m x ), (4.14)
for a given x ≤ x 0 . The amount of diffracted data related to this point will be given then by the integral

g(x) = ω M (x) ωm(x) σ 2D (π -ω) cos(ω) dω. (4.15)
An optimal value of x can be easily found that maximise the (concave) function g. It derivative is

g ′ (x) = σ 2D (π -ω M (x)) cos(ω M (x))ω ′ M (x) -σ 2D (π -ω m (x)) cos(ω m (x))ω ′ m (x)
. (4.16) Then the x position of a pixel generating more data is the x verifying g ′ (x) = 0. A simple look for the root of this equation numerically gives x = 142 pixels for our case that depends on lengths of both the detector and the hole.

Therefore, the choice of the distance x 0 to place the object from the detector is determined by the region of interest that we want to recover and not much bigger than 142 pixels due to missing data from large angles. For example, if we want to recover pixels at 40 µm (20 pixels) inside the object, x 0 = 100 pixels (200 µm) should be a good choice. The functions g as a function of x is showed in Figure 4.11 where the maximum is at 142 pixels.

4.4.2.c Results and conclusions

Five energy resolutions of detectors are considered in simulations, namely ∆E ∈ {1, 5, 25, 50, 100} eV, as well as an ideal detector not resolved in energy but in angle. Different distances x 0 from the phantom to the detector are also considered.

First, an ideal detector resolved angularly is considered. Simple reconstruction is performed with 16 horizontal translations of 8 µm having no overlaps between them and allowing the imaging of the full sample. Edge artifacts are present due to different scale rates between edge pixels and centric pixels. The bordering horizontal pixel lines suffer for a huge decay in scale, this phenomenon is also present in the classical Radon transform.

These artifacts are eliminated by performing vertical translations of a smaller length than the complete ones, in such way that the top and bottom pixel lines can be dropped, i.e., we perform 32 translations of 4 µm in our case. This solution shown in Figure 4.12 will be employed in the rest of the section for all reconstructions. Figure 4.13 presents reconstructions of the copper sample with different distances x 0 to the detector, with 4 µmtranslations and removing bordering pixels for each slice. The effect of the hole is appreciated in the scale variations with respect to Figure 4.12 that has almost the correct scale variation. The horizontal artifacts presented in these reconstructions are due to the missing data due to the hole. As it has been said, pixels closer to the detector suffer more from this missing data, as it is seen on Figure 4.13a where x 0 = 1 pixel, then it has the full data from positions x > 142 where the artifacts are no longer present as discussed previously. Figure 4.15 presents the profile for 1 pixel slice of reconstructions, it is clear there the difference in scales when the simulation do not have a hole in the detector. The correlations coefficients plotted in Figures 4. [START_REF] Guerrero | Three-dimensional imaging of flat natural and cultural heritage objects by a compton scattering modality[END_REF] choice of x 0 will depend then in the target feature of the object studied.

A finite spectral resolved detector was considered afterwards. Figure 4.17 shows full reconstructions of the copper phantom varying the distance x 0 to the detector that is resolved in energy with a uniform resolution of 50 eV. One can see a correct identifications of x-positions of different electronic density. The regions of the sample that are further to the detector presents less contrasting features between different materials. This is due to the lack of data related to larger angles, this issue can be addressed by expanding the This effect is indeed reduced by placing the object further at the price of losing information corresponding to large angles. Figure 4.18 shows reconstruction results of the same sample placed at x 0 = 50 pixels, that is, as in Figure 4.17d with different spectral resolutions, namely 1, 5, 25 and 100 eV where we can appreciate the blurring in reconstructions affected by a larger spectral resolution. Actually, with a 5 eV resolution we dispose of more than 700 projections where the largest angular channel spread corresponds to 0.035 radians. This angular channel corresponds to the energy channel covering from E π to E π + ∆E, that is, the channel related to low angular values. With a resolution of 50 eV we have 74 projections and the largest angular channel is 0.171 radians wide while for 100 eV, we have 37 projections and a maximum angular channel of 0.242 radians width. Notice that reconstruction in 4.13 from angular resolved data were obtained with an angular step of 0.005 radians.

An ideal detector resolved in angle is obviously not realistic and gives better results in terms of correlation coefficient as the angular sampling is uniform and then no numerical artifacts are generated due to interpolations from spectral images to energy densities. Without considering the Doppler broadening effect, reconstructions depend exclusively on the energy resolution of the detector and the quality of interpolations. Therefore, an energy resolution broader than 100 eV (corresponding to 37 non uniform angular projections) will not provide sufficient observations. As mentioned, nearest neighbour interpolation of the spectral image could be improved by applying a C 1 interpolation to avoid discontinuity effects in densities and better estimation of their derivatives.

The farther the sample is placed, the less data related to big angles will be obtained. However, if the sample is placed closer to the detector, the effect of the hole will be bigger as more data will be lost at low angular values. This is seen in the variance ratio for 1 and 5 eV that gets better by placing the object further to the detector. When dealing with larger spectral resolutions, the interpolations errors imply a huge decay in scale and then in the variance ratio. Translations of half the size of the beam cross-sections allow to eliminate edge pixels that are not reconstructed a the correct scale generating huge edge artifacts. However, it will require twice the time to complete the experiment than the case when no intersections between translations are performed.

Correlations coefficients and variance ratios are presented in table 4.4. Of course higher spectral resolution implies better coefficients. We can appreciate the huge lost in scale by looking the resulting variance ratios due to all the data that is lost in summation per energy channel together with the data not measure in the hole region. The correlation coefficients suggest that the different materials of the sample are well identified in depth but with an important lost in their morphology in deeper regions of the sample due also to the limited detector length.

The results presented in this chapter according to the different setup configurations, sampling parameters and the numerical implementations will help us to design numerically the three dimensional device that is presented in the next chapters. The different artifacts in reconstructions are better understood now as well as the different errors caused by numerical approximations of integrals and interpolations. As the 3D scheme is more computationally expensive, this understanding will allow us to set the different parameters in simulations. A GPU parallelization of both the direct and inverse problems is also developed reducing the time complexity. The different algorithms are exposed in Appendix A.

Backward Compton Scattering Tomography: 3D framework

In this chapter, a three-dimensional Compton scattering mathematical model is developed both for the direct and inverse problems. As in the 2D framework studied in chapters 3 and 4, the system to be modeled receives an input function defined now on some domain of R 3 representing the electronic density of an object and returns an output function that stands for the measured spectral image. First, the direct problem is studied in a complete physical context, that is in a continuous and also a discrete framework. It is expressed with the help of the so-called conical Radon transform, that is a transform where the output function is an integral over circular cones surfaces with a fixed axis direction. Then, the inverse problem is solved with a filtered back-projection procedure of the output spectral image leading to the estimation of the original 3D distribution of the electronic density of the object. This chapter is adapted from the reference [START_REF] Guerrero | Three-dimensional imaging of flat natural and cultural heritage objects by a compton scattering modality[END_REF], partly published in Journal of Electronic Imaging in 2017.

The direct problem

Given an input function f , the imaging system returns an output function I, the energy resolved spectral image, defined on some bounded domain of R 2 × [E m , E 0 [ where E m is the lowest energy of a photon after a Compton interaction (E m is indeed near ω = π in the Compton equation) and E 0 the incident energy of photons.

In the case where photons are scattered by free electrons at rest, the conservation of energy and momentum implies a diffeomorphism between the scattered resulting energy of photons and the scattering angle given by the Compton equation (1.1). The conical Radon transform performs conical projections of the input object where all photons scattered at a given energy lay over a circular conical surface.

Actually, as electrons are not at rest and move with a certain momentum distribution, the Doppler broadening effect must be taken into account. Moreover, as the electrons are not free, their ionization energy determines whether or not an interaction occurs leading to the so-called binding effect [START_REF] Salvat | PENELOPE, a code system for Monte Carlo simulation of electron and photon transport[END_REF]. In a realistic scenario where electrons belong to a given material, the Doppler broadening and binding effects must be taken into account in the image formation process. Then, the input function will need to be defined on a larger domain that must contain information about the material in addition to the electronic density.

The first section of this chapter deals with the mathematical definition of the conical Radon transform, then a model of the spectral image formation by this Compton scattering modality is given. In the second part of this chapter, the associated inverse transform that can approximate a reconstruction from a spectral image resolved in energy is studied.

The conical Radon transform

The conical Radon transform was first defined in [START_REF] Cree | Towards direct reconstruction from a gamma camera based on compton scattering[END_REF] without an inverse formula but the imaging application based on Compton scattering of gamma rays was proposed. In [START_REF] Nguyen | Radon transforms on a class of cones with fixed axis direction[END_REF], a first inversion formula using circular component analysis, with applications to emission imaging was introduced. Generalisations of this kind of transforms to higher dimensions spaces with related inverse formulae were presented in [START_REF] Haltmeier | Exact reconstruction formulas for a radon transform over cones[END_REF]. A back-projection inversion algorithm for a conical Radon transform in R 3 was developed in [6].

The transform output function integrates a function over circular cone surfaces. It is the 3D analogue of the V-line Radon transform due to the the fact that the axis direction is fixed and the vertices lay on a 2D fixed plan. Here, a cone is uniquely determined by the vertex position and the opening angle, while its parametrization parameters are the azimuthal angle and the slant heigh or the distance of a point on the surface to the vertex. Such cone has been already depicted in Figure 1.4 in the introduction.

Assume that the vertices of the conical surfaces lay on the x0y-plane and that the input function is defined for z > 0. A (ψ, r)-parametrization of the circular cone surface is considered for a fixed opening angle, where ψ is the azimuthal angle and r the euclidean distance to the vertex. Such parametrization for a cone with a vertex located at D = (ζ, ξ, 0) and an opening angle ω reads C ω,D (ψ, r) = D + (r sin ω cos ψ, r sin ω sin ψ, r cos ω),

(5.1)

where ψ ∈ [0, 2π[, r ∈ R + , ω ∈]0, π 2
[, and (ζ, ξ) ∈ R 2 . The conical surface may also be characterized by a cartesian equation valid for all (x, y, z

) ∈ R 2 × R + by z 2 tan 2 ω = (x -ζ) 2 + (y -ξ) 2 .
(5.

2)

The conical Radon transform applied on a electron density function f returns an output function, called Cf on R 2 ×]0, π 2 [ that performs the integral of f over a cone surface parametrized as in (5.1).

Definition 8. The conical Radon transform is a mapping (5.4)

C : S(R 2 × R + ) → S(R 2 ×]0, π 2 [) f → Cf, (5.3
The measure on the conical surface r sin ω dψ dr times the photometric spreading of radiation factor 1 r 2 must be included in the last integral. Note that, for simplicity reasons, the factor sin ω does not explicitly appear and will be included later in the image formation model.

It may be useful to express (5.4) as an integral with respect to the z-axis with the change of variable z = r cos ω and also to introduce a new variable t such that t = tan ω.

If notation Cf is conserved, (5.4) may be expressed as

Cf (ζ, ξ, t) = ∞ 0 2π 0 f (ζ + tz cos ψ, ξ + tz sin ψ, z) dψ dz z , (5.5)
where t ∈ R + .

5.1.1.a Kernel and Point Spread Function

The point spread function (PSF) is the image of the input impulse function

f 0 (x, y, z) = δ(x -x 0 , y -y 0 , z -z 0 ), for (x 0 , y 0 , z 0 ) ∈ R 2 × R + . (5.6)
Its conical Radon transform can be derived as [6] Cf

0 (ζ, ξ, ω) = 1 z 2 0 tan ω δ(z 0 tan ω -(x 0 -ζ) 2 -(y 0 -ξ) 2 ). (5.7)
The support of the PSF in the (ζ, ξ, ω) domain is given by concentric circles of radius z 2 0 tan 2 ω centered at (x 0 , y 0 ) for all ω. This arctangent type support is represented in Figure 5.1 for an impulse function at (0, 0, 1). Note also that the kernel of the conical Radon transform is then given by the expression

k(ζ, ξ, ω | x, y, z) = 1 z 2 tan ω δ(z 2 tan 2 ω -(x -ζ) 2 -(y -ξ) 2 ). ( 5

Image formation model

The aim here is to estimate the output spectral image returned by this CST system from a given input object. First, a rough model is going to be derived when the input object consists of a function representing the electronic density distribution where electrons are free and at rest. Then, the function will also include different material characteristics to model the Doppler broadening and binding effects.

Define dI(E, D|M ) as the recorded scattered photon flux density from a scattering position M (number of photons of energy laying in dE recorded per unit time at a detecting region D coming from M ). The recording point D can be characterized by the cartesian coordinates of its middle point (ζ, ξ, 0) (see figure 1.4 already depicted in the introduction for a description of the geometrical configuration). The value dI(E, D|M ) then incorporates the following parameters:

-I 0 : the incident photon flux density just before the scattering event at M.

-σ(ω, E): the Compton scattering double differential cross-section (DDCS) at a scattering angle ω and with a scattered energy E.

-f (M): the electron density at M.

-dE the elementary variation of E.

-Ω(M, D): the solid angle from M to D.

-dM: the volume element around M.

The solid angle Ω(M, D) can be seen from (5.10)

The Compton scattering differential cross-section (DCS)

The Compton differential cross-section characterises each interaction between photon and matter as a function of the photon state variables, that is, as a function of the scattering angle and scattered energy. When dealing with free-electrons at rest the differential crosssection (DCS) is a function of only the scattering angle leading to a diffeomorphism between the angle and energy. Contrary to the 2D framework where the free-electrons hypothesis was the only one considered, we will here also consider the hypothesis that the Compton differential crosssection is a double differential cross section (DDCS) defined as a function of both the energy angle and also the scattering energy of the particle. This more complete formulation makes it possible to include both the Doppler broadening and the electron binding energy into the photon inelastic scattering model.

The general model will thus not perform integrations over a conical surface but over the whole input domain, that is

dI(E, D) = R 2 ×R + dI(E, D|M) = τ I 0 dE R 2 ×R + 1 r 2 M cos ωM σ(ω M , E) f (M) dM.
(5.11)

5.1.3.a The Compton scattering DCS for free electrons

As in the 2D case, the Klein-Nishina differential cross-section σ KN could be considered in the image formation process when electrons are free. It is a function of only the scattering angle ω giving the probability of a photon to be scattered in a given direction ω when the azimuthal angle is uniformly distributed in the interval [0, 2π[. This differential crosssection assumes an exact relationship between scattered energy and angle given by the Compton equation (1.1). It can be seen as a Dirac probabilistic law that reads σ(ω, E) = σ KN (ω)δ(E -E ω ), (5.12) where E ω is the energy verifying the Compton equation and σ KN the Klein-Nishina function already defined in chapter 3 limited to a 2D plane and recalled here for the general 3D case by

σ KN (ω) = 1 2 r 2 e E ω E 0 2 E ω E 0 + E 0 E ω -sin 2 ω .
(5.13)

The Klein-Nishina DCS indicates that the probability of a photon to be Compton scattered is higher for scattering angles close to 0 or π than for angles close to π 2 . The support of this DCS in the ω -E plane is given by the line verifying the Compton equation, that is a non-linear decreasing function of the scattering angle. Both curves are presented in Figure 5.2. In this particular and simpler case, as the DCS depends only on the scattering angle, the scattered flux energy density dI(E, D) in (5.11) is simplified onto an integral over a conical surface thanks to (5.12) by doing (5.14) where C ω,D is the cone with an opening angle ω and a vertex located at D.

dI(E ω , D) = τ I 0 dE ω M: E=Eω 1 r 2 M cos ωM σ KN (ω M ) f (M) dM = τ I 0 cos ω σ KN (ω)dE ω C ω,D 1 r 2 M f (M) dM,
In terms of the conical Radon transform of f previously defined the energy density can also be expressed as

dI(E ω , D) = τ I 0 sin ω cos ω σ KN (ω)Cf (ζ, ξ, ω) dE ω .
(5.15)

Notice that the integration measure over the cone surface dM = r sin ω dψ dr explains the presence of the factor sin ω in this expression. Notice also that an angular density can be derived by writing

dI(ω, D) = τ I 0 sin ω cos ω σ KN (ω)Cf (ζ, ξ, ω)
dE ω dω dω.

(5.16)

5.1.3.b The general Compton scattering DDCS in heterogeneous objects

As already said, electrons are not at rest but move with a certain momentum distribution p. This phenomenon gives rise in the spectral image to the Doppler broadening effect. Furthermore, electrons are not free and then interactions are allowed only if the energy transfer is larger than the ionisation energy U i of the atomic electron shell i where the interaction has occurred, i.e., if E 0 -E > U i , where E is the energy of the scattered photon. This is known as the binding effect. Hence, some electrons will not participate in the inelastic scattering, depending on E 0 and on their binding energy. In this general case, the Compton Equation (1.1) and the Klein-Nishina differential cross-section given in (5.13) are not longer valid. However, they still remain a good approximation for incident energies up to the order of few MeVs [START_REF] Salvat | PENELOPE, a code system for Monte Carlo simulation of electron and photon transport[END_REF].

The corresponding Compton DDCS takes into account the material influence generating the Doppler broadening and binding effects. It is a function of both the scattering angle and the energy and it is sampled in Monte Carlo simulations, as done in [START_REF] Brusa | Fast sampling algorithm for the simulation of photon compton scattering[END_REF][START_REF] Salvat | PENELOPE, a code system for Monte Carlo simulation of electron and photon transport[END_REF] by

σ(ω, E) = σ KN (ω)F (p z ) i Z i J i (p z )H(E 0 -E -U i ) dp z dE .
(5.17)

In this expression, p z is the projection of the initial momentum p of the electron on the direction of the scattering vector and H the Heaviside step function. J i (p z ) is the Compton profile of the i-th atomic shell, Z i is the number of electrons in the i-th shell and the function F is a rejection sampling function [START_REF] Brusa | Fast sampling algorithm for the simulation of photon compton scattering[END_REF]. In the Penelope code [START_REF] Salvat | PENELOPE, a code system for Monte Carlo simulation of electron and photon transport[END_REF], the functions F and J i are approximated by an analytical expression. J i also needs the value J i (0) per material, that is sampled and tabulated.

For a given ω, this DDCS is bell-shaped and symmetrical about the Compton line E = E w . The spreading of the distribution varies according to the incident energy, the angle and of course the material. This DDCS is ploted in Figure 5.3 as a function of E for two scattering angles and two elements, hydrogen and aluminium. The spreading of this density increases with the scattering angle and with the atomic number. The Compton line is also represented showing that the Klein-Nishina DCS is an approximation of the DDCS by a Dirac law located at its expectation value.

Limited spectral resolved images

As in the 2D case, let C i = [E i , E i+1 [ be a given energy channel related to the detector spectral resolution. The Compton Scattering Tomography system will collect the spectral image I(C i , D), or the summation of the scattered photon flux density over the energy channel C i at a detecting site D, that is

I(C i , D) = E∈C i dI(E, D), (5.18)
where dI(E, D) is the energy density defined in (5.11).

In the case of free-electrons at rest when the DDCS is given by the Klein-Nishina DCS and the Dirac function over the Compton line, the last integral is indeed an integral over such Delta function that is given by the Heaviside step function H, that is

I(C i , D) = τ I 0 R 2 ×R + 1 r 2 M cos ωM σ KN (ω M ) E∈C i δ(E -E ω ) dE f (M) dM = τ I 0 R 2 ×R + 1 r 2 M cos ωM σ KN (ω M ) H(E -E i+1 ) -H(E -E i ) f (M) dM. (5.19)
The last model of the direct problem differs from the one in the 2D framework where the equivalent here is a summation of the angular density on every related angular channel C i as

I(C i , D) = τ I 0 ω∈C i sin ω cos ωσ KN (ω)Cf (ζ, ξ, ω)dω = C i ∞ 0 2π 0 f (ζ + r sin ω cos ψ, ξ + r sin ω sin ψ, r cos ω)σ KN (ω) cos ω 1 r 3 r 2 sin ωdψdrdω, (5.20)
where the elementary volume element dM = r 2 sin ωdψdrdω = dxdydz appears and then we can perform this integral on a cartesian coordinate system. We need then the characteristic function over the volume element laying between two conical surfaces δ C i ,D (x, y, z) given by

δ C i ,D (x, y, z) = 1 for cos ω = z √ (x-ζ) 2 +(y-ξ) 2 +z 2 ∈ cos C i , 0 otherwise.
(5.21)

Therefore, we can express the last integral as

I(C i , D) = R 2 ×R + f (x, y, z) δ C i ,D (x, y, z) σ KN (ω M ) cos ωM r 3 M dM. (5.22)
The numerical implementation of this transform is detailed in the next chapter. It is based either on the first expression (5.4) or on the last cartesian integral (5.22) more suited to a parallel implementation due to its cartesian nature. This will be detailed in the coming chapter 6 and can also be visible in the annex A.

In the more general case when electrons are not free, the full DDCS in (5.17) must be applied in the last integrals per each material, angle and energy. This may not be possible numerically, then we can approximate this DDCS by a probability law f ω of the scattered energy per each angle as σ(ω, E) ≈ f ω (E), (5.23) having its expectation value at the Compton energy E ω as the DDCS is a spreading of the Compton line. See Figure 5.3. Notice that in this case we can not perform the integral over an angular channel instead of the energy channel C i because they are not equivalent due to the diffeomorphism angleenergy is not longer available. We must then perform integrals with boundaries as in (5.18) and not as in (5.20).

If F is the cumulative distribution function of this law, the integral (5. [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF]) is rewritten as

I(C i , D) = τ I 0 R 2 ×R + 1 r 2 M cos ωM F (E -E i+1 ) -F (E -E i ) f (M) dM.
(5.24)

In the next chapter where the problem is treated numerically, we will explain the approximation of the DDCS by a Lorentzian law.

The inverse problem

As in the 2D case, we present here the related inverse formula to the conical Radon transform (5.4) and then a filtered back-projection reconstruction for the output spectral image. First, the Hankel transform is introduced as a preliminary step to derive the inverse formula.

Auxiliary functions

The inverse expression of the Hankel transform, defined through the Bessel function, is used to solve the conical transform in (5.4) due to the relationship between the 2D (ζ, ξ)-Fourier transform of the input function f and the Hankel transform of the conical function Cf . The definitions of the Bessel and Hankel functions are recalled below. Definition 9. The first kind Bessel function of order 0, denoted J 0 , is defined for x ∈ R by the identity

J 0 (x) = 1 2π
2π 0 e -ix cos(α-φ) dα, (5.25) for all φ ∈ [0, 2π]. The Hankel transform, or the Fourier-Bessel transform of order 0, denoted h 0 , is then defined for a function g as a weighted integral of Bessel functions by h 0 g(k) = R+ g(r)J 0 (kr)r dr, (5.26) for all k ∈ R + . Its related inverse transform is

g(r) = R+ h 0 g(k)J 0 (kr)k dk.
(5.27)

A first inverse formula

A first analytical inverse formula, announced in [6], by means of the Hankel transform for the conical Radon transform with a fixed axis direction comes in the following theorem.

Theorem 9. Let f in S(R 2 × R + ) and C the conical Radon transform defined in (5.4).

The related inverse formula can be expressed as

f (x, y, z) = 2πz 2 R 2 (u 2 + v 2 )e 2πi(ux+vy) π/2 0 sin ω cos 3 ω Cf (u, v, ω)J 0 (2πz tan ω u 2 + v 2 ) dω dudv.
(5.28)

Proof. In order to derive the inverse transform, we start writing the bidimensional (ζ, ξ)-Fourier transform of the conical Radon transform expressed with t = tan ω in (5. (5.31)

Last inner integral is actually the (x, y)-Fourier transform of f at (u, v) noted f (u, v, •). Then, the last expression writes as:

Cf (u, v, t) = ∞ 0 f (u, v, z) 2π 0 e 2πitz(u cos ψ+v sin ψ) dψ dz z , (5.32)
We can then identify the first kind Bessel function of order 0 by switching to polar coordinates: u = q cos β and v = q sin β:

Cf p (q, β, t) = ∞ 0 fp (q, β, z) 2π 0 e 2πiqtz cos(β-ψ) dψ dz z = 2π ∞ 0 fp (q, β, z)J 0 (2πqtz) dz z .
(5.33)

Note the use of notation Cf p and fp to point out the use of polar coordinates. Last integral turns out to be the Hankel transform of the function

g : z → 1 z 2 fp (q, β, z), (5.34)
then (5.33) writes as: Cf p (q, β, t) = 2πh 0 g(2πqt).

(5.35)

One can therefore apply the Hankel inverse formula to get

g(z) = 1 2π ∞ 0 Cf p (q, β, t)J 0 (2πqzt)2πqt d(2πqt), (5.36)
where by definition of the function g, we have:

fp (q, β, z) = 2πz 2 q 2 ∞ 0
Cf p (q, β, t)J 0 (2πqzt) dt.

(5.37)

Finally, taking the inverse Fourier transform in polar coordinates, f can be recovered for all (x, y, z) ∈ R 2 × R + as:

f (x, y, z) = 2πz 2 2π 0 ∞ 0 ∞ 0
Cf p (q, β, t)J 0 (2πqtz)t dt q 3 e 2πiq(x cos β+y sin β) dqdβ.

(5.38)

Last integral can be expressed in cartesian coordinates related to the Fourier domain (u, v) and in terms of ω that verifies t = tan ω , and then the theorem holds.

Remark 12. A relationship between 3D Fourier transforms of the function f and 2D (ζ, ξ)-Fourier transforms of Cf can be written. It follows by applying Fubini's theorem to (5.31). If F 3 denotes the 3D Fourier transform of f , we then have:

Cf (u, v, t) = 2π 0 F 3 f (u, v, -t(u cos ϕ + v sin ϕ)) dϕ.
(5.39)

This identity can be seen as a version of the Fourier slice theorem for the conical transform and it can be exploited in the numerical computing of conical projections in the direct problem.

The inverse formula (5.28) is clearly the bidimensional inverse Fourier transform of the last inner integral times the factor (u 2 + v 2 ). A filtered back-projection inversion can then be developed.

A filtered back-projection inversion

Back-projection inversions are written down by means of the adjoint transform and the corresponding filter obtained from the inverse formula in (5.28). To define the adjoint transform in the same way as in the 2D case, let first introduce some basic notations: C is defined as an operator from

X = S(R 2 × R + ) to Y = S(R 2 × [0, π 2
]), associated with the L 2 inner products:

(f 1 , f 2 ) X = R R 2 f 1 (x, y, z)f 2 (x, y, z) dxdydz, (g 1 , g 2 ) Y = R 2 π/2 0 g 1 (ζ, ξ, ω)g 2 (ζ, ξ, ω) dωdζdξ.
(5.40) Definition 10. The adjoint transform C † of C is defined as the unique transform from Y to X that satisfies:

∀f ∈ X, ∀g ∈ Y, (g, Cf ) Y = (C † g, f ) X .

(5.41)

In order to derive an expression for C † , we compute

(g, Cf ) Y = R 2 π/2 0 g(ζ, ξ, ω)Cf (ζ, ξ, ω) dωdζdξ = R 2 π/2 0 g(ζ, ξ, ω) ∞ 0 2π 0 1 r f (ζ + r sin ω cos ψ, ξ + r sin ω sin ψ, r cos ω) dψdr dωdζdξ.
(5.42)

With the change of variables x = ζ + tz cos ψ, y = ξ + tz sin ψ and by using Fubini's theorem, we have:

(g, Cf ) Y = R R 2 1 z f (x, y, z) π/2 0 2π 0 g(x -z tan ω cos ψ, y -z tan ω sin ψ, ω) dψdω dxdydz, (5.43)
from where one is able to extract the adjoint transform C † of C in the form:

C † g(x, y, z) = 1 z π/2 0 2π 0 g(x -z tan ω cos ψ, y -z tan ω sin ψ, ω) dψdω.
(5.44)

The adjoint transform is interpreted as a back-projection procedure scaled by a factor 1 z

. Indeed, one assigns to (x, y, z) the values of projections starting at this point and forming a cone towards the detector with an opening angle ω.

For a first rough reconstruction, f (x, y, z) can be approximated through the adjoint transform or equivalently in this case, through a simple back-projection of the output image as f (x, y, z) ≈ C † Cf (x, y, z).

(5.45)

The adjoint transform C † is not the inverse operator of C. Actually, the resulting expression of C † Cf (x, y, z) and the inversion formula of C (5.28) given in section 5.2.2 differs only by some factor, or filter that must be applied to this adjoint transform before the back-projection procedure.

Theorem 10. Let (u, v) ∈ R 2 and let the projections C * f be defined by filtering Cf in the 2D (u, v)-Fourier space by the filter

C * f (u, v, ω) = (u 2 + v 2 ) sin ω cos 3 ω Cf (u, v, ω). (5.46) Then f can be recovered for all (x, y, z) ∈ R 2 × R + by f (x, y, z) = z 2 π/2 0 2π 0 C * f (x -z tan ω cos ψ, y -z tan ω sin ψ, ω) dψdω.
(5.47)

Proof. C † Cf (x, y, z) is obtained from (5.41) by inserting the 2D (ζ, ξ)-Fourier transform of C to have

C † Cf (x, y, z) = 1 z π/2 0 R 2 Cf (u, v, ω)e 2πi(ux+vy)
2π 0 e -2πiz tan ω(u cos ψ+v sin ψ) dψ dudvdω.

(5.48)

This last inner ψ-integration can be performed via polar coordinates as in (5.32-5.33) allowing us to introduce the first kind Bessel function of order 0 as:

C † Cf (x, y, z) = 2π z R 2 e 2πi(ux+vy) π/2 0 Cf (u, v, ω)J 0 (2πz tan ω u 2 + v 2 ) dω dudv. (5.49)
Last expression only differs from the inverse formula (5.28) by factors z 3 (u 2 +v 2 ) sin ω cos 3 ω . These factors are seen as a filter that we need to apply to data in order to have an exact back-projection inversion procedure.

Finally, we are also able to write the inverse formula of the conical Radon transform by applying the inverse Fourier transform on the last expression giving us the filtered projections C * and through the adjoint transform C † :

f (x, y, z) = z 3 C † C * f (x, y, z), (5.50)
and the theorem yields.

Reconstruction from the spectral image

The scattered photon flux density for all detecting regions belonging to a bounded plane detector forms the spectral image or the output of the imaging system. This spectral image is resolved in energy whereas we need an angular sampling to perform back-projections. The procedure is exactly the same as in the 2D case, that is, from the spectral image Finally, the angular density must be corrected by the only approximation of the DDCS we dispose, that is the Klein-Nishina DCS, the solid angle factor and the sin factor of the measure over conical surfaces not included in the conical transform, that is we have the corrected angular density as

I(C i , D)
1 τ I 0 σ KN (ω) sin ω cos ω dI dE (E, D) dE dω dω.
(5.52)

The conical back-projection inversion assumes that the Dirac approximation of the Compton scattering Klein-Nishina DCS holds. We will see in the next chapter that the Doppler broadening and binding effects will actually blur these reconstructions and another correction method is actually needed for examle by a deconvolution of signal caused by the Doppler broadening and also iterative reconstructions.

Conclusion

The direct and inverse problem for a Compton Scattering Tomography imaging system using an external source has been mathematically investigated. An approximation of the direct problem in which the input function represents the electronic density distribution of free electrons at rest on a bounded region of R 3 has first been presented. In such case, the output spectral image can be modeled by means of the conical Radon transform and the Klein-Nishina differential cross section (DCS). To be physically more sound, and in particular when the heterogeneity of the input object is considered, the direct problem must take into account a more realistic differential cross-section that varies as a function of the scattering angle, scattered energy and the material composition, termed doubledifferential cross-section (DDCS).

The inverse problem has also been addressed and has led to an inverse formula of the conical Radon transform in terms of a filtering back-projection procedure. The Klein-Nishina DCS hypothesis is maintained in this case, as it provides the expectation value of the real DDCS and as it allows to perform conical back-projections. Of course Doppler and binding effects will have a negative (blurring) impact in reconstructions, having larger incidence when the scattering angle is higher or the element is heavier. These effects will be detailed in the numerical simulations proposed in the next chapter. This section presents numerical simulations of both direct and inverse problems in a 3D framework. The discrete formulation of both problems is given with its related spatial and spectral sampling details used in simulations. The corresponding setup is designed in Figure 1.4 and described in section 1.3.

Simulations of the direct problem are first performed by means of the conical transform in (5.4) to compute the scattered photon flux per energy channel and per detecting site given in (5.18). This method assumes an exact relationship between the scattering energy E ω and the scattering angle ω according to (1.1). Through Doppler broadening effect and electron binding energy variations due to the chemical properties of the different materials presented in an input object, this relationship is indeed a non-Dirac probabilistic law. In this work it has been verified that the broadening effect can be approximated with a Cauchy distribution. That is, for a given ω, this Cauchy density distribution has the location parameter at the Compton energy line E ω and the scale parameter depending on the material. Simulations through this law are presented as well.

Realistic Penelope Monte Carlo simulations considering photoelectric absorption, Rayleigh scattering and multiple Compton scattering are presented to assess the consequences of our approximation on image formation and volume reconstruction.

Reconstruction or inversion of the Conical Radon transform is performed by the filtered back-projection procedure summarized in (5.50). An iterative reconstruction process which, whilst more expensive computationally, enable us to correct some artefacts produced by the mentioned effects is a nearly perspective work.

The discrete formulation

For the Compton scattering imaging system considered, the input object, its approximative reconstruction and the output image studied in a continuous point of view in the previous 95 chapter need to be sampled over some grid to perform simulations. Both the direct and inverse problem need also to be written in a discrete version. The following sections provides this formulation.

The discrete input and output spaces

The input object of the system is indeed a 3D array containing sampled information about electronic densities. The array has dimensions of I × J × K. It will be uniformly sampled over a grid {f

(x i , y j , z k ), 0 ≤ i < I, 0 ≤ j < J, 0 ≤ k < K}, (6.1)
where the cartesian coordinates (x i , y j , z k ) assigns the middle of a voxel belonging to the discrete input object of volume ∆x∆y∆z. Note that z k > 0 for all k.

The reconstructed object is sampled on the same grid. The output image or the spectral image is obtained first from conical projections sampled both in opening angles and vertex positions, i.e., through a sampling of the (ζ, ξ, ω) domain, in the form

{Cf (ζ m , ξ n , ω l ), 0 ≤ m < M, 0 ≤ n < N, 0 ≤ l < L}, (6.2)
for a discrete M ×N detector and L opening angles of the cone. The 2D detector is sampled uniformly where each pixel D mn has dimensions ∆ζ∆ξ with cartesian coordinates of its middle point (ζ m , ξ n , 0). The angular sampling step is also uniform and will be noted ∆ω. Finally the output spectral image is resolved in energy with a sampling rate ∆E and then it is an array

{I(D mn , C i ), 0 ≤ m < M, 0 ≤ n < N, 0 ≤ l < L ′ }, (6.3)
where L ′ is the number of energy channels of a detector.

The discrete direct problem

From an input discrete object, the discrete conical transform or projections has to be computed first. It can be written as a sum approximation of the conical transform in (5.4), i.e., as

Cf (ζ m , ξ n , ω l ) = ∆ψ∆r R-1 i=0 P -1 j=0 1 r i f (ζ m +r i sin ω l cos ψ j , ξ n +r i sin ω l sin ψ j , r i cos ω l ), (6.4)
where ∆ψ = π 2P is the azimuthal angular sampling rate into P points, and the parameter r is discretized into R points with a sampling rate ∆r.

From the uniform sampling (6.1) of the input object, the sampled values of f inside the last sum f (ζ m + r i sin ω l cos ψ j , ξ n + r i sin ω l sin ψ j , r i cos ω l ) (

are obtained by a trilinear interpolation of the object as these positions do not coincide in general with those of the grid (6.1). Nearest neighbour interpolation may be an other interesting solution as it would reduce the computational complexity of the computing of the projections.

From the discrete conical projections {Cf (ζ m , ξ n , ωk )}, the spectral image is obtained from integral (5.18), it is approximated, for a given pixel D mn in the 2D detector having its middle point coordinates (ζ m , ξ n , 0), by the sum

I(C i , D mn ) = τ I 0 ∆ω k: Eω k ∈C i sin ωk cos ωk σ KN (ω k )Cf (ζ m , ξ n , ωk ), (6.6)
where E ω k is the energy given by the Compton equation and ω = min{ω, π -ω} whether a forward or backward setup is considered. The Compton DDCS σ(ω, E) is approximated by σ KN (ω k ) if the electrons are supposed to be free. In the next section we approximate this DDCS by a non-Dirac law.

The last computation of the spectral image is performed with an uniformly angular sampling on every energy channel. It means that we need to interpolate the input object this number of angles per energy channel, per number of detecting sites over the grid sampling the conical surface. Moreover it works only when the hypothesis of free-electrons is considered because of the diffeomorphism angle-energy that allows us to sum projections over energy or angular channels equivalently.

In the next section another numerical scheme is developed where the summation on the energy channels is performed via the cumulative distribution function of a law the model the DDCS and then no sampling over each energy channel is performed. Clearly in the free-electron case this law is modeled by the Dirac distribution over the Compton line where the related cumulative function is the Heaviside step function applied on the boundaries of each energy channel.

The differential cross-section for heterogeneous objects

The analytic formulas of the double differential cross-section σ(ω, E) obtained from [START_REF] Salvat | PENELOPE, a code system for Monte Carlo simulation of electron and photon transport[END_REF] and plotted in section 5.1.3.b are going to be approximated by a Cauchy distribution, also called Lorentzian law. This DDCS is a function of the scattering angle and energy and having its maximum over the curve given by the Compton equation ω = E ω . For a fixed ω, it is bell shaped as a function of E having its expectation value at E ω . This 1D curve can be modelled by a Lorentzian law where the spreading factor γ, called the scale of the law, which specifies the half-width of the curve depends on the fixed ω and on the material. γ is then the objective variable in a discrete non-linear optimization problem that minimize the mean square error between the law and the analytic DCS per ω. The Lorentzian law fits better than the Gaussian law due to its more picked bell shape.

The density function of a Lorentzian law is given by

f l (x; x 0 , γ) = 1 πγ 1 + x -x 0 γ 2 -1 . (6.7)
The amplitude of the law is 1 πγ located at its median value x 0 . The related cumulative distribution function

F l (x; x 0 , γ) = 1 π arctan( x -x 0 γ ) + 1 2 . ( 6.8)
The factor γ must be optimized per different element or material and per scattering angle of the photon. As seen in section 5.1.3.b, larger angles induces a larger spreading of the DDCS. Figures 6.1 shows the DDCS for aluminium compared to hydrogen in the ω -E plane where a higher spreading of the law is appreciated for aluminium and for hydrogen, for ω ∈ [ π 2

, π]. It can be verified that heavier elements presents higher spreading coefficients. The increase of the spreading according to the scattering angle is also appreciated. Therefore, for a given element or material the Compton double differential cross section approximated by the Lorentzian law centred at E ω and the scale γ depending on the scattering angle reads σ(ω, E) = Zσ KN (ω)f l (E; E ω , γ ω ), (

where Z is the number of electrons in the atom. The factor γ ω is estimated from the real DDCS by fitting a non-linear model minimizing the least square error between the actual DDCS and the Lorentz density f l (E; E ω , γ ω ). For all scattering angles ω in simulations, the related scale γ ω will be estimated for different materials presented in the test input phantom in section 6.3. It is verified that γ increases regularly with the scattering angle. For a given material, γ ω is estimated by fitting a square polynomial model γω = a + bω + cω 2 , (

where (a, b, c) are the coefficients minimizing the least squares error between all γω and all the γ obtained from the Lorentz distribution estimation for a given sampling of ω. Therefore, the input information to the system will be the electronic density per voxel, and also these coefficients (a, b, c) also per voxel. Figure 6.3 shows both γ and γ as a function of ω for the material used in the input phantom, namely, skin, paraffin wax, lead oxide and aluminium oxide.

Finally, for a given energy channel C i , the spectral image is modelled by (5.11) and (5.18) where the integral over the conical surface is no longer available due to the non Dirac Compton DDCS. For a given detecting site, the integral in (5.11) is then performed not π for aluminium and its approximation by a Lorentzian low wit spreading factor γ = 495 eV which specifies the half-width at half-maximum. E 0 = 50 keV. over the conical surface but over R 2 × R. Per each couple detecting site -scattering site, the scattering angle is computed and the DDCS evaluated consequently.

Therefore, for a given material, the output spectral image in integral (5.18) is rewritten as

I(C i , D mn ) = E∈C i dI(E, D) = τ I 0 R 2 ×R + 1 r 2 M cos ωM σ KN (ω M ) E∈C i f l (E; E ω M , γ ω M ) dE f (M) dM = τ I 0 R 2 ×R + 1 r 2 M cos ωM σ KN (ω M ) F (E i+1 ; E ω M , γ ω M ) -F (E i ; E ω M , γ ω M ) f (M) dM. (6.11)
Last integral is solved numerically as the sum of the contribution of each voxel M ijk having its centre at coordinates (x i , y j , z k ) in the input discrete grid to a given detector pixel D mn , that is, (6.12) where r ijk and ω ijk are the distance from D mn to M ijk and the related scattering angle respectively.

I(C i , D mn ) = τ I 0 i,j,k 1 r 2 ijk cos ωijk σ KN (ω ijk ) F (E i+1 ; E ω ijk , γ ω ijk ) -F (E i ; E ω ijk , γ ω ijk ) f (M ijk ),
The fundamental difference of last sum with sum (6.6) lies on the computation of the spectral image per each couple voxel-pixel directly from a marginalisation of the energy and then the use of the cumulative distribution of the approximative law of the DDCS. This law can either be a Dirac law in the Klen-Nishina case or a bell shaped law when the DDCS is employed. By doing this the integral of the object is based on the cartesian integral of the transform (5.22) while the sum (6.6) is based on the conical integral (5.4).

The discrete inverse problem

From the output spectral image I(D mn , C i ), the aim is to approximate the input object over the same grid in (6.1). As the inverse formula of the conical Radon transform and related filtered back-projections are going to be used, the Dirac relationship between the scattering angle and energy given by the Compton equation is going to be assumed as being the expectation value of the non-Dirac true relationship. Therefore, the Klein-Nishina Compton interaction cross-section is also assumed.

From the spectral image I(D mn , C i ) we need to estimate an angular density to perform the inverse process. As done in the 2D case, this is done by first considering this spectral image as a piecewise constant function over each energy channel, then we sample the scattering angle with an uniform angular step from where an energy sampling {E k } is obtained through the Compton equation.

The discrete energy density is obtained from the spectral image by nearest neighbour interpolation of this image over the energy sampling {E k }.

Finally the angular density is estimated from this energy density {E k } via the change of variable given the Compton equation as

{E dE dω } k . (6.13)
This is the angular density to be use in the inverse problem. Notice that the angular factors in the spectral image sin ω cos ω σ KN (ω) (6.14) need also to be corrected in the last estimated angular density before performing the backprojections.

6.1.4.a Filtering and back-projections

From some angularly resolved data, we can apply the filtering back-projection procedure related to the conical transform. We need the discrete version of the back-projecting operator and of the filtering process.

The discrete adjoint transform that performs back-projections can then be written by using the trapezoidal rule to solve integral (5.44). For some dataset {g(ζ m , ξ n , ωl )} the adjoint discrete operator reads

C † g(x i , y j , z k ) = 1 z k ∆ψ P -1 p=0 L ′ -1 l=0 ∆ω l g(x i -z k tan ωl cos ψ p , y j -z k tan ωl sin ψ p , ωl ) , (6.15)
where ∆ψ = π 2P is the azimuthal angular sampling rate into P points and ∆ω l is the angular step related to the energy channel C l .

Linear interpolation of the dataset {g(ζ m , ξ n , ωl )} is used again here to compute the values g(x i -z k tan ωl cos ψ p , y j -z k tan ωl sin ψ p , ωl ). (

Before back-projecting, the conical projections need to be filtered in the 2D Fourier domain related to variables (ζ, ξ) with the corresponding frequency variables (u, v) by the filters

(u 2 + v 2 ) sin ω cos 3 ω , (6.17)
and by the filter z 3 in the object domain. Note that due to the presence of the factor 1 z in the adjoint transform, this filter is simplified in numerical implementations to z 2 . Therefore, to perform numerically the inverse process from an output dataset I(D mn , C i ) of the imaging system, the procedure works as follows.

1. Compute the conical projections {Cf (ζ m , ξ n , ωl )} sampled angularly and corrected by factors (6.14).

2. Filter the projections in the Fourier domain to get C * f applying (5.46).

Perform an inverse

Fourier transform to obtain the filtered projections C * f . 4. Then the back-projection procedure is performed by means of the adjoint transform C † by (6.15).

5.

Finally, the factor z 3 is applied and we obtain an approximative reconstruction of the input object.

A discrete version of reconstruction formula (5.2.3) is then given, for all (x i , y j , z k ) by Figure 6.4: The filter (u 2 + v 2 ) and its smoothing by two cosine apodization windows.

f (x i , y j , z k ) = z 3 k C † C * f (x i , y j , z k ). ( 6

6.1.4.b Window functions

As in the 2D case, we can control numerically high frequencies dues to the divergence of the filter (u, v) → (u 2 + v 2 ). This effect amplifies contrasting features in reconstructions but also amplifies numerical noise in the spectral images. This effect can be controlled by an apodization window function (u, v) → w(u, v). The resulting filtered projection at ω now reads

C * f (u, v, ω) = w(u, v)(u 2 + v 2 ) sin ω cos 3 ω Cf (u, v, ω). (6.19)
Cosine related window functions are going to be used in our simulations, they are defined on its 2D form by

w(u, v) = cos n (π √ u 2 + v 2 ) , √ u 2 + v 2 ≤ 1 2 0 , √ u 2 + v 2 > 1 2 .
(6.20)

Higher values of n makes the function decrease faster and the filter has then more effect over the (u 2 + v 2 ) filter. On the other hand, a high value of n might over-smooth the reconstruction and therefore some contrasting details in images will be lost. The choice of n in the cosine window, depends on the image configuration and on the expected results. Figure 6.4 shows the behaviour of the filter (u 2 + v 2 ) as well as its smoothing by two cosine windows with n ∈ {1, 2}. Others 2D filters and more details about apodization windows and filtering can be found in [START_REF] Easton | Fourier Methods in Imaging[END_REF].

First simulations, the 3D Shepp-Logan phantom

For our first simulation work, the 3D version of the classical Shepp-Logan phantom containing different ellipsoids that represent the human head [START_REF] Shepp | Computerized tomography and nuclear magnetic resonance[END_REF] is going to be considered as the input object. It is an array of dimensions 64 × 64 × 64 voxels. We will used it to support the forward and inverse conical Radon transform and to see its limitations and numerical artifacts in the reconstructed phantom.

Therefore, neither energy nor spatial resolution are going to be taken into account, the sampling spatial step is given in voxels for the object and in pixels for the detector. The conical projections are uniformly sampled angularly where the Klein-Nishina DCS and the solid angle factor under-which the pixel sees the voxel are not included in the direct problem of these first simulations. Doppler broadening, binding effects, multiple scattering and attenuation are then neglected.

The 2D plane detector is located over the xOy-plane having dimensions of 192 × 192 pixels. The object is at a perpendicular distance of 10 pixels to the detector following the z-axis, i.e., the z-coordinate of the input Shepp-Logan phantom verifies z > 10. As we don't have a hole producing incomplete data, there is no need to place the phantom farther as seen in the 2D case. The beam has the same square cross-sectional area of the frontal surface of the phantom, then no translations are performed.

Sampling parameters

The discrete conical Radon transform in (6.4) requires a spatial sampling over the conical surface that is set here to ∆r = 1.2 times the voxel length for the slant heigh and ∆ψ = 0.01 for the azimuthal angular step in radians. Different cones are considered by sampling the polar angle ω by 78 uniform opening angles of cones from 0 to π 2 , i.e., ∆ω = 0.02 radians. With this sampling choice we have an output dataset of dimensions 192 × 192 × 79 to reconstruct an object over a grid of 64 × 64 × 64. A voxel in the phantom has then the same length to a pixel in the detector.

The same angular sampling steps are chosen to compute the discrete adjoint transform in (6.15). A soft apodization cosine window with n = 1 is applied in the filtering process due to the numerical noise and the non applying of random Poisson noise.

Simulation results

Figure 6.5 shows the conical projections of the phantom for three different opening angles, namely ω ∈ { 1 12 π, 1 4 π, 5 12 π} radians. Low angles data are almost fully covered by the detecting plane whilst data related to large angular values will be lost due to the limited detector area. The filtering process will act on the 2D Fourier transform of each plane for all ω fixed. Figure 6.6 shows the filtered projections for the same angles, we can appreciate the increase of contrasting details in projections due to the filter (u 2 + v 2 ). The quality of the filtered back-projection inversion for some noise-less data related to free electrons at rest depends clearly on the sampling and dimensions. Although contrasting features are well reconstructed at the correct positions. Reconstruction results are shown in Figure 6.7 for cross-sections given by planes parallels to the xOy-plane, that is, to the detector, planes related to a larger value of z presents then more artifacts due to the lack of data back-projected from bordering regions of the detector. Cross-section given by planes parallels to the xOz-plane are also presented in Figure 6.8 where we can also appreciate more artifacts in reconstructions in zones that are farther to the detector.

As noise was not added, higher values of n for the cosine apodization window are not necessary, these results are not presented as we are investigating an ideal scenario with a phantom constituted of elements showing regular boundaries as the Shepp-Logan phantom is. The correlation coefficient of 0.70 is obtained while the variance ratio is 1.27. Reconstructions with noise data and different apodization windows can be found in [6].

6.3 A 3D stratigraphical analytical phantom 6.3.1 The phantom A 3D phantom is constructed in a way that allows us to compare the deterministic approach of the imaging problem presented via the conical Radon transform to Monte Carlo simulations performed via the code Penelope [START_REF] Salvat | PENELOPE, a code system for Monte Carlo simulation of electron and photon transport[END_REF]. This code offers constructions of scenes based on potentially degenerated quadric surfaces, for phantoms as well as for detectors. The test phantom is constructed correspondingly.

The phantom, corresponding to a flattened stratigraphic sample is analytically constructed through planar surfaces with two layers of 16 × 16 µm 2 regarding its frontal area and about 64 and 48 µm thickness each layer. These layers are placed parallel to a 2D detector. The closest layer to the detector presents two spheres of 8 µm radius, while the further one contains one sphere of 14 µm of radius.

The object is placed at a perpendicular distance of 50 µm from the detector, the 2 layers are then limited by 3 planes parallel to the detector located at 50, 114 and 162 µm respectively from the detector following the z-axis. A rotation of an angle θ about the y-axis was done per each layer, namely by the angles θ = {3.2, -3, 2.8} degrees respectively.

Regarding the chemical composition of the phantom, the spheres represent inorganic grains within the two organic layers, each value in the phantom describes the electron density of the considered material. In this section, the electron densities are the number of electrons per cm 3 .

The layer closest to the detector represents an organic layer of skin, having an electron density of 3.638981 × 10 23 . The other layer is composed of paraffin wax, with electron density 3.207726 × 10 23 . The two grains in the first layer are grains of aluminium oxide and the third one, in the second layer, represents a grain of lead oxide. These materials have electron densities of 1.172404 × 10 24 and 2.314164 × 10 24 respectively. All these electronic densities are taken from [START_REF] Salvat | PENELOPE, a code system for Monte Carlo simulation of electron and photon transport[END_REF].

The sample is represented in Figures 6.12 and 6.16, giving cross sections parallels to the xOz and xOy planes respectively.

Dimensions and sampling details

The unit length considered in the 2D detector and in 3D objects is 1 µm. The spatial resolution of the detector is thus 1 µm 2 represented by a pixel in the spectral image and a voxel in objects represents 1 µm 3 . Dimension parameters are set according to simulations results performed in the 2D case.

The sample is therefore represented in a 16 × 16 × 120 array. The setup is configured as follows. A square parallel X-ray beam of 8 × 8 µm 2 of sectional area crosses a hole in the detector of 10 × 10 µm 2 (ζ 0 = ξ 0 = 10 µm in Figure 1.4) if a backward configuration is considered. The incident energy of the beam is 50 keV. As the simulated sample has the mentioned dimensions, then we need 1 translation of 8 µm following the 0x-axis and the same for the Oy-axis, to cover the full sample when no intersections between the two imaging process are performed. As noticed in the 2D case, a specific loss in scale is presented in boundaries voxels perpendicular to the detector, then we performed translations of 4 µm instead following the axis x and y and we take only the 4 × 4 × 120 central voxels of the reconstructed slice. The detector plane, located over the xOy-plane, has an area of 320 × 320 µm 2 including the hole. The number of detecting sites, having an area of 1 × 1 µm 2 each as mentioned, is 89900 considering the absence of detecting sites in the hole area.

Spectral image formation

Two numerical schemes are developed to generate the spectral image, the first is based on the definition of the transform given in (5.4) that follows the same principle employed in the 2D case where an uniform angular sampling is considered. The output spectral image is obtained by integrating over the energy channel as in 6.6. Due to the complexity of computations, this version is only suited for the free electrons case where the angular DCS given by the Klein-Nishina formula is employed. We need to perform interpolations of the phantom over each sampling of the cone per both each angle and each detecting site. If the DDCS is considered, another loop must be added related to sampling of the Lorentz law per scattering angle.

That was the reason to develop the second image formation scheme. Based now on the definition (5.22) of the conical transform, the spectral image is computed by the sum (6.3) for each discrete energy channel. No angular sampling is done and then no need of interpolations over conical surface. The sum is performed over every voxel in the input object for each pixel in the detector but instead of interpolations, only a distance and angle are calculated. Using this numerical scheme, a GPU parallelization is possible due to the cartesian nature of the integral. The DDCS can be approximated by a law where the cumulative distribution function is known, as the Lorentzian law for the non-free electrons case or the Dirac law for the free-electrons one.

For the first scheme, the conical one, cones in (6.4) are sampled by ∆r = 1.2 µm for the slant heigh and ∆ψ = 0.01 radians for the azimuthal angle. 157 cones are simulated per detecting site then the cone opening sampling rate is ∆ω = 0.01 radians from 0 to π 2 .

The spectral resolution of detectors is set to 5 and 50 eV in the simulations. Then we will dispose of 740 and 74 energy channels in a backward setup sampled from E π/2 to E π .

The central region of the phantom of 8 × 8 µm 2 of sectional area and 120 µm of depth is going to be analysed first. Figures 6.9, 6.10 and 6.11 show respectively the spectral image of the sample at three different energy channels obtained by the two schemas mentioned and a third result through Monte Carlo simulations with the code Penelope [START_REF] Salvat | PENELOPE, a code system for Monte Carlo simulation of electron and photon transport[END_REF]. The three different rings corresponding to each grain are appreciated in the first case due to the absence of broadening in the first scheme.

The second scheme gives the same spectral image to the first one if the Dirac law is applied instead of the Lorentzian within numerical errors, the image is not presented but it is remarkable the difference in computing time. Whilst the configurations are identical, the first scheme requires a couple of hours to proceed while the latter, using GPU acceleration, is finished within less than a minute. Monte Carlo simulations of about 1 × 10 11 incident photons to produce the thirds series of images is even more time consuming, the computation to get those took nearly a full month on a machine with 64 CPUs.

Object reconstruction and conclusions

The phantom is reconstructed through the back-projection procedure detailed in this section. First, filtered back-projections from angular sampled data are performed. No energy resolution is considered and then we can back-project filtered projections directly obtained in the detector without performing a summation into an energy channel. Therefore, the quality of reconstructions doesn't suffer from this averaging of projections and depends uniquely on the angular sampling rate, here ∆ω = 0.01 radians. Horizontal and vertical translations of 4 µm were performed to reduce the edge artifacts. As angular sampled projections are considered, we have for these reconstructions 157 projections.

Then, energy sampled data are obtained through the sum (6.6). The Klein-Nishina and solid angle factors are applied to these projections to generate the energy resolved image. These factors were corrected as explained in last section to the angular density estimated from the spectral image. Then the filtered back-projection algorithm is performed. A cosine window function in the form of (6.20) when n = 1 due to the non presence of noise was used in the Fourier domain of projections. Missing data not recorded in the hole area was filled by a numerical solution of a heat diffusion problem from values of energies captured at the hole edge as showed in Figure 6.20. This allows us to correct partially the specific loss in scale in reconstructions that is caused by the missing data that would have been recorded in the hole area. The effect of these missing data in a 50 eV resolved detector affects mainly the two lowest energy channels, that is, the channels covering energies within [E π , E π + 50[ and [E π + 50, E π + 100[ eV. In the first one, most of the data would be recorded in the hole area then the effect of the missing data is larger. Filling it with numerical heat propagation corrects partially the lost in scale but also the discontinuity in the hole boundary. From the third channel the effect of the hole is negligible in the free-electrons scenario. Of course if a spreading of the signal is present due to the Doppler broadening, even at larger energy channels, some signal will be lost in the hole. First reconstructions are presented for cross sections parallels to the xOz plane, where we can appreciate the quality of reconstructions in recovering the correct depth of different elements. Cross sectional reconstructions parallels to the xOy planes are presented then, that its, for planes parallels to the detector.

Reconstruction from angular sampled projections are presented in Figures 6.13 and 6.17. Figures 6.14 and 6.18 show reconstructions from an energy resolved spectral image by 5 eV and Figures 6.15 and 6.19 for 50 eV. The error measurements are presented in table 6.1 namely the correlation coefficient and the variance ratio for each configuration considered.

A smoothing in reconstructions is appreciated due to the sum of all projections belonging to a given energy channel. This summation also generates a loss in scale in reconstructions appreciated even more at large energy resolutions. We can appreciate that reconstructions from uniform angularly sampled projections are reconstructed almost at the correct scale, this is seen also in the related variance ratio. Given a low energy resolution as 5 eV, the contrasting details of the input object are well reconstructed as appreciated in figures and in their correlation coefficients. However the absolute scale is affected as in the 2D case because at low angles, the summation and the latter approximation of the angular density by nearest neighbours mostly eliminates the picks in the 1D angular spectra per detecting site. This is well appreciated in Figure 4.9. For larger energy resolution the correlation coefficient is also affected because of the loss in details when performing the nearest neighbour approximation of the angular density, mostly for low angles. As in the 2D case, the hole affects mostly in the scale because, first most of the signal would be detected there and second, it corresponds to low angles projections that will be latter averaged on a potentially large energy channel. For all three configurations, the hard grains are detected at their correct depth within the lower density matrix. Half translations help to correct the difference in scales from boundary voxels and inner ones. Reconstructions from non-overlapping translations are not presented since strong artifacts from one translation to another are present.

All these simulations and reconstructions were computed on the basis of the free electron model. When taking into account a DDCS model accounting for Doppler broadening and binding effects. The blurring of the spectral images showed in Figures 6.10 or 6.11 renders the direct application of the filtering back-projection process to those images totally ineffective. Other imaging alternatives may be available to face this issue, some of which are exposed in the conclusions in chapter 7. To illustrate the problematic, a reconstruction considering a realistic scenario including Doppler and binding effects accounted in the DDCS model simulated through a Lorentzian distribution is exposed in Figure 6.21. The reconstructed slice concentrates the signal at the area closest to the detector. The original 8 × 8 × 120 slice is also shown, the different grains and layer boundaries are not well identified.

Conclusion and perspectives

A X-ray imaging modality based on Compton scattering is presented in this work, the main objective is to perform a 3D mapping of flat heritage objects without the need of relative rotations between object, source and detectors. Both the image formation process by means of a conical Radon transform and the object reconstruction by a filtered back-projection procedure are exposed and supported by numerical simulations showing feasibility with real data. Results are already very encouraging considering the problem of non-destructive and non-invasive 3D imaging of samples supported by a deep or dense material. Clearly, the current method should still be worked on to eliminate some artifacts which are strong enough to blur out some of the contrast of the electron density. Yet our methodology enables us to strengthen our detailed understanding of the image formation process and shows how it is coupled to the volume reconstruction itself. This is proving instrumental in designing the physical device that will ultimately enable such imaging modality.

To be able to tackle the system computationally, we have first focused on flat objects and using incident energy that optimise the relative importance of a single inelastic scattering event versus more complex, multiple scatterings and absorption interactions. At the chosen incident energy the probability of multiple scattering is very small since the mean free path of the X-ray photons is in the same order of magnitude as the probed volume thickness. When dealing with single scattering and a monochromatic incident beam, it is easy to discriminate photons generated by elastic / inelastic scattering events. At the selected incident energy the mean free path attached to scattering is much shorter than that attached to the photoelectric absorption process, in other words the effects of photoelectric absorption is much lower than scattering effects. It is important to note that inelastically scattered photons have an energy which is in the same domain as the incident photon energy when it comes to relative importance of absorption and scattering.

For a mathematical understanding of an imaging system, the related direct and inverse problem must be modeled. Different physical parameters describing the device have to be chosen properly as for example spatial and spectral resolutions, dimensions or initial configurations. A good understanding of these parameters in both the output image formation process and the volume input reconstruction are crucial.

A simplified two-dimensional (2D) setup is considered first to answer this questions where the incident beam and diffracted photons are limited to a 2D plane. This limitation is first considered because to the equations are more easy to manipulate, simulations schemas are less complex and time consuming, and the azimuthal scattered angle is uniformly distributed on the interval [0, 2π[ in Compton events. We present the direct and inverse problem related to this 2D scenario based mainly on the V-line Radon transform and its filtered back-projection inversion [START_REF] Morvidone | On the V-line radon transform and its imaging applications[END_REF][START_REF] Guerrero | Modelling of a New Compton Imaging Modality for an In-depth Characterisation of Flat Heritage Objects[END_REF] having also the hypothesis of free electrons at rest. That means that the input object only represents the electron density at a given position, for the same material or element. Numerical simulations support these models where different physical parameters are tested. An interpretation of results according to these parameters is presented to be afterwards used in the complete three-dimensional (3D) framework.

The 3D configuration is presented next, the direct and inverse problem are modelled through mainly the conical Radon transform [6,[START_REF] Guerrero | Three-dimensional imaging of flat natural and cultural heritage objects by a compton scattering modality[END_REF]. The model first assumes the freeelectrons hypothesis and a more realistic scenario where the chemical composition of different materials are taken into account.

The image formation process is dominated by the angle-energy double differential crosssection (DDCS) of Compton scattering in this case. In the free-electrons hypothesis the DDCS is simply given by a Dirac distribution concentrated on the line given by the Compton equation, while in the non-free electrons case, a broadening of this line is presented due to the Doppler effect with some discontinuities due to binding effects of electrons of a given element. The DDCS of the latter case is approximated by a Lorentz distribution located on the Compton line and the broadening depending on both the element and the scattering angle. For a given element, the broadening factor of the Lorentz law varies regularly with the angle and then a square polynomial regression is performed to estimate the factor.

The Lorentz model gives a good approximation of the spectral image formation, this is supported by a realistic Monte Carlo simulation of the image obtained from the same input object through the code Penelope [START_REF] Salvat | PENELOPE, a code system for Monte Carlo simulation of electron and photon transport[END_REF].

The output spectral image obtained through the polynomial approximation is compared with the one given by the Dirac distribution. The broadening is too strong to directly perform a filtered back-projection inversion even if the spectral resolution is low enough.

We envisage some ideas to face this issue in a perspective work. A deconvolution of the output spectral image can be performed. The related convolution kernel induced by the Doppler broadening effect can be estimated directly from the relationship between the simulated spectral images of both the Dirac and Lorentz DDCS. This can then be followed by an iterative reconstruction process as the direct problem is well understood and, from a given model, it provides simulated data that is a close approximation to the observed data, such that an iterative approach to fit the simulated image to the observations should provide a good model. Another approach would be to consider the problem in the framework of a linear least-square problem, then a reconstruction through a pseudo-inverse of the related matrix of the direct problem can also not only address the Doppler broadening issue but also the incompleteness of the data. These works are already in progress.

Conclusion et perspectives

Une modalité d'imagerie par rayons X basée sur le rayonnement diffusé Compton est présentée dans ce travail, l'objectif principal est d'effectuer une cartographie 3D des objets plats du patrimoine sans mener des rotations relatives entre l'objet, la source et les détecteurs. Le processus de formation d'image basée sur la transformée de Radon conique et la reconstruction de l'objet par une procédure de rétro-projection filtrée sont exposés et supportés par des simulations numériques montrant la faisabilité avec des données réelles. Les résultats sont déjà encourageants au vu du problème complexe d'imagerie 3D non destructive et non invasive d'échantillons potentiellement supportés par un matériau dense. Clairement, la méthode actuelle doit encore être nourrie avec des méthodes de correction de certains artéfacts qui affectent le contraste de la densité électronique reconstruite. La méthodologie présentée nous permet déjà de renforcer notre compréhension détaillée du processus de formation d'image et montre comment il est couplé à la reconstruction du volume. Ceci fournira les outils nécessaires pour la conception de l'instrument physique qui permettra finalement cette modalité d'imagerie.

Pour aborder le système numériquement, nous nous sommes d'abord concentrés sur des objets plats avec une énergie incidente qui optimise l'importance relative d'un événement de diffusion inélastique par rapport à des interactions multiples et d'absorption plus complexes. À l'énergie incidente choisie, la probabilité de diffusion multiple est très faible du fait que le libre parcours moyen des photons est du même ordre de grandeur que l'épaisseur du volume sondé. Lorsqu'il s'agit d'une diffusion unique et d'un faisceau incident monochromatique, il est facile de discriminer les photons générés par les événements de diffusion élastique / inélastique. A l'énergie incidente considérée, le libre parcours moyen lié à la diffusion est beaucoup plus court que celui correspondant au processus d'absorption photoélectrique, en d'autres termes les effets de l'absorption photoélectrique sont beaucoup plus faibles que les effets de diffusion. Il est important de noter que les photons diffusés inélastiquement ont une énergie qui est dans le même domaine que l'énergie des photons incidents par rapport à l'importance relative de l'absorption et de la diffusion.

Pour La configuration 3D est présentée ensuite, les problèmes direct et inverse sont modélisés principalement au travers de la transformée de Radon conique [6,[START_REF] Guerrero | Three-dimensional imaging of flat natural and cultural heritage objects by a compton scattering modality[END_REF]. Le modèle suppose d'abord l'hypothèse d'électrons libres et ensuite un scénario plus réaliste où la composition chimique de différents matériaux est prise en compte.

Le processus de formation d'image est dominé par la double section efficace différentielle en énergie et en angle (DSED) de la diffusion Compton dans ce cas. Sous l'hypothèse des électrons libres, la DSED est simplement donné par une distribution de Dirac concentrée sur la ligne donnée par l'équation de Compton, alors que dans le cas des électrons non libres, un élargissement de cette ligne est présenté dû à l'effet Doppler avec quelques discontinuités dues aux effets de liaison des électrons d'un élément donné. La DSED de ce dernier cas est approché par une distribution de Lorentz située sur la ligne Compton et dont l'élargissement dépend à la fois de l'élément et de l'angle de diffusion. Pour un élément donné, le facteur d'élargissement de la loi de Lorentz variant régulièrement avec l'angle, une régression polynomiale de degré 2 est effectuée pour estimer ce facteur d'élargissement.

Le modèle de Lorentz donne une bonne approximation de la formation d'image spectrale. Ceci est confirmé par des simulations réalistes Monte Carlo de l'image obtenue à partir du même objet d'entrée avec le code Penelope [START_REF] Salvat | PENELOPE, a code system for Monte Carlo simulation of electron and photon transport[END_REF].

L'image spectrale de sortie obtenue par l'approximation polynomiale est comparée à celle donnée par la distribution de Dirac. L'élargissement est trop important pour effectuer une inversion par rétroprojection filtrée même si la résolution spectrale est très fine. Ainsi, le processus développé ici est insuffisant pour effectuer une reconstruction à partir de données entachées de cet élargissement.

Nous envisageons quelques idées pour faire face à ce problème dans un travail à venir.

Une déconvolution de l'image spectrale de sortie peut être effectuée.

A Numerical implementations

The different algorithms developed in this work are going to be detailed here. Simulations are implemented and performed in the R language [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF] together with some of its different libraries detailed below. The C-based programming language OpenCL was employed to implement a GPU parallelization of the algorithms. The MonteCarlo code for simulations of electron and photon transport through matter Penelope [START_REF] Salvat | PENELOPE, a code system for Monte Carlo simulation of electron and photon transport[END_REF] was employed as well to simulate alternatively the direct problem.

The main part of algorithms is the compute of conical projections and their backprojections. Only those parts are detailed here. Concerning the compute of projections regarding the original 2D Radon transform and its different versions in chapter 2 are not treated here as their are already well studied, for example in [START_REF] Toft | The Radon Transform: Theory and Implementation[END_REF]. The V-line transform implementation as being the 2D version of the conical transform is not provided here.

The different R packages used in codes are given in table A.1.

A.1 The conical projections

Given an input discrete object f , the conical projections are obtained from the sum (6.4). We recall that (ω, ψ) are the polar and azimuthal angle of a cone, (ζ, ψ) the vertex or detecting site and r the distance from the vertex to a point in the conical surface. Projections are obtained by first, defining a function cone.proj in R with inputs omega, psi, xyr.grid and phantom that computes all half-line projections with a same (ω, ψ)-direction for all detecting sites (ζ m , ξ n ) included in the input grid xyr.grid. This grid must also contain the samples r i of the half-line characterised by (ω, ψ). That is, the function returns the 2D array

(P(ζ m , ξ n )) m,n (ω, ψ) = R-1 i=0 1 r i f (ζ m + r i sin ω cos ψ, ξ n + r i sin ω sin ψ, r i cos ω). (A.1)
Projections are obtained by trilinear interpolation of the phantom over all the half-lines given by the couple (ω, ψ). Interpolations are performed with the R function approx3d. The sampling steps ∆ψ∆r will be included after all projections are computed, then the projection function is We will apply the cone.proj function to a uniform sampling of (ω, ψ) here, an uniform energy sampling needs a change of variables through the Compton equation and is not detailed here.

The projections related to samples on ω are computed in a CPU parallel process by means of the R function mcmapply. Last, due to a limited RAM memory, the samples of ψ are added one by one through a external for loop. The sampling steps ∆ψ∆r are applied in a last step. Then the conical function is From the conical projections, the spectral image is obtained if the Compton-Klein-Nishina (free electrons) hypothesis is assumed. For non-free heterogeneous objects, we need another loop in the cone.proj function that compute conical projections for a discrete Lorentzian law per each ω. This process is time consuming, hence a GPU parallelization through a cartesian parametrisation of the conical transform was developed in R using the OpenCL interface provided by the OpenCL package.

A.2 The conical back-projections

The conical back-projections are implemented in R by means of the function back.proj having as inputs the projections and the coordinates (x, y, z) of the voxel in the object where the projections will be back-projected and summed over all ω. This function computes per each (x, y, z) the discrete adjoint transform in (6.15) of projections g times the filter z 3 , that is, B(x i , y j , z k ) = z 2 k P -1 p=0 L-1 l=0 g(x i -z k tan ω l cos ψ p , y j -z k tan ω l sin ψ p , ω l ) .

(A.

3)

The sampling steps ∆ψ∆ω are applied outside the function. The function performs bilinear interpolations through the R function interp.surface of projections to estimate the values g(x i -z k tan ω l cos ψ p , y j -z k tan ω l sin ψ p , ω l ) per each ω l inside a for loop on ω. The function is then x.size <-dims [START_REF] Als-Nielsen | Elements of modern X-ray physics[END_REF] y.size <-dims [2] omega.size <-dims [START_REF] Bockwinkel | On the propagation of light in a biaxial crystal about a midpoint of oscillation[END_REF] psi.size <-dim(psi.vector) # detector sampling detector.x <-seq(from = -x.size*0.5+1, to = x.size*0.5, by = attr(projections, 'delta.x')) detector.y <-seq(from = -y.size*0.5+1, to = y.size*0.5, by = attr(projections, 'delta.y')) omega.vector <-attr(projections, 'omega.vector') To recover the full region of interest, the last function is applied over a grid (x, y, z) containing the voxel coordinates of the phantom to recover. This is done in parallel over the set of voxels with the function back.proj. The sampling steps must be applied as well.

# output array

Projections are first filtered in the Fourier domain through R functions as fft that computes the fast Fourier transform. This is not exposed here. The object is finally recovered, for an object of dimensions (Nx, Ny, Nz) and xs, ys, zs the vector of x, y, z coordinates respectively of the phantom with the process # cartesian grid to recover the phantom grid <-expand.grid(xs,ys,zs) colnames(grid) <-c("X","Y","Z") # parallel backprojections reconstruction <-mcmapply(back.proj, grid$X, grid$Y, grid$Z, mc.cores = detect.cores())

# organize the output array and sampling steps reconstruction <-array(reconstruction, dim = c(Nx,Ny,Nz))*delta.psi*delta.omega

A.3 A GPU cartesian parallelization

The conical projections sampled on the opening angle computed in section A.1 are based on the definition of the conical transform (5.4) where integrals are performed over conical surfaces. These surfaces are parametrised by (ψ, r) for different cones given by a fixed opening angle ω and detecting pixel. The algorithm is written to compute conical projections in parallel on ω for different cones. A cartesian parallelization over each pixel in the detector and a cartesian parametrisation of the conical surface written in (5.22) allows us to use a GPU parallelization in the calculation of projections.

If the double differential cross-section is approximated by a probabilistic law with cumulative distribution function F , the spectral image can be computed per each detecting pixel D mn by the sum

I(C i , D mn ) = τ I 0 i,j,k 1 r 3 ijk cos ωijk σ KN (ω ijk ) F (E i+1 ; E ω ijk , γ ω ijk ) -F (E i ; E ω ijk , γ ω ijk ) f (M ijk ). (A.4)
This is done by defining a OpenCL kernel direct.ocl that will be called by R and the OpenCL package which output is the photon flux density over a given energy sampling for each detecting pixel.

The kernel performs 4 for loops, 3 regarding the voxel index in the phantom and a last inner loop for each energy channel. For every voxel in the input object, the scattering angle is calculated, and then the energy recorded in the pixel related to this angle or voxel is determined by the double differential cross-section and the related solid angle. The DDCS is approximated by a Lorentz probability law for every material and angle as explained in chapter 6. The kernel follows. The kernel is executed in R with the R function OCLRun. The function input is the kernel doppler.ocl pre-compiled previously by the R function oclSimpleKernel and all the dimensions and sampling parameters. We define a function simulate.direct with inputs the kernel, input phantom and configuration samplings as Rarrays that reorganize all these arrays in the way that OCl request.

The simulation is divided by splitting the detector into some number of tiles due to the prescence of a watchdog timer in the GPU that kill a process if a simulations takes too much time. The function is If the DDCS including Doppler broadening and binding effects is considered, it is approximated by a Lorentz distribution that is integrated over a given energy channel. The result of the integral is given by the cumulative distribution function (6.8). In the free electrons case, considered in the previous section, the DDCS is done by a Dirac probabilistic distribution concentrated on the Compton line and having a null measure support. By marginalizing the angle by a given energy channel, the related cumulative distribution is the Heaviside step function that is used in a new free-electron OCL kernel which most significant change according to the last kernel is in lines 161 -175 that changes to These non-free electron simulations were obtained through a cartesian integration of the simulated sample using a GPU implementation. . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Paint cross-section showing a stratigraphical assemblage of The Anatomy Lesson of Dr. Nicolaes Tulp, 1632 by Rembrandt, Mauritshuis, The Hague. c Sample taken and prepared as cross-section by P. Noble during the conservation treatment of the painting in 1997, and re-photographed by A. van Loon, Mauritshuis, in 2010 for the Rembrandt Database.

Figure 1 . 2 :

 12 Figure 1.2: A flat fossil actinopterygian over a high thick X-ray absorbing support from the Kem Kem Beds in Morocco dated back to the Lower Cretaceous (95 million years ago). c P. Gueriau (MHNM/MNHN).

Figure 1 . 3 :

 13 Figure 1.3: A single Compton scattering event. A photon of energy E 0 is scattered at M and arrives at D with energy E ω .

Figure 1 . 4 :

 14 Figure 1.4: The general imaging configuration device. A scattering site M produces scattered radiation captured at a detecting site D. The incident photon is coming from a parallel beam with a square cross-section.

Figure 1 . 5 :

 15 Figure 1.5: The 2D imaging configuration. A scattering site M produces scattered radiation by an angle ω = π -ω captured at a detecting site D. The incident photon is coming from a parallel beam. Three vertical translations have been performed to irradiate the whole object.

[ 1 ]

 1 P. Guerrero, M. K. Nguyen, L. Dumas, S. X. Cohen. Three-dimensional imaging of flat natural and cultural heritage objects by a Compton scattering modality. Journal of Electronic Imaging 26(1):011026, 2017. http://dx.doi.org/10.1117/1.JEI.26.1.011026

. 4 ,

 4 un photon incident d'énergie E 0 , dans la direction de l'axe Oz, est diffusé par effet Compton par un électron situé à l'intérieur d'un objet au point M sous-tendant un angle ω qui varie dans l'intervalle π 2 < ω < π avec la direction de l'incidence. Le photon diffusé, d'énergie E ω approximée par (1.1), atteint un point de détection D = (ζ, ξ, 0) sur le détecteur 2D qui est situé sur le plan xOy. Comme représenté dans la figure 1.4, le détecteur sera placé entre la source et l'objet pour capturer des photons rétro-diffusés. Un trou sera donc présent au milieu du détecteur possédant une surface de 4ζ 0 ξ 0 pour deux nombres réels positifs (ζ 0 , ξ 0 ) pour permettre la circulation du faisceau. Par conséquent, nous aurons des données pour des valeurs de (ζ, ξ) sur le plan xOy vérifiant |ζ| > ζ 0 et |ξ| > ξ 0 . Des translations horizontales et verticales de l'échantillon pourront être nécessaires pour permettre la radiation de l'objet entier. Lorsque tous les électrons du matériau sont à la fois libres et au repos, l'équation de Compton établit dans (1.1) un difféomorphisme entre l'angle de diffusion et l'énergie diffusée, à condition que le faisceau incident soit monochromatique. Dans de telles circonstances, le processus de formation d'image associé à la TCD correspond mathématiquement à la transformation de Radon sur des surfaces coniques, appelée la transformée de Radon conique. Une configuration 2D va être étudiée préliminairement dans l'intérêt d'avoir des calculs simplifié dans les simulation numériques. Ceci est dû au fait que l'angle de diffusion azimutale est uniformément distribué dans un événement Compton, c'est-à-dire que la distribution angulaire des photons diffusés a une symétrie axiale autour de la direction d'incidence. Ce cadre est représenté dans la figure 1.5. Ici, le détecteur est situé sur l'axe Oy du plan cartésien de référence et noté D ζ = (0, ζ). Ce détecteur 1D est placé entre la source et l'objet pour capturer des photons rétro-diffusés comme dans le cas 3D et il aura un trou de longueur 2ζ 0 pour permettre au faisceau de traverser le détecteur. Nous aurons donc des données pour des valeurs de ζ vérifiant |ζ| > ζ 0 pour ζ 0 > 0. Des translations verticales de l'échantillon seront également nécessaires. Dans ce cas, l'intégration est effectuée sur deux lignes dans un demi-espace, conduisant à ce que l'on appelle la transformation de Radon en V.

  Data (projections).Reconstruction(c) Reconstruction from (b).

Figure 2 . 1 :

 21 Figure 2.1: The reconstruction problem related to a parallel source.

Figure 2 . 2 :

 22 Figure 2.2: Geometrical representation of the original (p, ϕ)-Radon transform.

Figure 2 . 3 :

 23 Figure 2.3: Geometrical representation of the (u, v) cartesian Radon transform. Definition (2.8) of the Radon transform in terms of coordinates (p, ϕ) ∈ R × [0, π[ can be reformulated with cartesian coordinates (u, v) ∈ R 2 via

3 Figure 2 . 4 :

 324 Figure 2.4: Geometrical representation of the half-space Radon transform for three different angles ω i and a given starting position ζ.

Definition 5 .

 5 Let f ∈ S(R + × R). The half-space Radon transform, denoted by H, is defined, for -∞ < ζ < ∞ and -

  Half-space Radon transform.

Figure 2 . 5 :

 25 Figure 2.5: Point Spread Function supports related to the three different Radon transforms for an impulse function located at (1, 1).

Figure 2 . 6 : 3 Figure 2 . 7 :

 26327 Figure 2.6: Test object, the Shepp-Logan phantom.

  Radon space (projections).

Figure 2 . 8 :

 28 Figure 2.8: The original Radon transform, reconstruction via the filtered back-projection procedure.

Figure 2 . 9 :

 29 Figure 2.9: The cartesian Radon transform, reconstruction via filtered back-projections.

Figure 2 . 10 :

 210 Figure 2.10: The half-line Radon transform, reconstruction via filtered back-projections.

π 2 [ 3 Figure 3 . 1 :

 2331 Figure 3.1: Geometrical representation of the V-line Radon transform.

. 14 )Figure 3 . 2 Figure 3 . 2 :

 143232 Figure 3.2 shows the support of PSF related to the V-line transform for an impulse input function at (20, 0). It gives an idea of the image shape in the Radon space of an input function obtained as a decomposition of delta functions.

Figure 3 . 3 :

 33 Figure 3.3: A single scattering event from M detected at D ζ with a solid angle dΩ.

) 3 .

 3 1.3.b Scattered photon flux angular density Denote dI(ω, ζ) the recorded scattered photon flux density, that is the number of photons of energy laying in ωdω recorded per unit time at D ζ .If a scattering site M is fixed, let's first write this density only for photons scattered at M, and call it dI(ω, ζ|M). It incorporates the following parameters:

  The solid angle Ω(M, D ζ ) can be seen from Fig. 3.3 to be Ω(M, D ζ ) = 2 arctan τ 2r cos ω , (3.18) where τ is the length of the detecting element located at D and r the Euclidean distance from M to D. If τ is small enough, then Ω(M, D ζ ) can be approximated by τ r cos ω. Therefore, the scattered photon flux density at D ζ , for a fixed site M, is given by dI(ω, ζ|M) = τ r I 0 cos ω σ 2D (ω) f (M) dω dM. (3.19) From the last equation, we can write the scattered photon flux density dI(ω, ζ) recorded at D ζ . It is the integral of (3.19) over all scattering sites laying on the V-line that starts at D ζ towards the object with opening angle ω, called V ω,ζ and represented in Figure 3.1. This density is hence given by the integral dI(ω, ζ) = M∈V ω,ζ dI(ω, ζ|M), (3.20) that can be expressed with the V-line Radon transform as dI(ω, ζ) = τ I 0 cos ω σ 2D (ω)Vf (ω, ζ) dω. (3.21)

Figure 4 . 1 :

 41 Figure 4.1: Parameters and setup description of the 2D backward Compton scattering tomography.

  : copper sample, with segmentation in 6 classes.

Figure 4 . 2 :

 42 Figure 4.2: First test objects, a simulated stratigraphic phantom (a) and a heritage corroded copper sample obtained with classical transmission tomography (b) segmented in 6 classes.

  Copper sample.

Figure 4 . 3 :

 43 Figure 4.3: Respective V-line transforms times the factor cos(ω)σ(ω) uniformly sampled on the opening angle ω of objects in Figure 4.2. Notice how most of the energy of the diffracted beam is detected at very low angles and close to the incident beam position.

Figure 4 . 4 :Figure 4 . 5 :

 4445 Figure 4.4: Respective V-line transforms times the factor cos(ω)σ(ω) in log scale of objects in Figure 4.2. The log scale allows to appreciate the response of the different elements in the input object, such as the grains inside the stratigraphical sample making lines in the shape of the PSF represented in Figure 3.2.

Figure 4 . 6 :

 46 Figure 4.6: Reconstructions of the real phantom in 4.2. Same layout as in Figure 4.5.

Figure 4 . 7 :

 47 Figure 4.7: A slice of 16 µm (8 pixels) of width of both phantoms.

Figure 4 . 8 :

 48 Figure 4.8: Direct problem, image formation process. Output angular density (left column) representing the limits of each energy channel by a vertical red line, followed by the summation of projections in each 50 eV energy channel and then the estimation of the angular sampled density by nearest neighbours interpolation of energy data (right column).Notice the non-uniforme spacing of the angles bounding the energy channels and that a large fraction of the information is contained in the first energy channel.

Figure 4 . 9 :

 49 Figure 4.9: Interpolation by nearest neighbours of the scattered photon flux energy density from the finite energy resolved data. Different spectral resolutions exposed for two detecting sites where we can appreciate the errors caused by interpolations when the energy resolution is larger.

  Figure 4.8. A correction by the factor cos ω σ 2D (ω) must be performed to the estimate angular

Figure 4 . 10 :

 410 Figure 4.10: Filtering process related to the stratigraphical phantom. From the angular density (a) the 1D Fourier transform of this density is obtained in (b) and of the same density by replacing interpolated values of the second channel to the first one (c). The application of the ramp filter (d), followed by the cosine apodization window (e) and the filtered projections through the inverse Fourier transform (f). The spectral resolution is fixed to be ∆E = 50 eV. All images are in logarithmic scale.

Figure 4 . 11 :

 411 Figure 4.11: The amount of data collected as a function of the x-coordinate of a pixel.

  validates the optimal distance x 0 near 140 when more data is recorded. The same Figure shows that the variance ratio is increasing as a function of x 0 while the correlation coefficient presents a global maximum. The right No intersection in translations.

  Intersecting translations such as 1 bordering pixel line per slice is removed.

Figure 4 . 12 :

 412 Figure 4.12: Edge artifacts in reconstructions. Per each slice, bordering pixels suffers from a huge scale diminution (a), these pixels are taken out in (b).

  x 0 = 40 pix.

  x 0 = 80 pix.

  x 0 = 200 pix.

Figure 4 . 13 :

 413 Figure 4.13: Reconstructions placing the copper sample at different distances to the detector, resolved angularly with incomplete data. The effect of the hole missing data is appreciated following the different distances x 0 from phantom to detector.

  Stratigraphical sample, the slice considered in reconstructions at the white line is chosen to got through the two red grains. Copper sample.

Figure 4 .

 4 Figure 4.14: 1 pixel horizontal line profiles to be analysed.

Figure 4 . 15 :

 415 Figure 4.15: Reconstruction profile from a single slice of phantoms located at the white lines in Figure 4.13 according to different distances between phantom and detector, in pixels. Reconstructions from angular densities. Notice when the hole is present the main difference is the decay in scale, due to the lack of low angles data.

Figure 4 . 16 :

 416 Figure 4.16: Error measures related to the copper phantom as a function of the distance x 0 .

Figure 4 .

 4 Figure 4.17b where an important decay in scale for deeper regions of the object is presented. This effect is indeed reduced by placing the object further at the price of losing information corresponding to large angles.Figure4.18 shows reconstruction results of the same sample placed at x 0 = 50 pixels, that is, as in Figure4.17d with different spectral resolutions, namely 1, 5, 25 and 100 eV where we can appreciate the blurring in reconstructions affected by a larger spectral resolution. Actually, with a 5 eV resolution we dispose of more than 700 projections where the largest angular channel spread corresponds to 0.035 radians. This angular channel corresponds to the energy channel covering from E π to E π + ∆E, that is, the channel related to low angular values. With a resolution of 50 eV we have 74 projections and the largest angular channel is 0.171 radians wide while for 100 eV, we have 37 projections and a maximum angular channel of 0.242 radians width. Notice that reconstruction in 4.13 from angular resolved data were obtained with an angular step of 0.005 radians.An ideal detector resolved in angle is obviously not realistic and gives better results in terms of correlation coefficient as the angular sampling is uniform and then no numerical artifacts are generated due to interpolations from spectral images to energy densities. Without considering the Doppler broadening effect, reconstructions depend exclusively on the energy resolution of the detector and the quality of interpolations. Therefore, an energy resolution broader than 100 eV (corresponding to 37 non uniform angular projections) will not provide sufficient observations. As mentioned, nearest neighbour interpolation of the spectral image could be improved by applying a C 1 interpolation to avoid discontinuity effects in densities and better estimation of their derivatives.The farther the sample is placed, the less data related to big angles will be obtained. However, if the sample is placed closer to the detector, the effect of the hole will be bigger

  x 0 = 1 pix, ∆E = 50 eV.

  x 0 = 25 pix, ∆E = 50 eV.

  x 0 = 50 pix, ∆E = 50 eV.

  x 0 = 100 pix, ∆E = 50 eV.

  x 0 = 150 pix, ∆E = 50 eV.

Figure 4 . 17 :

 417 Figure 4.17: Reconstructions placing the copper sample at different distances to a detector with an energy resolution of 50 eV. The effect of the low angles missing data due to the presence of a hole is appreciated following the different distances x 0 from phantom to detector.

  x 0 = 50 pix, ∆E = 25 eV.

  x 0 = 50 pix, ∆E = 100 eV.

Figure 4 . 18 :

 418 Figure 4.18: Reconstructions placing the copper sample at 50 pixels from the detector with low angles incomplete data. Spectral resolutions of 1, 5, 25 and 100 eV are presented. The blurring in reconstructions is due to the missing data generating after performing summation of projections in each energy channel. Horizontal artifacts are due to the vertical translations of the object.

f

  (ζ + r sin ω cos ψ, ξ + r sin ω sin ψ, r cos ω) dψ dr r .

Figure 5 . 1 :

 51 Figure 5.1: Point spread function related to the conical Radon transform for an impulse function located at (0, 0, 1).

Fig. 1 2 M 2 M

 122 τ is the area of the detecting region D and r M the Euclidean distance from M to D. ω M is the scattering angle from M to D. As ω M ∈ [0, π], then, ωM is defined to be bounded in [0, π/2], namely ωM = min(ω M , π -ω M ). If τ is small enough, Ω(M, D) can be approximated by τ r cos ωM . Consequently, the scattered photon flux energy density at D, surrounding a scattering site M is given by dI(E, D|M) = τ I 0 1 r cos ωM σ(ω M , E) f (M) dM dE.

  Support of the Klein-Nishina DCS given by the Compton equation.

Figure 5 . 2 :

 52 Figure 5.2: The Klein-Nishina differential cross section (DCS)

Figure 5 . 3 :

 53 Figure 5.3: Normalized Compton scattering DDCS with respect to atomic number and incident energy. Two materials and two scattering angles are considered for E 0 = 50 keV.

2 f 2 f

 22 (ζ + tz cos ψ, ξ + tz sin ψ, z)e -2πi(uζ+vξ) dζdξ dψ dz z (5.29) With a change of variables translating the vertex of the cone to the origin, that is, with x = ζ + tz cos ψ and y = ξ + tz sin ψ, (x, y, z)e -2πi(ux+vy) dxdy e 2πitz(u cos ψ+v sin ψ) dψ dz z .

  we first estimate an energy density dI dE (E, D) by interpolating the spectral image into a thin enough uniform energy sampling.Then, through a change of variables given by the Compton equation, we obtain the required angular

  Aluminium, Z = 13.

Figure 6 . 1 :

 61 Figure 6.1: DDCS for hydrogen and aluminium for ω ∈ [ π 2 , π]. The spreading increases with the angle and is bigger for heavier elements. Notice that the maximum value per each ω is located at the Compton line, the spreading is due to the Doppler broadening. Discontinuities are due to the binding effect. E 0 = 50 keV.

Figure 6 .

 6 Figure 6.2 shows the real and fitted DDCS for aluminium at ω = 2 3 π. This was done through the R function nls.

Figure 6 . 2 : 3 π

 623 Figure 6.2: Compton scattering DDCS at ω = 2 3 π for aluminium and its approximation by a Lorentzian low with spreading factor γ = 495 eV which specifies the half-width at half-maximum. E 0 = 50 keV.

Figure 6 . 3 :

 63 Figure 6.3: Compton scattering DDCS at ω = 2 3 π for aluminium and its approximation by a Lorentzian low wit spreading factor γ = 495 eV which specifies the half-width at half-maximum. E 0 = 50 keV.

v 2 ( 2 + v 2 )

 222 a) The (u 2 + v 2 ) filter.(c) Cosine window with n = 2.

Figure 6 . 5 :

 65 Figure 6.5: Conical projections of the 3D Shepp-Logan phantom at a fixed ω.

Figure 6 . 6 :

 66 Figure 6.6: Conical filtered projections of the 3D Shepp-Logan phantom at a fixed ω. See how contrasting details are amplified with respect to the non filtered projections.

  : Phantom and reconstruction for z = 32.

  : Phantom and reconstruction for z = 40.

Figure 6 . 7 :

 67 Figure 6.7: 3D Shepp-Logan phantom of dimensions 64 × 64 × 64 voxels. Planes presented having equations z = z 0 with z 0 ∈ {24, 32, 40}. Left column, original phantom. Right column, reconstruction. Notice the ellipsoidal artifacts more accentuated when z is bigger due to less data from large angular values.

  : Phantom and reconstruction for y = 32.

  : Phantom and reconstruction for y = 40.

Figure 6 . 8 :

 68 Figure 6.8: 3D Shepp-Logan phantom of dimensions 64 × 64 × 64 voxels. Planes presented having equations y = y 0 with y 0 ∈ {24, 32, 40}. Left column, original phantom. Right column, reconstruction. Again, notice the influence of limited data related to large angular values as z increases.

  Mean energy E = 42541 eV. Mean energy E = 43291 eV.

Figure 6 . 9 :

 69 Figure 6.9: Spectral image for the analytical phantom at different energy channels. The free-electron case obtained by a conical integration of the object using the CPU parallelization implementation.

Figure 6 . 10 :

 610 Figure 6.10: Spectral image for the analytical phantom at different energy channels. These non-free electron simulations were obtained through a cartesian integration of the simulated sample using a GPU implementation.

  Mean energy E = eV.

  Mean energy E = 42541 eV.

  Mean energy E = 43291 eV.

Figure 6 . 11 :

 611 Figure 6.11: Spectral image for the analytical phantom at different energy channels. Results obtained using a complete Monte Carlo simulation performed with the code Penelope.

Figure 6 . 16 :

 616 Figure 6.16: The constructed phantom of dimensions 16×16×120 voxels. Planes presented, from left to right: z = {66, 98, 138} µm.

Figure 6 . 17 :Figure 6 . 18 :Figure 6 . 19 :

 617618619 Figure 6.17: Reconstruction from angular resolved data performed with half translations. Same layout as in Figure 6.16.

  The hole in the detector, its has no effect from the third in the free-electrons model.

  The hole is filled with numerical heat propagation.

Figure 6 . 20 :

 620 Figure 6.20: The hole in the detector is filled with heat propagation. Images presented for the three lowest energy channels, namely the 3 channels defined by the 4 energy boundaries (41816, 41866, 41916, 41966) eV.

Figure 6 .

 6 Figure 6.21: A slice of the phantom and its filtered back-projection reconstruction when the complete DDCS model is considered that accounts for Doppler and binding effects.

  cone.proj <-function(phantom, xyr.grid, omega, psi) { # CONICAL PROJECTIONS WITH LINEAR INTERPOLATION # inputs: # phantom: 3D array with attributes: xs, ys, zs # (coordinates of voxels middle points) # xyr.grid: grid (detector.x, detector.y, r) having colnames = c("xd","yd","r") # and attributes: x.size, y.size, r.size with each sampling sizes # and delta.x, delta.y, delta.r with each sampling steps # detector.x is the x sampling of the detector and # r the sampling of distances over the conical surface from the vertex # omega, psi: 2 real numbers, the polar and azimuthal angle # of a half-line in the cone # output # 2D array of dimensions (x.size, y.size) with projections for (omega, psi) # setup variables x.size <-attr(xyr.grid, 'x.size') y.size <-attr(xyr.grid, 'y.size') r.size <-attr(xyr.grid, 'r.size') xs <-attr(phantom, 'xs') ys <-attr(phantom, 'ys') zs <-attr(phantom, 'zs') # interpolating positions s1 <-sin(omega) xout <-xyr.grid$xd+xyr.grid$r*s1*cos(psi) yout <-xyr.grid$yd+xyr.grid$r*s1*sin(psi) zout <-xyr.grid$r*cos(omega) # taking only on the region of interest # x-coordinate index.x <-which(xout >= min(xs) & xout <= max(xs)) xout <-xout[index.x] yout <-yout[index.x] zout <-zout[index.x] # y-coordinated index.y <-which(yout >= min(ys) & yout <= max(ys)) xout <-xout[index.y] yout <-yout[index.y] zout <-zout[index.y] # linear interpolations result.1 <-approx3d(x = xs, y = ys, z = zs, f = phantom, xout = xout, yout = yout, zout = zout) # remove nan result.1[is.na(result.1)] <-0 result.2 <-0*index.x result.2[index.y] <-result.1 result <-rep(0, nrow(xyr.grid)) result[index.x] <-result.2 # rearranging result in an array (xd,yd,r) and * by 1/r result <-array(data = result/xyr.grid$r, dim = c(x.size, y.size, r.size)) # sum over r result <-rowSums(result, dim = 2) return(result) } Finally, the discrete conical Radon transform is obtained from last function as the 3D array (Cf (ζ m , ξ n , ω l )) m,n,l = ∆ψ∆r P -1 j=0 (P(ζ m , ξ n )) m,n (ω l , ψ j ), (A.2)

  cone.transform <-function(phantom, xyr.grid, omega.vector , psi.vector) { # CONICAL TRANSFORM # inputs: # phantom and xyr.grid as in the function cone.proj # omega.vector: uniform sampling of opening angles, attribute: delta # psi.vector: uniform sampling of azimuthal angles, attribute: delta # output: # 3D array of dimensions (x.size, y.size, omega.size) with attributes # delta.x, delta.y, omega.vector # setup variables x.size <-attr(xyr.grid, 'x.size') y.size <-attr(xyr.grid, 'y.size') delta.r <-attr(xyr.grid, 'delta.r') delta.psi <-attr(psi.vector, 'delta') omega.vector <-attr(omega.vector, 'delta') # initial array to save projections projections <-rep(0, x.size * y.size * omega.size) # for loop in azimuthal angle psi for (j in psi.vector) { # apply the function for all omega in parallel projections <-projections + mcmapply(FUN = function(x) cone.proj(phantom = phantom, xyr.grid = xyr.grid, omega = x, psi = j), omega.vector, mc.cores = detect.cores()) } # sampling steps projections <-projections * delta.psi * delta.r attributes(projections)$delta.x <-attr(xyr.grid, 'delta.x') attributes(projections)$delta.y <-attr(xyr.grid, 'delta.y') attributes(projections)$omega.vector <-omega.vector return(projections) }

  back.proj <-function(projections, x,y,z, psi.vector) { # CONICAL BACK-PROJECTIONS # input: # projections: 3D array of dimensions (x.size, y.size, omega.size) # with attributes: delta.x, delta.y, delta.omega # x,y,z: voxel coordinates # psi.vector: uniform sampling of azimuthal angles # output: # real number of the cumulate back-projections on (x,y,z) # setup variables dims <-dim(projections)

  result <-matrix(nrow = omega.size, ncol = psi.size) # 2D interpolation per omega for (i in 1:omega.size) { # object to interpolate obj <-list(x = detector.x, y = detector.y, z = projections[,,i]) # locations where the object is interpolated zeta <-x-z*tan(omega.vector[i])*cos(psi.vector) xi <-y-z*tan((omega.vector[i])*sin(psi.vector) loc <-cbind(zeta, xi) # interpolation result[i,] <-interp.surface(obj, loc) } # double omega -psi integration and factor z^2 result <-sum(result, na.rm = TRUE)*z^2 return(result) }

  two from R OpenCL module __global float* output, // Output const unsigned int count, // Size (in index) of the output array __global float* inds, // The various size, see later __global float* phantom, // Density of the volume, indexed by voxel __global float* energies, // 1D arrays of boundaries of E channels (contiguous) // should be ascending order __global float* doppler // array containing the 3 coefficients to //simulate doppler broadening per voxel ) { // Number of items on certain dimensions : size_t det_x, det_y, en_c, pht_x, pht_y, pht_z; det_x = get_global_size(0); // Copying the results to the output buffer : for (size_t i=0; en_c != i; ++i) { output[out_base_ind + i] = detector[i]; } }

  simulate.direct <-function(kernel, phantom, doppler, energies, E0, num.tiles, det.size, det.pixel.edge.size) { start.time <-Sys.time() cat("Starting date/time : ") print(start.time) ## Preparing the input to the OpenCL kernel p <-clFloat(phantom) d <-if (missing(doppler)) clFloat(c(0, 0, 0)) else clFloat(doppler) e <-clFloat(energies) counts <-prod(det.size) * (length(energies)-1) inds <-c(c(0, 0), length(energies)-1, dim(phantom), attr(phantom, 'entry.z'), det.pixel.edge.size, attr(phantom, 'edge.size'), E0) ## Preparing the data-object to store the results res <-array(NA, dim=c(num.tiles[1]*det.size[1], num.tiles[2]*det.size[2], (length(energies)-1))) ## Sampling the tiles : tile.sampling <-expand.grid(1:num.tiles[1], 1:num.tiles[2]) ## print(tile.sampling) ## Computing the min.min (bottom left) corner : bot.left = (num.tiles -1)/2 * det.size * -1 for (ti in 1:nrow(tile.sampling) ) { ## Setting properly the offset : inds[1:2] <-(c(tile.sampling[ti, 1], tile.sampling[ti, 2]) -1) * det.size + bot.left cat(sprintf("Computing tile %d (%d, %d), offset (%.2f, %.2f) ..", ti, tile.sampling[ti, 1], tile.sampling[ti, 2], inds[1], inds[2])) flush(stdout()) ## Running the OpenCL kernel res.tile <-oclRun(kernel, size=counts, inds, p, e, d, wait = TRUE, dim = det.size) res[(1:det.size[1]) + (det.size[1] * (tile.sampling[ti, 1]-1)), (1:det.size[2]) + (det.size[2] * (tile.sampling[ti, 2]-1)),] <-aperm(array(res.tile, c((length(energies)-1), det.size)), c(2, 3, 1)) cat((det.size[1] * (tile.sampling[ti, 1]-1)), (det.size[2] * (tile.sampling[ti, 2]-1))) cat(".. done\n") flush(stdout()) } attributes(res)$energies <-energies attributes(res)$E0 <-E0 attributes(res)$z0 <-attr(phantom, 'entry.z') attributes(res)$det.pixel.edge.size <-det.pixel.edge.size end.time <-Sys.time() cat("Ending date/time : ") print(end.time) cat("Total (wall-clock) executing time : ") print(end.time -start.time) return(res) }

/ 3 π

 3 / Heaviside step function per energy channel: float E_min, E_max; E_min = energies[0]; for (size_t E_ind = 0; en_c != E_ind; ++E_ind) { E_max = energies[E_ind+1]; detector[E_ind] += (compton_energy <= E_min || compton_energy > E_max) ? 0.0f : (vox_dens * solid_angle_factor * KN_factor); E_min = E_max; } 3.2 Point Spread Function related to the V-line transform for an impulse function located at (20, 0). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.3 A single scattering event from M detected at D ζ with a solid angle dΩ. . . . 48 4.1 Parameters and setup description of the 2D backward Compton scattering tomography. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.2 First test objects, a simulated stratigraphic phantom (a) and a heritage corroded copper sample obtained with classical transmission tomography (b) segmented in 6 classes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.3 Respective V-line transforms times the factor cos(ω)σ(ω) uniformly sampled on the opening angle ω of objects in Figure 4.2. Notice how most of the energy of the diffracted beam is detected at very low angles and close to the incident beam position. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.4 Respective V-line transforms times the factor cos(ω)σ(ω) in log scale of objects in Figure 4.2. The log scale allows to appreciate the response of the different elements in the input object, such as the grains inside the stratigraphical sample making lines in the shape of the PSF represented in Figure 3.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.5 Reconstructions of the simulated phantom . . . . . . . . . . . . . . . . . . . 65 4.6 Reconstructions of the real phantom . . . . . . . . . . . . . . . . . . . . . . 66 4.7 A slice of 16 µm (8 pixels) of width of both phantoms. . . . . . . . . . . . . 67 4.8 Direct problem, image formation process. Angular density, the related spectral image and the estimation of the angular density. . . . . . . . . . . . . . 68 4.9 Interpolation by nearest neighbours of the scattered photon flux energy density from the finite energy resolved data. Different spectral resolutions exposed for two detecting sites where we can appreciate the errors caused by interpolations when the energy resolution is larger. . . . . . . . . . . . . . . 69 4.10 Filtering process related to the stratigraphical phantom. . . . . . . . . . . . 70 4.11 The amount of data collected as a function of the x-coordinate of a pixel. . 72 4.12 Edge artifacts in reconstructions. Per each slice, bordering pixels suffers from a huge scale diminution (a), these pixels are taken out in (b). . . . . . 73 4.13 Reconstructions placing the copper sample at different distances to the detector, resolved angularly with incomplete data. The effect of the hole missing data is appreciated following the different distances x 0 from phantom to detector. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.14 1 pixel horizontal line profiles to be analysed. . . . . . . . . . . . . . . . . . 74 4.15 Reconstruction profile from a single slice of phantoms from angular projections. 74 4.16 Error measures related to the copper phantom as a function of the distance x 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.17 Reconstructions placing the copper sample at different distances to a detector energy resolved by 50 eV . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.18 Reconstructions placing the copper sample at 50 pixels from the detector with low angles incomplete data. Spectral resolutions of 1, 5, 25 and 100 eV are presented. The blurring in reconstructions is due to the missing data generating after performing summation of projections in each energy channel. Horizontal artifacts are due to the vertical translations of the object. 5.1 Point spread function related to the conical Radon transform for an impulse function located at (0, 0, 1). . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Klein-Nishina differential cross section (DCS) of Compton scattering. . . . . 5.3 Normalized Compton scattering DDCS with respect to atomic number and incident energy. Two materials and two scattering angles are considered for E 0 = 50 keV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Double differential cross-sections for hydrogen and aluminium . . . . . . . . 6.2 Compton scattering DDCS at ω = 2 3 π for aluminium and its approximation by a Lorentzian low with spreading factor γ = 495 eV which specifies the half-width at half-maximum. E 0 = 50 keV. . . . . . . . . . . . . . . . . . . 6.3 Compton scattering DDCS at ω = 2 for aluminium and its approximation by a Lorentzian low wit spreading factor γ = 495 eV which specifies the half-width at half-maximum. E 0 = 50 keV. . . . . . . . . . . . . . . . . . . 6.4 The filter (u 2 + v 2 ) and its smoothing by two cosine apodization windows. . 6.5 Conical projections of the 3D Shepp-Logan phantom at a fixed ω. . . . . . . 6.6 Conical filtered projections of the 3D Shepp-Logan phantom at a fixed ω. See how contrasting details are amplified with respect to the non filtered projections. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.7 3D Shepp-Logan phantom of dimensions 64 × 64 × 64 voxels. Planes presented having equations z = z 0 with z 0 ∈ {24, 32, 40}. Left column, original phantom. Right column, reconstruction. Notice the ellipsoidal artifacts more accentuated when z is bigger due to less data from large angular values.106 6.8 3D Shepp-Logan phantom of dimensions 64 × 64 × 64 voxels. Planes presented having equations y = y 0 with y 0 ∈ {24, 32, 40}. Left column, original phantom. Right column, reconstruction. Again, notice the influence of limited data related to large angular values as z increases. . . . . . . . . . . . . 6.9 Spectral image for the analytical phantom at different energy channels. The free-electron case obtained by a conical integration of the object using the CPU parallelization implementation. . . . . . . . . . . . . . . . . . . . . . . 6.10 Spectral image for the analytical phantom at different energy channels.

Table 2 .

 2 1: Simulations parameters for the Shepp Logan phantom reconstruction.

	Transform	rotation steps detector sampling correlation
	(p, ϕ)-Radon	314 values,	340 values,	0.99
	(u, v)-cartesian	ϕ ∈ [0, π[ 314 values,	p ∈] -169, 170] 1024 values,	0.92
	(ω, ζ)-half line	v ∈ [-5, 5] 314 values,	u ∈] -512, 512] 1024 values,	0.90

Table 4 .

 4 1: Setup parameters represented in Figure 4.1 for an ideal detector resolved angularly. 1 pixel represents 4 µm 2 .

	E 0	50 keV.
	dω	0.005 radians
	6 µm. ζ m 2 11 and 2 12 µm. ζ 0
	l	8 µm.
	τ	2 µm.

Table 4 .

 4 2: Error measurements between the object and reconstructions according to the different sampling parameters. The dimensions and sampling affect most the correlation coefficient and the related angle as some contrasting features are correctly recovered according to the sampling. The variation is scale is not really remarkable here due to most of the signal is captured at low angles.

	Phantom	L	M correlation angle variance ratio
	Stratigraphical 2 10 100	0.906	25.04	0.71
			200	0.911	24.36	0.70
		2 11 100	0.949	18.38	0.78
			200	0.952	17.82	0.78
	Copper	2 10 100	0.956	17.06	0.79
			200	0.963	15.63	0.77
		2 11 100	0.972	13.59	0.90
			200	0.980	11.47	0.87

Table 4 .

 4 3: Setup parameters represented in Figure 4.1 for a finite spectral resolved detector. 1 pixel represents 2 × 2 µm 2 .

	E 0	50 keV.
	dE 1, 5, 25, 50 and 100 eV.
	ζ 0 ζ m	10 µm. 2 11 µm.
	l	16 µm.
	τ	2 µm.

Table 4 .

 4 4: Error measurements between the object and reconstructions according to the different setup parameters. For larger spectral resolutions the variance ratio is affected mostly due to large angular interval corresponding to the first energy channel.

	∆E (eV.) x 0 (pixels) correlation variance ratio
	1	25	0.97	0.77
		50	0.97	0.84
	5	25	0.99	0.85
		50	0.99	0.87
	25	25	0.97	0.36
		50	0.97	0.33
	50	1	0.94	0.21
		25	0.95	0.19
		50	0.96	0.17
		100	0.96	0.15
		150	0.96	0.13
	100	25	0.93	0.09
		50	0.94	0.08

Table 6 .

 6 1: Error measurements related to different configurations.

	angular sampling/energy resolution correlation variance ratio
	∆ω = 0.001 radians		0.97	0.84
	∆E = 5 eV.		0.94	0.72
	∆E = 50 eV.		0.79	0.18
	4
	x
	-2
	60	80	100	120	140	160
			z
	0	5	10	15	20
	4
	x
	-2
	60	80	100	120	140	160
			z
	0	2	4	6	8

  une compréhension mathématique d'un système d'imagerie, le problème direct et 123 le problème inverse doivent être modélisés. Différents paramètres physiques décrivant le dispositif doivent être étudies comme par exemple les résolutions spatiales et spectrales, les dimensions ou les configurations initiales. Une bonne compréhension de ces paramètres à la fois dans le processus de formation d'image de sortie et dans la reconstruction du volume d'entrée est cruciale. Une configuration bidimensionnelle (2D) simplifiée est considérée dans un premier temps pour répondre aux cas où le faisceau incident et les photons diffusés sont limités à un plan 2D. Cette limitation est d'abord envisagée car les équations sont plus faciles à manipuler, les schémas numériques sont moins complexes et moins couteaux en temps de calcul, et aussi car l'angle azimutal d'un événement de diffusion Compton est uniformément distribué sur l'intervalle [0, 2π[. Nous présentons le problème direct et inverse lié à ce scénario 2D basé principalement sur la transformée de Radon sur des lignes en V et son inversion reposant sur le principe de la rétro-projection filtrée [21, 15], en considérant aussi l'hypothèse d'électrons libres et au repos. Cela signifie que l'objet d'entrée est uniquement représenté par la densité électronique à une position donnée et pour un même matériau. Des simulations numériques supportent ce modèle où différents paramètres physiques sont testés. Une interprétation des résultats selon ces paramètres est présentée pour être ensuite utilisée dans le cadre tridimensionnel (3D) complet.

  Le noyau de convolution induit par l'effet d'élargissement Doppler peut être estimé directement à partir de la relation entre les images spectrales simulées avec l'approximation de Dirac et de Lorentz pour la DSED. Ceci peut ensuite être suivi d'un processus de reconstruction itératif car le problème direct est bien compris et la simulation mise en place ici donne une bonne approximation des données mesurées et permet ainsi d'assurer une excellente similitude entre les observations et le modèle. Une autre approche consisterait à traiter le problème sous une formulation de moindre-carrée linéaire: une reconstruction à travers un pseudoinverse de la matrice associée du problème direct peut non seulement traiter le problème d'élargissement Doppler mais aussi l'incomplétude des données. Ces travaux sont en cours de réalisation.

Table A .

 A 1: The R packages and their application.

	data.table fast import of large dataset, used to import Penelope simulations data
	fields	spatial data tools, used in 2D interpolations and surfaces plotting
	oce	oceanographic data analysis, used in 3D interpolations
	OpenCL	Interface between R and OpenCL, used in computing forward and
		backward conical projections on a GPU parallel implementation
	parallel	used for CPU parallelization
	ptw	parametric time warping, used to perform zero padding on signals
	SynchWave	wavelet transform, used in shifting the Fourier transform to apply filters
		127

Introduction

indicating the same probability for a photon to be scattered on one or the other branch of

Backward Compton Scattering Tomography: 2D framework

Numerical simulations for the 2D framework

Remerciements

The conical back-projection can also be computed in a GPU parallel architecture. To do that, bilinear interpolations of the angular projections on every angular plane has to be performed to compute the values on the detector corresponding to such angle. That is, for a fixed ω, we create a GPU kernel to obtain the values g(x i -z k tan ω cos ψ p , y j -z k tan ω sin ψ p , ω), (A.5) for all coordinates (x i , y j , z k ) of a voxel in the object and for all samples ψ p of the azimuthal angle. This has to be done inside a for loop on ω.

The interests of performing GPU computations lies specially on the good parallel structure of the GPU where the price to pay is one more dimension in the cartesian integration with respect to the conical integration. Interpolations are not longer necessary and in the inner loop only the distance and angle between pixel and voxel are calculated.

The direct and inverse problems under the first schema, the CPU parallelization, take about 4 hours to be completed on 24 virtual 2.93 GHz CPU cores, providing an aggregated performance of 288 GFLOPS in single precision. In the other hand, using a 3GB GPU card with a theoretical performance de 2.87 TFLOPS, it is done in about 1 minute for the same imaging configuration. Abstract : Three-dimensional characterization of flat ancient material objects has remained a challenging activity to accomplish by conventional X-ray tomography methods due to their anisotropic morphology and flattened geometry. To overcome the limitations of such methodologies, an imaging modality based on Compton scattering is studied in this work. Classical X-ray tomography treats Compton scattering data as noise, while in Compton scattering tomography the conditions are set such that Compton data become the principal image agent. Under these conditions, we are able to avoid relative rotations between the sample and the imaging setup enabling the 3D imaging of flat samples. Mathematically this problem is addressed by means of the conical Radon transform. A model of the direct problem is presented where the output of the system is the spectral image obtained from an input object when illuminated by a parallel and monochromatic beam. The inverse problem is then addressed to estimate the 3D distribution of the electronic density of the input object from the spectral image. The feasibility of the proposed imaging system is supported by numerical simulations.

List of Figures

Universit é Paris-Saclay

Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France