8 February 2017 Three-dimensional imaging of flat natural and cultural heritage objects by a Compton scattering modality
Author Affiliations +
J. of Electronic Imaging, 26(1), 011026 (2017). doi:10.1117/1.JEI.26.1.011026
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 dense material 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.
Guerrero Prado, Nguyen, Dumas, and Cohen: Three-dimensional imaging of flat natural and cultural heritage objects by a Compton scattering modality





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.

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 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).


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 ω.

Fig. 3

Compton backscattering fundamentals. A photon of energy E0 is Compton scattered by an angle ω=πω¯ generating a new photon of energy Eω.


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 E0 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 E0 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×104 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, 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


where mec2=511  keV is the rest mass energy of the electron.

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 R3 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.


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) f:R2×R+R+.

Essentially, as shown in Fig. 4, an incident photon of energy E0, 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.

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.


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 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 imaging of the full object.


Image Formation


Conical Radon Transform

The conical Radon transform2,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


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


The factor 1/r comes from the integration measure on the cone rsinωdψdr and the approximation of the solid angle 1r2cosω, 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=rcosω 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:

  • I0: the incident photon flux time density just before the scattering event at M.

  • σ(ω): the Klein–Nishina differential cross-section9 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


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 1r2τ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 by9


where re is the classical electron radius. The factor 12re2 comes from the diametrical transversal area of an electron πre2 and the uniformity of the azimuthal angle 12π.

Consequently, the scattered photon flux density at D, given a scattering site M, reads


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).


Object Reconstruction


Bessel Function and Hankel Transform

We recall the Bessel function of the first kind of order 0, and the Hankel transform of a function f of order 0 noted, respectively, J0 and h0,





The inverse expression of the Hankel transform reads




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=ζ+tzcosψ and y=ξ+tzsinψ, Eq. (13) becomes


where we recognize the (x,y)-Fourier transform of f, then the last expression reads


where we can identify the Bessel function of the first kind of order 0 by switching to polar coordinates via u=qcosβ and v=qsinβ by


Note the use of notation Cfp^ and f^p to point out the use of polar coordinates.

The last integral turns out to be the Hankel transform of the function g:z1z2f^p(q,β,z), then Eq. (16) reads



One can, therefore, apply the Hankel inverse equation to have


Then going back to f^p from g, we have


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




Adjoint Transform

Let us introduce some basic notations. C is defined as an operator from XL2(R2×R+) to YL2(R2×[0,π2]), and the respective inner products of those spaces are defined by


The adjoint transform C of C, closely related to the backprojection operator, is defined as the transform from Y to X that verifies


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


and apply a change of variables in the form x=ζ+tzcosψ, y=ξ+tzsinψ, and Fubini’s theorem to have



From which one is able to extract the adjoint transform C of C from identity (23) in the form


The way it is defined, the adjoint transform can be interpreted here as a backprojection procedure of projections g(xztanωcosψ,yztanω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




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 CCf(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 CCf(x,y,z) in accordance to such an inverse equation.

Therefore, to compute CCf(x,y,z) from Eq. (23), we insert the bidimensional (ζ,ξ)-Fourier transform of C 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 z3(u2+v2)sinωcos3ω. 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


and, finally, we are able to write the inverse equation of the CRT in the form




Simulation Results

A numerical object corresponding to a flattened stratigraphic sample is created, reconstructed by means of inversion equation (31), and presented in Figs. 5Fig. 6Fig. 78 for different cross sections of a phantom detailed as follows. The phantom is built with three layers of 512×512  μm2 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±0.1, 6.0±0.5, and 2.0±0.3 for each layer, respectively.

Fig. 5

Stratigraphical phantom of dimensions 128×128×32 voxels. Plane presented: z=4. (a) original phantom, (b) direct inversion, (c) inversion after normalizing values following z-axis, and (d) inversion after removing voxels located at edges of single reconstructions of dimensions 4×4×32 voxels.


Fig. 6

Stratigraphical phantom of dimensions 128×128×32 voxels. Plane presented: z=30. Same layout as in Fig. 5.


Fig. 7

Stratigraphical phantom of dimensions 128×128×32 voxels. Plane presented: y=22. From top to bottom: original phantom, direct inversion, inversion after normalizing values following z-axis, inversion after removing voxels located at edges of single reconstructions of dimensions 4×4×32 voxels.


Fig. 8

Stratigraphical phantom of dimensions 128×128×32 voxels. Plane presented: y=25. Same layout as in Fig. 7.



Discretization Parameters

The unit length considered in discretizations is 4  μm thus, each pixel in the images represents 16  μm2.

A square parallel x-ray beam of 256  μm2 of sectional area crosses the hole in the detector of 576  μm2 (ζ0=ξ0=12  μm). The medium has dimensions of 512×512×128  μm3 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  μm2 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.


Numerical Reconstruction

To numerically compute inversion Eq. (31) from data Cf(ζ,ξ,ω) generated by Eq. (3), the procedure works as follows:

  • 1. Filter the projections in the Fourier domain to get C*f^(ζ,ξ,ω) applying Eq. (30).

  • 2. Perform an inverse Fourier transform to obtain the filtered projections C*f(ζ,ξ,ω).

  • 3. The backprojection procedure is performed by means of the adjoint transform C by Eq. (26).

  • 4. Finally, applying the factor z3, 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 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


M. Haltmeier, “Exact reconstruction formulas for a radon transform over cones,” Inverse Probl. 30(3), 035001 (2014).INPEEY0266-5611http://dx.doi.org/10.1088/0266-5611/30/3/035001Google 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


S. J. Norton, “Compton scattering tomography,” J. Appl. Phys. 76, 2007–2015 (1994).JAPIAU0021-8979http://dx.doi.org/10.1063/1.357668Google Scholar


M. Morvidone et al., “On the V-line radon transform and its imaging applications,” Int. J. Biomed. Imaging 2010, 11 (2010).http://dx.doi.org/10.1155/2010/208179Google Scholar


T. T. Truong and M. K. Nguyen, “On new V-line radon transforms in R2 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


O. Klein and Y. Nishina, “The scattering of light by free electrons according to Dirac’s new relativistic dynamics,” Nature 122, 398–399 (1928).http://dx.doi.org/10.1038/122398b0Google 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.

Patricio Guerrero Prado, Mai K. Nguyen, Laurent Dumas, Serge X. Cohen, "Three-dimensional imaging of flat natural and cultural heritage objects by a Compton scattering modality," Journal of Electronic Imaging 26(1), 011026 (8 February 2017). http://dx.doi.org/10.1117/1.JEI.26.1.011026
Submission: Received 1 July 2016; Accepted 29 November 2016

Compton scattering



3D image processing

Cultural heritage

Image processing

Image acquisition

Back to Top