Three-dimensional imaging of flat natural and cultural heritage objects by a Compton scattering modality

Characterization and interpretation of flat ancient material objects, such as those found in archaeology, paleoenvironments, paleontology, and cultural heritage, have remained a challenging task to perform by means of conventional x-ray tomography methods due to their anisotropic morphology and flattened geometry. To overcome the limitations of the mentioned methodologies for such samples, an imaging modality based on Compton scattering is proposed in this work. Classical x-ray tomography treats Compton scattering data as noise in the image formation process, while in Compton scattering tomography the conditions are set such that Compton data become the principal image contrasting agent. Under these conditions, we are able, first, to avoid relative rotations between the sample and the imaging setup, and second, to obtain three-dimensional data even when the object is supported by a densematerial by exploiting backscattered photons. Mathematically this problem is addressed by means of a conical Radon transform and its inversion. The image formation process and object reconstruction model are presented. The feasibility of this methodology is supported by numerical simulations. © The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI. [DOI: 10.1117/1.JEI .26.1.011026]


Context
Material characterization of ancient flat object encountered in natural or cultural heritage studies remains notably challenging with nowadays x-ray imaging methods.When one deals with heritage objects, the noninvasiveness and nondestructiveness properties of inspections are a requirement that x-ray imaging methodologies provide, enabling two-dimensional (2-D) imaging of the sample.However, one can easily be facing samples presenting a flattened geometry, i.e., samples presenting a large ratio between its front area and thickness.The challenge is then to perform a three-dimensional (3-D) probing without using a relative rotation between the sample and the imaging setup as 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, e.g., studies of conservation/restoration of easel paintings requiring 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 study is shown in Fig. 1 which was performed by means of an invasive method.Figure 2 presents another example of objects possessing this morphology, namely in paleontology with the Lagerstätten fossils 1 which are mechanically flattened during the fossilization process and stand on one side of a thick sedimental slab which cannot 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-rays either due to its density or its thickness.
Two alternatives are currently 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 bidimensional analysis limited to the front surface of the sample, for example with synchrotron x-ray fluorescence spectral raster-scanning as performed in Ref. 1, none of which enabling a 3-D study of the sample.Furthermore, because of the opaque supporting material, transmission, and forward scattering data, i.e., with a scattering angle inferior to π∕2, is impossible to collect.This motivates the proposal of a modality using backscattered data, i.e., data collected with a scattering angle ω comprised between π∕2 and π, as shown in Fig. 3.In the following, we will call ω ¼ π − ω the supplementary angle of ω.
Identity (1) introduces the Compton equation, a diffeomorphism between the scattering angle and the scattered energy under the hypotheses that electrons are both free and at rest and that the incident beam is monochromatic, i.e., a single value E 0 of an incident energy irradiating the object.[4] 1.2 Compton Scattering Tomography For the energy range from 50 eV up to 1 MeV, the dominant photon-matter interaction processes are the photoelectric absorption, Rayleigh scattering which is both elastic and coherent, and Compton scattering which conversely is both inelastic and incoherent.In this last one, an incident photon of energy E 0 is absorbed by a target electron, which re-emits a secondary photon scattered by an angle ω relative to the direction of the original photon.The scattered photon has then an energy E ω which is related to the scattering angle ω by Eq. ( 1), the Compton equation, presented later.
In classical x-ray imaging and tomography, the Compton scattered signal is considered as noise added to photoelectric absorption and coherent scattered data.This is because the x-ray transmission signal is dominated by the photoelectric absorption while coherent scattering may produce significant amplitude variations at low scattering angles because of constructive and destructive interference effects due to the coherent nature of this scattering.However, depending on the material, if the incident radiation has a superior 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).
Classical tomographic imaging modalities developed and used in most all applications in the last half century include transmission computed tomography, single photon emission tomography, and positron emission tomography.All of them regard primary radiation and perform 3-D mappings leaning on relative rotations between the object and the imaging setup.In such a framework, Compton scattering events are adding a nonuniform background to the observation, a systematic bias which leads to artifacts in the reconstruction if it is not accounted for.As the relative importance of the Compton scattering effect over the other two processes mentioned above is increasing with an increase on incident energy, its effect is even more important when using higher energy γ-ray imaging.
The idea of exploiting scattered radiation from the Compton effect in imaging techniques has been introduced and studied simultaneously and it has given birth to CST, 2,4-8 which focuses on reconstructing the electron density map of the object.
To compensate the information provided by multiple projection angles in classical tomography, CST exploits the energy loss of the scattered photons.This energy loss is related to the scattering angle via the Compton equation given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 3 0 1 where m e c 2 ¼ 511 keV is the rest mass energy of the electron.
The first CST scanner was proposed in 1994 5 through a Radon transform over an arc of circles starting at a γ-ray source and ending at a fixed detecting site.Radon transforms over conical surfaces having fixed axis directions and variable opening angle are studied in Ref. 2 where a first analytic inversion equation is proposed using circular component analysis, with applications to emission imaging based on Compton scattered radiation.Generalizations of this kind of transform to higher dimensions spaces with related inverse formulas in a filtered backprojection type are presented in Ref. 3. A backprojection inversion algorithm for a conical Radon transform in R 3 was recently developed in Ref. 4.  This paper is organized as follows.Section 2 describes the imaging configuration.Section 3 defines the conical Radon transform and presents the image formation process.   2 Proposed Setup CST aims to reconstruct the electron density map of the studied object.In this work, the electron density is represented mathematically by a nonnegative bump function (both smooth and of bounded support) Essentially, as shown in Fig. 4, an incident photon of energy E 0 , belonging to a monochromatic parallel beam with a square section, is Compton scattered by an electron situated inside the object at M subtending an angle ω ( π 2 < ω < π) with the direction of incidence.The scattered photon, of energy E ω approximated by Eq. ( 1), reaches a detecting site D ¼ ðζ; ξ; 0Þ on the 2-D detector that is located over the xOy-plane.The parallel beam is centered at the Oz-axis.
The imaging configuration is described then as follows: as mentioned, a synchrotron radiation setup with a parallel monochromatic x-ray beam (about 50 keV) and a spaceenergy resolved detector are considered.As represented in Fig. 4, the detector will be placed between the source and the object to capture backscattered 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 jζj > ζ 0 or jξj > ξ 0 .Horizontal and vertical translations of the sample will be needed to allow imaging of the full object.
3 Image Formation

Conical Radon Transform
The conical Radon transform 2,3 integrates a function over circular cone surfaces; its backprojection reconstruction procedure is developed in Ref. 4.
A ðψ; rÞ-parametrization of the lateral surface of a circular cone is considered, where φ is the azimuthal angle and r is the distance of a point in the cone surface from its vertex.Such parametrization reads, for a cone with vertex at D ¼ ðζ; ξ; 0Þ and opening angle ω, noted C ω;D E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 7 1 9 C ω;D ¼ ðζ þ r sin ω cos ψ; ξ þ r sin ω sin ψ; r cos ωÞ: (2) Therefore, the CRT applied on an electron density function f, noted Cf, defined as the integral of this function over the surface of a cone parametrized as Eq. ( 2) reads E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 6 4 2 The factor 1∕r comes from the integration measure on the cone r sin ωdψdr and the approximation of the solid angle 1 r 2 cos ω, the factor sin ω cos ω not appearing in the last integral will be taken into account later in the scattered photon flux time density in Eq. ( 9).
It may be useful to express Eq. ( 3) as an integral with respect to z ¼ r cos ω and introducing the variable t ¼ tan ω.If notation Cf is conserved, Eq. ( 3) is also expressed as

Scattered Photon Flux Time Density
Let IðE ω ; DÞ represent the recorded scattered photon flux density (number of photons of energy E ω recorded per unit time at D).It incorporates the following parameters: • I 0 : the incident photon flux time density just before the scattering event at M. • σðωÞ: the Klein-Nishina differential cross-section 9 at an angle ω.
• fðMÞ: the electron density at M.
• dΩðM; DÞ: the solid angle from M to D.
• dM: the area element around M over the cone surface.
The solid angle dΩðM; DÞ can be seen from Fig. 4 to be E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 3 2 6 ; 2 3 2 dΩðM; where τ is the area of the detecting element located at D and r is the Euclidean distance from M to D.
If τ is small enough, then dΩðM; DÞ can be approximated by 1 r 2 τ cos ω.The Klein-Nishina differential cross-section is a function of 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πÞ.It is given by 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 3 2 6 ; 9 7 σðωÞ ¼ where r e is the classical electron radius.The factor 1 2 r 2 e comes from the diametrical transversal area of an electron πr 2 e and the uniformity of the azimuthal angle 1 2π .Consequently, the scattered photon flux density at D, given a scattering site M, reads E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 6 3 ; 6 9 4 dIðE ω ; DjMÞ ¼ I 0 σðωÞdΩðM; DÞfðMÞdM: The total scattered flux time density IðE ω ; DÞ recorded at D is the integral over all scattering sites laying over the surface of the cone C ω;D .It is hence given by the integral From the last integral, we can extract the conical Radon transform, and we are able to express the scattered photon flux density for a recorded energy E as Note the use of ω E and ωE , the scattering angle and its supplementary, related to the recorded energy E approximated by an inversion of identity (1).
The inverse expression of the Hankel transform reads E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 6 3 ; 3 1 8 fðrÞ ¼

Inverse Formula
We start writing the bidimensional ðζ; ξÞ-Fourier transform of Eq. ( 4) as and with a change of variables (translation of the vertex to the origin) x ¼ ζ þ tz cos ψ and y ¼ ξ þ tz sin ψ, Eq. ( 13) becomes where we can identify the Bessel function of the first kind of order 0 by switching to polar coordinates via u ¼ q cos β and v ¼ q sin β by Finally take the inverse Fourier transform in polar coordinates to recover f as The last integral may be expressed in Cartesian coordinates related to the Fourier domain ðu; vÞ and in terms of ω, then we are able to write the inverse equation of Eq. ( 3) as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 1 ; 3 2 6 ; 2 9 3 fðx; y; zÞ

Adjoint Transform
Let us introduce some basic notations.C is defined as an operator from , and the respective inner products of those spaces are defined by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 3 2 6 ; 1 6 5 The adjoint transform C † of C, closely related to the backprojection operator, is defined as the transform from Y to X that verifies E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 3 ; 6 3 ; 7 5 2 ðg; CfÞ Y ¼ ðC † g; fÞ X : In order to derive an expression for C † , we start from the left side of the last identity and introduce the definition of C given in Eq. ( 3) writing E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 4 ; 6 3 ; 6 9 7 ðg; and apply a change of variables in the form From which one is able to extract the adjoint transform C † of C from identity (23) in the form E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 6 ; 6 3 ; 4 7 8 C † gðx; y; zÞ ¼ The way it is defined, the adjoint transform can be interpreted here as a backprojection procedure of projections gðx − z tan ω cos ψ; y − z tan ω sin ψ; ωÞ to the position ðx; y; zÞ times the factor 1∕z, i.e., one assigns to ðx; y; zÞ the values of projections starting at this point, forming a cone toward the detector with opening angle ω times the mentioned factor.For a first rough reconstruction, fðx; y; zÞ can be approximated through the adjoint transform or equivalently in this case, through a simple backprojection of data as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 7 ; 6 3 ; 3 0 7 fðx; y; zÞ ≈ C † Cfðx; y; zÞ: (27)

From the Adjoint Transform to Filtered Backprojections
To show that the adjoint transform, C † is not the inverse operator of C, one only needs to compute C † Cfðx; y; zÞ and compare the result with the inversion equation of C Eq. ( 21) given in Sec.4.2.Then filters must be added to projections in order to adjust the result of C † Cfðx; y; zÞ in accordance to such an inverse equation.
Therefore, to compute C † Cfðx; y; zÞ from Eq. ( 23), we insert the bidimensional ðζ; ξÞ-Fourier transform of C to have E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 8 ; 6 3 ; 1 6 2 e −2πiz tan ωðu cos ψþv sin ψÞ dψdudvdω: (28) Last inner ψ-integration can be performed via polar coordinates as in Eqs. ( 15) and ( 16) allowing us to write the Bessel function of the first kind of order 0 as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 9 ; 3 2 6 ; 7 4 1 The last expression differs from inverse equation (21) only by factors z 3 ðu 2 þ v 2 Þ sin ω cos 3 ω .These factors are seen as filters that we need to apply onto the data to have an exact backprojection inversion procedure.
Therefore, let us define the projections C Ã fðζ; ξ; ωÞ filtered in the Fourier space by the filters E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 0 ; 3 2 6 ; 6 3 1 and, finally, we are able to write the inverse equation of the CRT in the form ; t e m p : i n t r a l i n k -; e 0 3 1 ; 3 2 6 ; 5 7 0 fðx; y; zÞ ¼ z 3 C † C Ã fðx; y; zÞ: (31)

Simulation Results
A numerical object corresponding to a flattened stratigraphic sample is created, reconstructed by means of inversion equation (31), and presented in Figs.5-8 for different cross sections of a phantom detailed as follows.The phantom is built with three layers of 512 × 512 μm 2 regarding its frontal area and about 40 μm of thickness per layer, with electron densities of 0.9, 1.1, and 1.0, respectively.Some grains of random diameters and positions are inserted inside the layers, having densities of 1.3 AE 0.1, 6.0 AE 0.5, and 2.0 AE 0.3 for each layer, respectively.

Discretization Parameters
The unit length considered in discretizations is 4 μm thus, each pixel in the images represents 16 μm 2 .A square parallel x-ray beam of 256 μm 2 of sectional area crosses the hole in the detector of 576 μm 2 (ζ 0 ¼ ξ 0 ¼ 12 μm).The medium has dimensions of 512 × 512 × 128 μm 3 thus we need 32 translations following the 0x-axis and the same for the Oy-axis, to cover the full medium.The detector plane, located over the xOy-plane at a perpendicular distance of 4 μm to the sample, has an area of 256 × 256 μm 2 including the hole.Regarding the number of detecting sites, having an area of 1 pixel each, and considering the absence of detecting sites in the hole, we have 4060 detecting sites.

Numerical Image Formation
Numerical resolution of integral ( 3) is performed via the trapezoidal rule with a discretization spatial step dr ¼ 4 μm and an azimuthal angular step dψ ¼ 0.01 rad.64 opening angles of the cone from 0 to π∕2 are considered.
3. The backprojection procedure is performed by means of the adjoint transform C † by Eq. ( 26). 4. Finally, applying the factor z 3 , we compute the inversion Eq. (31).

Inversion Procedure Details
Direct application of Eq. ( 31) is first presented in figures, just after the original sample, and one can see edge artifacts due to horizontal and vertical translations of the sample and the different scale rates between voxels located at the edge of a single reconstruction and those located at the center.First, the obvious solution presented is to normalize scales following the z-axis direction, where we can still see these edge artifacts, and conclude that problem is not only about scaling.In the last solution presented, we removed these voxels located at edges, where we can appreciate an encouraging 3-D reconstruction of the sample.A cosine window function was used in the Fourier domain of projections to control high frequencies.Data not recorded in the hole area were filled by a numerical solution of a heat diffusion problem from values of energies captured at the hole edge.

Conclusion and Perspectives
An x-ray imaging modality based on Compton scattering is presented; the interest is to perform a 3-D mapping of flat heritage objects without relying on relative rotations among object, source, and detectors.Both image formation by means of a conical Radon transform and object reconstruction by a filtered backprojection procedure are exposed and supported with numerical simulations showing feasibility with real data.Results are already very encouraging considering the problem of nondestructive and noninvasive 3-D 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 how it is coupled to the volume reconstruction itself.This is proving instrumental in designing the physical instrument that will ultimately enable such imaging modality.
The current method assumes an exact relationship between the diffusion energy E ω and the scattering angle ω according to Eq. ( 1).Through the Doppler effect and electron binding energy variations this relationship is indeed a non-Dirac probabilistic law.We will use a combination of Monte-Carlo simulations and variational analysis to assess the consequences of our approximation on image formation and volume reconstruction in a current work.
To be able to tackle the system computationally, we have first focused on flat objects and using incident energy that optimize the relative importance of single event inelastic scattering 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 versus 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, i.e., the effect of photoelectric absorption is much smaller than that of 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 the relative importance of absorption and scattering.
Still we plan to work on an iterative reconstruction process which, while more expensive computationally, will enable us to also account for the smaller effects, photoelectric absorption and multiple scattering, in the produced reconstruction.
the image of the Rembrandt painting, and the field workers who collected the fossil actinopterygian.Patricio Guerrero would like to thank the French Fondation des Sciences du Patrimoine (PATRIMA LabEx) for providing a PhD grant to support his research work.The same author would also like to thank Javier Cebeiro for his crucial help and discussions.
Patricio Guerrero Prado has been a PhD candidate in computed tomography, Radon transforms, and applied mathematics since November 2014 at the European research platform on ancient materials IPANEMA under the direction of Serge Cohen, Mai Nguyen, and Laurent Dumas.His thesis title is "Compton 2D imaging for 3D reconstruction."Mai K. Nguyen received her PhD in signal and image processing from Grenoble National Polytechnic Institute, France, in 1988.She has been a professor in the Department of Computer Sciences at the University of Cergy-Pontoise since 2005.Her research interests include inverse problems, generalized Radon transforms, and their applications in imaging science, scattered ionizing radiation imaging, biomedical imaging, nondestructive evaluation.She is an IEEE senior member.
Laurent Dumas received his PhD in applied mathematics from Paris Diderot University in 1995.He has been a professor at Versailles University since 2010.His research interests include numerical analysis, derivative free optimization, uncertainty quantification, and have been applied in various engineering and medical fields.
Serge X. Cohen earned his PhD in biophysics at EMBL Grenoble.After a postdoctorate period at NKI (Amsterdam) working on shape recognition and statistical learning applied to structural biology, he is now a mathematician at CNRS.Since 2010, leading information extraction and analysis scientists at IPANEMA, he has focused on theoretical statistics aspects for spectromicroscopy and synchrotron radiation tomography for ancient materials.

Fig. 1
Fig. 1 Paint cross-section showing a stratigraphical assemblage of The Anatomy Lesson of Dr. Nicolaes Tulp, 1632 by Rembrandt, Mauritshuis, The Hague.© Sample taken and prepared as cross-section by P. Noble during the conservation treatment of the painting in 1997 and rephotographed by A. van Loon, Mauritshuis, in 2010 for the Rembrandt Database.

Fig. 2 A
Fig. 2 A flat fossil actinopterygian over a thick, highly x-ray absorbing support from the Kem Kem Beds in Morocco dated back to the Lower Cretaceous (95 million years ago).© P. Gueriau (MHNM/MNHN).

Fig. 3
Fig. 3 Compton backscattering fundamentals.A photon of energy E 0 is Compton scattered by an angle ω ¼ π − ω generating a new photon of energy E ω .

Fig. 4
Fig. 4 Imaging configuration proposed.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.

1
Bessel Function and Hankel TransformWe recall the Bessel function of the first kind of order 0, and the Hankel transform of a function f of order 0 noted, respectively, J 0 and h 0 , E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 6 3 ; 4 1 2 J 0 ðxÞ ¼ T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 6 3 ; 3 6 7 h 0 fðkÞ ¼ Z ∞ 0 fðrÞJ 0 ðkrÞrdr:

5 c 2 e
E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 6 3 ; 1 4 Cfðu; v; tÞ ¼ Z ∞ fðx; y; zÞe −2πiðuxþvyÞ dxdy × e 2πitzðu cos ψþv sin ψÞ dψdz;(14)where we recognize the ðx; yÞ-Fourier transform of f, then the last expression reads E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2πitzðu cos ψþv sin ψÞ dψdz;(15)

Ed
Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 Cf p ðq;β;tÞJ 0 ð2πqtzÞtdtq 3 × e 2πiqðx cos βþy sin βÞ dqdβ: and Fubini's theorem to have E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 tan ω cos ψ;y − z tan ω sin ψ;ωÞdψdωdxdydz:(25)