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 fossils1 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 , is impossible to collect. This motivates the proposal of a modality using backscattered data, i.e., data collected with a scattering angle comprised between 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 of an incident energy irradiating the object. Equation (1) provides a reliable approximation of a real scattering scenario and with it, the scattering image formation process corresponds to a Radon transform over the surface of a cone with a fixed axis direction.23.–4
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 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 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 electron-volts, Compton scattering becomes the dominant phenomenon in the process, even more when detection is performed outside the direct transmission area ().
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,45.6.7.–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
The first CST scanner was proposed in 19945 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 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. Section 4 details object reconstructions via a backprojection inversion from Ref. 4. In Sec. 5, we present a numerical scheme and simulation results for flat objects both for energy resolved image formation and sample reconstruction. Finally, Sec. 6 closes this work with conclusions and perspectives.
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 , belonging to a monochromatic parallel beam with a square section, is Compton scattered by an electron situated inside the object at subtending an angle () with the direction of incidence. The scattered photon, of energy approximated by Eq. (1), reaches a detecting site on the 2-D detector that is located over the -plane. The parallel beam is centered at the -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 space-energy 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 for some two positive real numbers to allow the beam to go through. Therefore, we will have data for values of on the -plane verifying or . Horizontal and vertical translations of the sample will be needed to allow imaging of the full object.
Conical Radon Transform
A -parametrization of the lateral surface of a circular cone is considered, where is the azimuthal angle and is the distance of a point in the cone surface from its vertex. Such parametrization reads, for a cone with vertex at and opening angle , noted
It may be useful to express Eq. (3) as an integral with respect to and introducing the variable . If notation is conserved, Eq. (3) is also expressed as
Scattered Photon Flux Time Density
Let represent the recorded scattered photon flux density (number of photons of energy recorded per unit time at ). It incorporates the following parameters:
• : the incident photon flux time density just before the scattering event at .
• : the Klein–Nishina differential cross-section9 at an angle .
• : the electron density at .
• : the solid angle from to .
• : the area element around over the cone surface.
The solid angle can be seen from Fig. 4 to be
If is small enough, then can be approximated by .
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 . It is given by9
Consequently, the scattered photon flux density at , given a scattering site , reads
Bessel Function and Hankel Transform
We recall the Bessel function of the first kind of order 0, and the Hankel transform of a function of order 0 noted, respectively, and ,
The inverse expression of the Hankel transform reads
We start writing the bidimensional -Fourier transform of Eq. (4) as
The last integral turns out to be the Hankel transform of the function , then Eq. (16) reads
One can, therefore, apply the Hankel inverse equation to have
The last integral may be expressed in Cartesian coordinates related to the Fourier domain and in terms of , then we are able to write the inverse equation of Eq. (3) as
Let us introduce some basic notations. is defined as an operator from to , and the respective inner products of those spaces are defined by
From which one is able to extract the adjoint transform of from identity (23) in the form
For a first rough reconstruction, can be approximated through the adjoint transform or equivalently in this case, through a simple backprojection of data as
From the Adjoint Transform to Filtered Backprojections
To show that the adjoint transform, is not the inverse operator of , one only needs to compute and compare the result with the inversion equation of Eq. (21) given in Sec. 4.2. Then filters must be added to projections in order to adjust the result of in accordance to such an inverse equation.
Therefore, to compute from Eq. (23), we insert the bidimensional -Fourier transform of to have
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
The last expression differs from inverse equation (21) only by factors . 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 filtered in the Fourier space by the filters
A numerical object corresponding to a flattened stratigraphic sample is created, reconstructed by means of inversion equation (31), and presented in Figs. 5Fig. 6Fig. 7–8 for different cross sections of a phantom detailed as follows. The phantom is built with three layers of regarding its frontal area and about 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 , , and for each layer, respectively.
The unit length considered in discretizations is thus, each pixel in the images represents .
A square parallel x-ray beam of of sectional area crosses the hole in the detector of (). The medium has dimensions of thus we need 32 translations following the -axis and the same for the -axis, to cover the full medium. The detector plane, located over the -plane at a perpendicular distance of to the sample, has an area of 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 and an azimuthal angular step rad. 64 opening angles of the cone from 0 to are considered.
To numerically compute inversion Eq. (31) from data generated by Eq. (3), the procedure works as follows:
1. Filter the projections in the Fourier domain to get applying Eq. (30).
2. Perform an inverse Fourier transform to obtain the filtered projections .
3. The backprojection procedure is performed by means of the adjoint transform by Eq. (26).
4. Finally, applying the factor , 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 -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 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 authors would like to thank Pierre Gueriau and Mathieu Thoury for their support, discussions, and help to acquire images of heritage objects. Annelies van Loon for providing 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.
P. Gueriau et al., “Trace elemental imaging of rare earth elements discriminates tissues at microscale in flat fossils,” PLoS One 9(1), e86946 (2014).POLNCL1932-6203http://dx.doi.org/10.1371/journal.pone.0086946Google Scholar
M. K. Nguyen, T. T. Truong and P. Grangeat, “Radon transforms on a class of cones with fixed axis direction,” J. Phys. A: Math. Gen. 38(37), 8003–8015 (2005).JPHAC50305-4470http://dx.doi.org/10.1088/0305-4470/38/37/006Google Scholar
J. Cebeiro, M. Morvidone and M. K. Nguyen, “Back-projection inversion of a conical Radon transform,” Inverse Probl. Sci. Eng. 24(2), 328–352 (2015).http://dx.doi.org/10.1080/17415977.2015.1034121Google Scholar
T. T. Truong and M. K. Nguyen, “On new V-line radon transforms in and their inversion,” J. Phys. A: Math. Theor. 44(7), 075206 (2011).JPAMB51751-8113http://dx.doi.org/10.1088/1751-8113/44/7/075206Google Scholar
T. T. Truong and M. K. Nguyen, “New properties of the V-line radon transform and their imaging applications,” J. Phys. A: Math. Theor. 48(40), 405204 (2015).JPAMB51751-8113http://dx.doi.org/10.1088/1751-8113/48/40/405204Google Scholar
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.