Recent Developments on Compton Scatter Tomography: Theory and Numerical Simulations

1. Introduction

Progress in nuclear physics acquired during World War II has naturally led many scientists to devote their research activities to the field of tomographic imaging techniques using ionizing radiation. This was, in the post-war booming economy, particularly of great importance for medical diagnostic as well as for industrial non destructive evaluation (NDE). As the objective is to extract information on the inner part of objects of interest, penetrating radiation was the most appropriate agent for this purpose. With the availability of high quality X-ray and gamma-ray sources (either directly originated from nuclear transitions or from pair annihilation) and the emergence of sensitive detectors, three types of imaging have been introduced and developed throughout half a century. Nowadays they have emerged as mature standard investigation methods for several domains of application.

These are: a) the X-ray transmission Computed Tomography (CT), which exploits the physical law of radiation attenuation in matter, b) the Single Photon Emission Tomography (SPECT), which uses the possibility of implanting radiation sources inside objects, c) the Positron Emission Tomography (PET), which uses the possibility of implanting positron sources in objects and exploits the properties of electron-positron pair annihilation. Milestones of their sensational evolution throughout decades are vividly recalled in recent reviews, see e.g. [35, 47].

It is observed that all three cited imaging methods deal exclusively with primary (or non-deviated) radiation. The physical quantity, which is non-uniformly distributed over an object and responsible for the imaging process, is respectively: the linear attenuation coefficient for CT, the $\vartheta$-ray radio-activity density inside the object for SPECT, and the ${ϐ}^{+}$-ray (or positron) radio-activity density for PET.

However right at the start, scientists were also attracted by the idea of using scattered radiation by Compton effect (the scattering of X- or gamma photons by an electron) to image the inner parts of an object by reconstructing its electron density map. While this last imaging technology has not yet reached the same level of maturity as the three quoted above, it has gone through revolutionary conceptual steps which have led to what is nowadays known as Compton scatter tomography (CST). Yet it has stirred continuous interest in numerous applications, see e.g. [2], [4], [11], [18], [5], [22], [29], [1], [23]. The aim of this chapter is to recount the past episodes of research, to give a comprehensive account of what has been accomplished in CST and to describe some new ideas which have arisen recently, see [42, 54]. The emphasis of the discussion shall be placed at the theoretical level. Negative effects on imaging such as beam attenuation and multiple scattering, which complicate enormously the analytic treatment, shall be dealt with conventional retrieval or compensation methods. The crux of the matter is to see how fertile ideas evolve in time and generate new fruitful concepts.

2. Compton scattering

To appreciate the role of scattered radiation for imaging purposes, it is useful to recall some of the key points on Compton scattering. When a thin pencil of X- or gamma radiation shines through a medium, its intensity weakens as its traverses matter. One of the cause for this attenuation process is Compton scattering: deflected rays by scattering will not reach a detector placed along the incident direction.

From the geometry and kinematics of this scattering, see e.g. [6], one may infer that:

• the scattered photon flux density in a given spatial direction is given by the differential Compton scattering cross-section $d{ϡ}_C/d\mathrm{\chi }=Ϟr_e^{2}P\left(ϧ\right)$, where $r_e=2.82\times 10^{-15}$ m is the classical electron radius and $P\left(ϧ\right)$, the Klein-Nishina scattering probability under a scattering angle $ϧ$.

• the scattered photon energy $E$ is directly connected to the scattering angle $ϧ$ by the so-called Compton relation

$$E=E\left(ϧ\right)=E_0\frac{1}{1-ϓ\cos ϧ},$$E1

where $E_0$ is the energy of incident photons and $ϓ$ the ratio of $E_0$ to the electron rest energy.

Then the number of particles $d^2{N}_{\mathit{sc}}$ scattered in a solid angle $d{\mathrm{\chi }}_{\mathit{sc}}$ along a direction making an angle $ϧ$ with the incident direction follows from the definition of the differential scattering cross section $d{ϡ}_C/d\mathrm{\chi }$, if the following quantities are given: a) ${ϳ}_{\mathit{in}}$, the incident photon flux density, b) $n\left(\mathbf{M}\right)$, the electron density at the scattering site $\mathbf{M}$, surrounded by volume element $d\mathbf{M}$. Thus, for a given incident energy, the angular distribution of scattered photons around the scattering site $\mathbf{M}$ is no longer isotropic. The final form of this number of scattered photons in the direction given by the angle $ϧ$ is

$$d^2{N}_{\mathit{sc}}={ϳ}_{\mathit{in}}n\left(\mathbf{M}\right)d\mathbf{M}Ϟr_e^{2}P\left(ϧ\right)\mathrm{}d{\mathrm{\chi }}_{\mathit{sc}}\mathrm{.}$$E2

Figure 1.

Compton scattering

But this is still not realistic. For radiation emitted from a point source and incident on site $\mathbf{M}$, see Fig. 1, attenuation effects in matter before and after scattering are to be taken into account, as well as beam spreading due to straight line propagation, in the evaluation of the detected photon flux density. Equation 2 reads now for a point source emitting isotropically $I_0$ photons per second and per steradian

$$d^2{N}_{\mathit{sc}}=\frac{I_0}{4Ϟ}{m}_{ϡ}\left(\mathit{SM}\right)\mathrm{}{\mathcal{A}}_{\mathit{in}}\left(\mathit{SM}\right)\mathrm{}n\left(\mathbf{M}\right)d\mathbf{M}Ϟr_e^{2}P\left(ϧ\right)\mathrm{}{\mathcal{A}}_{\mathit{out}}\left(\mathit{MD}\right)\mathrm{}{m}_{ϡ\mathrm{\text{'}}}\left(\mathit{MD}\right)\mathrm{}d{\mathrm{\chi }}_{\mathit{sc}},$$E3

where the attenuation factors on the traveled distances $\mathit{SM}$ and $\mathit{MD}$ are given by

$${\mathcal{A}}_{\mathit{in}}\left(\mathit{SM}\right)=\exp \left(-{\cup }_0^{\mathit{SM}}\mathrm{‎}\mathit{ds}Ϛ(\mathbf{S}+s\mathbf{M})\right),\mathrm{}\mathit{and}{\mathcal{A}}_{\mathit{out}}\left(\mathit{MD}\right)=\exp \left(-{\cup }_0^{\mathit{MD}}\mathrm{‎}\mathit{ds}Ϛ(\mathbf{M}+s\mathbf{D})\right)\mathrm{.}$$E4

In equation 4, $Ϛ\left(\mathbf{M}\right)$ is the matter linear attenuation coefficient at site $\mathbf{M}$ and the respective beam spreading factor $m_{ϡ}\left(r\right)$ is of the form

$$m_{ϡ}\left(r\right)={\left(\frac{1}{\mathit{Ϟϡ}}{\tan }^{-1}\frac{ϡ}{2r}\right)}^2,$$E5

with $ϡ$ the linear size of scattering volume and $r$ the traversed distance. Note that for $ϡ\to 0$, $m_{ϡ}\left(r\right)\mathrm{~}1/r^2$, a well known photometric factor, which shall be used later. Equation 3 is fundamental to the image formation by scattered radiation.

When $E_0$ is chosen in such a way that competing events such as photoelectric absorption and pair creation are virtually absent, then Compton scattering becomes the main phenomena to be considered for imaging, and the relevant physical quantity is $n\left(\mathbf{M}\right)$, the electron density in matter at site $\mathbf{M}$, from which other quantities such as chemical composition may be deduced. As $n\left(\mathbf{M}\right)$ appears also in the attenuation process, its determination becomes quite involved since it must be retrieved from a non-linear expression in $n\left(\mathbf{M}\right)$. It is clear that, in CST imaging, $n\left(\mathbf{M}\right)$ plays the role of the attenuation map in CT, the gamma-ray activity density in SPECT and the beta-plus-ray activity density in PET.

3. The earlier CST modalities

Equation 3 has inspired many of the earlier CST modalities which shall be described below.

3.1. Point by point scanning CST

The simplest and earliest procedure is schematically described by Fig. 2. Both point source and detector are equipped with an axial collimator (lead cylindrical tube). Their axis lie in a plane and the detector is connected to a multichannel analyzer. When the axis of a collimated point-like detector is made to intersect the incident pencil beam, the intersection site is actually a scattering site of the Compton effect. So by measuring the scattered photon flux density at a given energy, by estimating the strength of the attenuation factors, and by evaluating the beam spreading factors one can obtain $n\left(\mathbf{M}\right)$. The process is then repeated for all sites in a trans-axial slice. It is interesting to note that equation 3 has been recast in the framework of an inverse problem by E. M. A. Hussein et al.[33], who have set up a discretization scheme to solve it.

Figure 2.

CST point by point scanning

At the turn of the twentieth century, the idea of using scattered gamma rays for investigating hidden structures in tissues seems to have been proposed for the first time by F. W. Spiers [52]. In 1959, P. G. Lale [40], realizing that radiographic images of organs fail to reveal their inner structure, has suggested a technique using Compton scattered radiation from a thin pencil of X-rays. Using equation 1, it is possible to move a collimated detector in space to measure the scattered radiation flux density and deduce the electron density at a precise point on the incoming pencil of X-rays. When this measurement can be performed in a planar slice of an object, it said that an image of the slice is obtained by Compton scatter tomography (CST). Of course this is clearly a point by point procedure which is quite time consuming.

An analogous problem arises in nuclear industry. There heat transfer research in circulating two-phase fluids requires knowledge and measurement of their densities without disturbing or arresting their flow. It was realized that penetrating gamma-radiation would be the most suitable for this purpose. As matter density is responsible for traversing radiation attenuation (since most of the attenuation occurs as a result of Compton scattering) measurements of radiation attenuation along linear paths crossing the section of a pipe can be easily done. D. Kershaw [36] has shown how the matter density distribution can be reconstructed from these measurements. In fact this problem is practically identical to the medical computed tomography which was developed simultaneously at that time and Kershaw has actually performed the mathematical inversion of the classical Radon transform, but most probably without being aware of the seminal works of J. Radon [49] and A. M. Cormack [15]. But it was N. N. Kondic [Kondic Hahn, 1970], who has first suggested the use of scattered radiation for obtaining the two-phase fluid density in a pipe.

In 1973, R. L. Clarke et al. [13] used both transmitted and scattered gamma-rays to measure bone mineral content. In [14], they have introduced the term gamma-tomography for their apparatus which is made of a fixed collimated pencil source and four focussing collimated detectors positioned so as the scattering angle is about $45^0$. The investigated object is placed on a moving platform which allows a raster scanning motion. The detected number of photons at each point (in fact a small "sensitive volume") is directly mapped onto the screen of an oscilloscope which yields a density image in a slice.

Later on, S. R. Gautam et al. [21] presented an imaging system called Compton Interaction Tomography, which consists of a pencil incident beam falling on a large object, the backscattered pencil is registered by a collimated detector situated on a line parallel to the large object. This is also a point by point determination of the electron density, but with the idea that the object is on one side of the line source-detector. Numerous works have been done by the NDE community in the 80’s, see e.g. [10], [27], [30], including some reviews such as [31],[32],[8], and the use of dual sources by [28].

3.2. Line by line CST

In 1971, F. T. Farmer & M. P. Collins [19] advocated the use of a wide angle collimated detector limited by two plates parallel to the incident radiation pencil. Fig. 3 shows how schematically how this procedure works. The detector is coupled to a multichannel analyzer. This design allows to obtain rapidly the electron density profile along the probing incident radiation pencil. This is the reason why it is called line scanning technique. A global image of an object slice is realized when parallel lines are put together. This technique has been refined in [20].

Figure 3.

CST with line by line scanning

Variants of this design have been developed by the NDE community to combat the low efficiency of the single voxel technique. In [25], G. Harding described a Compton scatter imaging system for NDE with horizontal scanning and incident pencil beam perpendicular to the scanning direction. Scatter detectors are disposed on both sides of the incident beam. Fig. 4 shows a sketch of the COMSCAN (Compton Scatter Scanner), developed by the Philips Research Laboratories in Hamburg, Germany which was quite successful (see [26]).

Figure 4.

Philips Research Laboratory Compton Scatter Scanner (COMSCAN)

3.3. Plane by plane CST

Then emerged several planar scatter imaging systems as extensions of linear scatter imaging systems in NDE [30], [53] as well as in medicine [24], which work at fixed scattering angle with relative success, since most of the reconstruction methods are still mathematically ill-defined and the technological problems not yet under control. Fig. 5 shows a plane by plane scanning at fixed scattering angle $ϧ$ using an movable collimated gamma camera. However this technique seemed to suffer from lack of sensitivity.

Figure 5.

CST with plane by plane scanning

By now, Compton scatter tomography as set up by these workers seemed to be well established with spatial resolutions of better than $1$cm and tissue density resolutions of better than $5/100$ for a radiation dose of less than $1$ rad. Yet two problems arose: reduction in image contrast and impaired tissue density resolution due to photon multiple scattering and attenuation artifacts. J. J. Battista et al. [7] has offered a way to cure for an improved Clarke scanner.

In the meantime the inversion problem has been formulated as a matrix inverse problem in [3] and [33], which was then solved numerically. It was realized that CST has considerable advantages over transmission CT in the low energy and low Z (atomic number) region, i.e. in medical range, and that CST has wide ranging NDE applications: viewing the interior of munitions, or structural material such as concrete, however with poor resolution because scanning was performed by hand moving the source-detector unit from pixel to pixel across the object. At this point, a remarkable leap forward, which provides a sound theoretical framework for CST, came into play.

4. Recent CST modalities

As early as 1978, a very original idea came from N. N. Kondic [38], who advocated the use of wide-angle collimators for both source and detector, instead of ray like collimators. Up to now, most proposed devices work with thin pencil-like beam source and with sharply (or widely) collimated detector coupled to a multichannel-analyzer. Kondic described a new way of irradiating a slice of an object by using a source’s collimator with a wide planar angle as well as recording scattered radiation with a detector’s collimator with wide angle. Thus the source takes the form of a fan beam and the wide angle collimator of the detector would delimitate an area in which an object is to be placed. In Fig. 6, we show a sketch of this proposal. Kondic made the crucial observation that when such a detector is connected to a multichannel analyzer, each measured energy channel, because of the Compton relation, is due to the sum of all scattering sites located on a circular arc starting from the source and ending at the detection site. Kondic called such an arc an isogonic line, wrote down the expression of the photon flux density registered at given scattering energy (including attenuation factors) and for the first time stated the electron density reconstruction problem in terms of integral data. Unfortunately a theoretical solution was not yet at hand. Expansion of this method and its variants were made later in [39]. Kondic cited also the advantages of using the isogonic lines. As collimation of the initial and scattered radiation in earlier modalities allows the detection of relatively few scattered photons, the counting statistics are very poor. The use of wide angle collimators would allow to work with higher counting rates which improves the data statistics and reduces the corresponding error, see [22]. Moreover Kondic did stress the need of varying the source - detector relative positions and (see item 6 on page 1146 of [39]), and proposed a rotation of the tomographic assembly in order to scan the whole field. Thus the photon flux density collected in an energy channel of a multichannel analyzer is due to the contribution not from a single scattering site but to the whole set of scattering sites corresponding to the same value of the scattering energy (or equivalently the same scattering angle $ϧ$). From elementary geometry, it can be seen that these sites lies on a circular arc passing through the source and the detection sites and subtending an inscribed angle of $(Ϟ-ϧ)$. The measured quantity is essentially an integral of the electron density over a curve in a first approximation whereby attenuation and other effects are set aside.

In 1985, elaborating on Kondic’s ideas, D. E. Bodette & A. M. Jacobs [9] were the first to realize that the scattered radiation data at constant energy may be expressed as an integral transform of the electron density on the isogonic line of Kondic. They pointed out the analogy with CT for which a measurement is also proportional to the integral of the linear attenuation coefficient along a straight line. They also conceded that the analytic derivation of a solution would be out of hand and offered a numerical treatment to extract a few first moments of the electron density from which several collective properties can be determined.

In 1993, starting from these ideas, T. H. Prettyman et al. [48] conceived a device which he called an energy dispersive projective scatterometer. Fig. 7 represents a sketch of this 1993 scatterometer device with two heavy wide angle collimators positioned at the gamma-ray source and at the detector. He then produced a numerical algorithm to deduce the electron density from the obtained measurements. He also pointed out that more complete data can be obtained by rotating the object, but one could have as well rotate the scatterometer. The coupling of the detector to a multichannel analyzer allows to collect the object projection data at constant scattering energy $E\left(ϧ\right)$ as an integral on an isogonic arc of circle. More importantly, he observed that "different projections can be obtained by rotating the sample.", hinting that, as in the case of CT, a large number of projections, generated by a spatial rotation, would be necessary for image reconstruction. However the concept of an integral transform is still missing in the formulation of this forward problem.

In a sectional slice, the electron density $n\left(\mathbf{M}\right)$ is a function of two coordinates. The necessary scatterometer data to recover $n\left(\mathbf{M}\right)$ should depend also on two variables. One of them obviously would be $ϧ$ the scattering angle (or alternatively the scattering energy $E\left(ϧ\right)$, which can be read off on a multichannel analyzer). But there are many ways to choose the second variable. We shall discuss three of them. The situation is very reminiscent of Synthetic Aperture Radar imaging which admits two approaches, one of them makes use of integrals of the ground reflectivity function over circles on a planar ground. Hence this approach turns the reconstruction problem into a problem of integral geometry in the sense of I. M. Gelfand, which consists of reconstructing a function on a two-dimensional space by its circular arc integrals.

In the sequel, we shall discuss three modalities based on three classes of circular arcs in the plane. They are illuminating examples of mathematical contributions to imaging science. Moreover as pointed out in S. J. Norton [44], the existence of an analytic inverse formula for an integral transform could provide the basis for a computationally efficient and robust imaging algorithm, which should be also resistant to noise.

Figure 6.

Sketch of Kondic’s 1978 proposal

Figure 7.

Sketch of Prettyman’s 1993 Scatterometer

The original idea of Kondic has led Bodette and Jacobs to conclude that a true CST should to be founded on a generalized Radon transform in the same way CT, SPECT and PET imaging have as mathematical foundation the classical Radon transform. In reality, SPECT and PET imaging need a more complicated transform called the attenuated Radon transform which includes an attenuation factor in the form of the exponential of a line integral to account for a realistic functioning. But the presence of such a factor has made the inversion problem extremely complicated even if the attenuation map is known a priori. After decades of effort, the inversion of this transform was finally and successfully performed by R. G. Novikov [46], thanks to an ingenious tour de force in complex analysis. Here we are facing the formidable inversion problem of an attenuated Radon transforms on circular arcs in the plane.

If a straight line in a plane is defined by two parameters, a circle requires three parameters: the two coordinates of its center and its radius. Thus to reconstruct a function of two variables such as $n(\mathbf{M}$), a condition should be imposed in order to bring down the circle parameters to two. But there exists an infinite number of ways to define a two-parameter circle in the plane. In the CST context, the scattering angle $ϧ$, or any related function of it, appears to be a necessary parameter. The second parameter may be taken from the rotational symmetry of the problem and is simply a polar angle $ϳ$ specifying the geometric position of the circle with respect to a reference direction. The integral transform may be then written down by summing equation 3 over the scattering sites $\mathbf{M}$. The result $\stackrel{\text{̠}}{n}(ϧ,ϳ)$ is called the attenuated Radon transform of $n\left(\mathbf{M}\right)$ on a circular arc. Then, for fixed $\mathbf{S}$ and $\mathbf{D}$, we have

$$\stackrel{\text{̠}}{n}(ϧ,ϳ)={\cup }_{{\mathbb{R}}^2}\mathrm{‎}d\mathbf{M}\frac{I_0}{4Ϟ}{m}_{ϡ}\left(\mathit{SM}\right)\mathrm{}{\mathcal{A}}_{\mathit{in}}\left(\mathit{SM}\right)\mathrm{}n\left(\mathbf{M}\right)\mathrm{}Ϟr_e^{2}P\left(ϧ\right)\mathrm{}{\mathcal{A}}_{\mathit{out}}\left(\mathit{MD}\right)\mathrm{}{m}_{ϡ\mathrm{\text{'}}}\left(\mathit{MD}\right)\mathrm{}\Upsilon \left(\mathrm{}\mathit{Circ}\mathrm{.}\mathit{Arc}\mathrm{}\right),$$E6

where $\Upsilon \left(\mathrm{}\mathit{Circ}\mathrm{.}\mathit{Arc}\mathrm{}\right)$ symbolizes the Dirac distribution concentrated on the chosen circular arc.

In the rest of the chapter, we shall discuss CST modalities, which are based on integral measurements, i.e. for which data are integrals of the electron density on arcs of circle linking the source point to the detection point in the plane. At present there are three of such CST modalities, corresponding to three choices of circle families in the plane. To remain on the track of the basic ideas, attenuation and beam spreading factors shall not be taken into account at first since solving such a general problem would be out of reach. Consequently we shall first concentrate on the problem of integral transform inversion. Attenuation and propagation spreading would be dealt with later. We shall see that in some cases the beam spreading effect can be included in the exact solution but attenuation must be either corrected or compensated for.

4.1. First CST modality (Norton 1994)

The first CST scanner which integrates the notion of Radon transform on Kondic isogonic arcs was proposed by S. J. Norton [44]. At that time it was known that A. M. Cormack [16] had established that the Radon transform on circles intersecting a fixed point in the plane is invertible. Norton proposed to use this result in the conception of his CST scanner, which was patented in 1995, see [45]. The Norton scanner has a gamma-ray source at fixed point $\mathbf{S}$ and a detector $\mathbf{D}$ movable on a line intersecting the source site, as sketched in Fig. 8.

Figure 8.

Fisrt CST modality

The point source $\mathbf{S}$ emits primary radiation towards an object, of which $\mathbf{M}$ is a running point. A point detector $\mathbf{D}$ moves along an $\mathit{Ox}$-axis and collects, at given energy $E\left(ϧ\right)$, scattered radiation from the object. The physics of Compton scattering demands that the registered radiation flux energy $\stackrel{\text{̠}}{n}\left(\mathbf{D}\right)$ at site $\mathbf{D}$ is due to the contribution of all scattering sites $\mathbf{M}$ lying on an arc of circle from $\mathbf{S}$ to $\mathbf{D}$ subtending an angle $(Ϟ-ϧ)$, where $ϧ$ is the scattering angle corresponding to the outgoing energy $E\left(ϧ\right)$, as given by equation 1.

Mathematically, $\stackrel{\text{̠}}{n}\left(\mathbf{D}\right)$ is essentially the integral of the object electron density $n\left(\mathbf{M}\right)$ on such arc of circles, when radiation attenuation and beam spreading effects due to radiation propagation are neglected. If polar coordinates with $\mathbf{S}$ as origin are used, then a running point $\mathbf{M}$ on a circle of diameter $p$, with center at the point of polar coordinates $(p/2,ϳ)$, is given by $(r,\varpi )$ with $r=p\cos (\varpi -ϳ)$. Thus, calling $\vartheta =(\varpi -ϳ)$ and recalling that the circle arc element is $\mathit{ds}=\mathit{pd}\vartheta$, we have $\stackrel{\text{̠}}{n}\left(\mathbf{D}\right)=\stackrel{\text{̠}}{n}(p,ϳ)$ with

$$\stackrel{\text{̠}}{n}(p,ϳ)={\cup }_{\mathit{arcSD}}\mathrm{‎}\mathit{dsn}(r,\vartheta +ϳ)=p{\cup }_{-}^{ϳ}\mathrm{‎}d\vartheta n(p\cos \vartheta ,\vartheta +ϳ),$$E7

where, for ease of notations, we have absorbed in the definition of $\stackrel{\text{̠}}{n}\left(\mathbf{D}\right)$ the Compton differential cross-section and the emitted flux density of the source and assumed no attenuation and beam spreading factors on the paths $\mathit{SM}$ and $\mathit{MD}$. Equation 7 is not really the circular Radon transform of $n\left(\mathbf{M}\right)$ in A. M. Cormack [16], since the integral goes over only the upper circular arc, which physically corresponds to the scattering angle $ϧ$ and the detected energy $E\left(ϧ\right)$. Note that the lower arc $\mathit{SD}$, is related to the scattering angle $(Ϟ-ϧ)$ and the scattering energy $E(Ϟ-ϧ)$, which is according to equation 1 not equal to $E\left(ϧ\right)$. Yet Norton argued that for an object situated above the line $\mathit{SD}$, which means that the support of $n(r,\varpi )$ is situated in the upper half-plane, the inversion procedure of A. M. Cormack should work. A closed form inverse formula, which mathematically well defined, has been given in [17] as

$$n(r,\varpi )=\frac{1}{2{Ϟ}^{2}r}{\cup }_0^2\mathrm{‎}\mathit{dϳ}{\cup }_0^{\mathrm{\propto }}\mathrm{‎}\mathit{dp}\frac{\complement \stackrel{\text{̠}}{n}(p,ϳ)}{\complement p}\frac{1}{r/p-\cos (\varpi -ϳ)}\mathrm{.}$$E8

In his work of 1994 [44], Norton has also derived an alternative inversion formula via the radial Fourier transform of $n(r,\varpi )$.

We now give some details on how equation 8 is derived. Assuming that both $n(r,\varpi )$ and $\stackrel{\text{̠}}{n}(p,ϳ)$ can be represented by their angular Fourier series

$$n(r,\varpi )=\coprod _l\mathrm{‎}{n}_l\left(r\right)\mathrm{}{e}^{\mathit{il}\varpi },\mathrm{}\mathit{with}\stackrel{\text{̠}}{n}(p,ϳ)=\coprod _l\mathrm{‎}{\stackrel{\text{̠}}{n}}_l\left(p\right)\mathrm{}{e}^{\mathit{ilϳ}},$$E9

then equation 7 takes the form of a Chebyshev

This name comes from the fact that the Chebyshev polynomial of the first kind is

𝑇𝑙(cos𝑥)=cos𝑙𝑥. transform, i.e.

$${\stackrel{\text{̠}}{n}}_l\left(p\right)=2{\cup }_0^p\mathrm{‎}\mathit{dr}\frac{\cos l\left({\cos }^{-1}(r/p)\right)}{\sqrt{1-(r/p{)}^2}}{n}_l\left(r\right),$$E10

which can be inverted using the following identity, (see [16),

$${\cup }_s^t\mathrm{‎}\frac{\mathit{dx}}{x}\frac{\cosh \left(l{\cosh }^{-1}\left(s/x\right)\right)}{\sqrt{{\left(s/x\right)}^2-1}}\frac{\cos l\left({\cos }^{-1}\right(t/x\left)\right)}{\sqrt{1-(t/x{)}^2}}=\frac{Ϟ}{2}\mathrm{.}$$E11

The inverse formula for $n_l\left(r\right)$ is given in [17] as

$$n_l\left(r\right)=\frac{1}{Ϟr}\left({\cup }_0^r\mathrm{‎}\mathit{dp}\frac{d{\stackrel{\text{̠}}{n}}_l\left(p\right)}{\mathit{dp}}\frac{{\left((r/p)-\sqrt{(r/p{)}^2-1}\right)}^l}{\sqrt{(r/p{)}^2-1}}-{\cup }_r^{\mathrm{\propto }}\mathrm{‎}\mathit{dp}\frac{d{\stackrel{\text{̠}}{n}}_l\left(p\right)}{\mathit{dp}}{U}_{l-1}(r/p)\right),$$E12

where $U_{l-1}\left(\cos x\right)=\sin \mathit{lx}/\sin x$. Then $n(r,\varpi )$ in equation 8 is obtained by resumming the series in equation 9.

4.2. Second CST modality (Nguyen-Truong 2010)

There are two scanning modes in this modality. If the object is "small", or can be put inside a circle of radius $p$ (a mechanical adjustable parameter), one can make an "internal" scanning inside this circle. If the object is "large" or situated far way from the observer, one can use the "external" scanning mode. The two modes shall be discussed separately.

4.2.1. Internal scanning

The second modality has originated from an proposal in [41]. At the time only numerical simulation and reconstruction using the Singular Value Decomposition method were performed on a turbine blade with encouraging results. One may conceive this modality as an evolved Prettyman’s scatterometer, which can rotate around an axis perpendicular to its plane, realizing what Prettyman had foreseen long ago to generate more useful data for his numerical reconstruction method. In fact, as will be shown the rotational motion is necessary for the reconstruction the electron density in a trans-axial slice.

Figure 9.

Second CST modality: internal scanning

The apparatus is is sketched in Fig. 9. An emitting radiation point source $\mathbf{S}$ is placed at a distance $2p$ from a point detector $\mathbf{D}$. The segment $\mathit{SD}$ joining them rotates around its middle point $\mathbf{O}$. At site $\mathbf{D}$ is collected the single-scattered radiation flux density from the scanned object for a given angular position of the line $\mathit{SD}$ and at a given scattering energy $E\left(ϧ\right)$, (equivalently at scattering angle $ϧ$). Then at fixed $ϳ$, a multichannel analyzer records the photon counts in each energy channel. Data acquisition is performed for every angular position $ϳ$ of $\mathit{SD}$.

Thus, thanks to the physics of the Compton effect, the detected radiation flux density is proportional to the integral of the electron density $n\left(\mathbf{M}\right)$ on a class of circular arcs sharing a chord of fixed length $2p$, which rotates about its fixed middle point $\mathbf{O}$.

With a polar coordinate system centered at $\mathbf{O}$, the equation of a circular arc lying inside a circle of center $O$ and radius $p$ reads (see [42])

$$r=r\left(\cos \right(\varpi -ϳ\left)\right)=p\left(\sqrt{1+{Ϣ}^2{\cos }^2(\varpi -ϳ)}-Ϣ\cos (\varpi -ϳ)\right),$$E13

which is defined by two parameters $Ϣ$ and $ϳ$: a) $ϳ$ is the angle made by its symmetry axis with the reference direction $\mathit{Ox}$ and b) $Ϣ$ is related to the scattering angle $ϧ$ by $Ϣ=\cot ϧ\mathrm{.}$ Note that $Ϣ$ is positive for $0<ϧ<Ϟ/2$ and the range of $\varpi$ is $(ϳ-Ϟ/2)<\varpi <(ϳ+Ϟ/2)$.

Again to simplify the notation by assuming that attenuation and beam spreading are neglected as well as by absorbing the Compton kinematic factor into one single function, we can say that the detected photon flux density $\stackrel{\text{̠}}{n}(Ϣ,ϳ)$ is the Radon transform of the electron density $n(r,\varpi )$ along arcs of circle given by equation 13. Thus using the auxiliary angle $\vartheta =(\varpi -ϳ)$ as in the previous subsection, we may express $\stackrel{\text{̠}}{n}(Ϣ,ϳ)$ with $\mathit{ds}$, the integration arc element given by

$$\mathit{ds}=\mathit{rd}\vartheta \sqrt{\frac{1+{Ϣ}^2}{1+{Ϣ}^2{\cos }^2\vartheta }}=\mathit{dr}\frac{\sqrt{1+{Ϣ}^2}}{Ϣ\sin \vartheta }\mathrm{.}$$E14

Now introducing the angular Fourier components of $\stackrel{\text{̠}}{n}(Ϣ,ϳ)$ and $n(r,\varpi )$ as given by equation 9, we see that they are related by

$$\frac{Ϣ{\stackrel{\text{̠}}{n}}_l\left(Ϣ\right)}{\sqrt{1+{Ϣ}^2}}=2{\cup }_p^{(}\mathrm{‎}\frac{\mathit{dr}}{\sqrt{1-\frac{1}{4{Ϣ}^2}{\left(\frac{p}{r}-\frac{r}{p}\right)}^2}}\cos \left[l{\cos }^{-1}\left(\frac{1}{2Ϣ}\left(\frac{p}{r}-\frac{r}{p}\right)\right)\right]n_l\left(r\right)\mathrm{.}$$E15

At first this Chebyshev transform looks a bit hopeless. However if one introduces a new variable $g$ defined by

$$g=\frac{1}{2}\left(\frac{p}{r}-\frac{r}{p}\right),$$E16

one can put it under the form

$$\frac{Ϣ{\stackrel{\text{̠}}{n}}_l\left(Ϣ\right)}{\sqrt{1+{Ϣ}^2}}=2{\cup }_0^{Ϣ}\mathrm{‎}\mathit{dg}\frac{\cos l{\cos }^{-1}\left(\frac{g}{Ϣ}\right)}{\sqrt{1-\frac{g^2}{{Ϣ}^2}}}\frac{p(\sqrt{1+g^2}-g)}{\sqrt{1+g^2}}{n}_l\left(p\right(\sqrt{1+g^2}-g\left)\right)\mathrm{.}$$E17

Then defining new functions by

$${\stackrel{\text{̠}}{N}}_l\left(Ϣ\right)=\frac{Ϣ{\stackrel{\text{̠}}{n}}_l\left(Ϣ\right)}{\sqrt{1+{Ϣ}^2}}\mathit{and}{N}_l\left(g\right)=\frac{p(\sqrt{1+g^2}-g)}{\sqrt{1+g^2}}{n}_l\left(p\right(\sqrt{1+g^2}-g\left)\right),$$E18

we obtain

$${\stackrel{\text{̠}}{N}}_l\left(Ϣ\right)=2{\cup }_0^{Ϣ}\mathrm{‎}\mathit{dg}\frac{\cos l{\cos }^{-1}\left(\frac{g}{Ϣ}\right)}{\sqrt{1-\frac{g^2}{{Ϣ}^2}}}{N}_l\left(g\right),$$E19

which is precisely of the form of equation 10, obtained in the Radon problem on circles passing through a fixed point. Hence an inversion formula exists in the $(Ϣ,g)$ variables, i.e.

$$N_l\left(g\right)=\frac{1}{Ϟg}\left({\cup }_0^g\mathrm{‎}\mathit{dp}\frac{d{\stackrel{\text{̠}}{N}}_l\left(Ϣ\right)}{dϢ}\frac{{\left((g/Ϣ)-\sqrt{(g/Ϣ{)}^2-1}\right)}^l}{\sqrt{(g/Ϣ{)}^2-1}}-{\cup }_g^{\mathrm{\propto }}\mathrm{‎}dϢ\frac{d{\stackrel{\text{̠}}{N}}_l\left(Ϣ\right)}{dϢ}{U}_{l-1}(g/Ϣ)\right)\mathrm{.}$$E20

A closed form of the inversion formula can be deduced from equation 8

$$N(g,\varpi )=\frac{1}{2{Ϟ}^{2}g}{\cup }_0^2\mathrm{‎}\mathit{dϳ}{\cup }_0^{\mathrm{\propto }}\mathrm{‎}dϢ\frac{\complement \stackrel{\text{̠}}{N}(Ϣ,ϳ)}{\complement Ϣ}\frac{1}{g/Ϣ-\cos (\varpi -ϳ)}\mathrm{.}$$E21

Finally going back to the original functions $n(r,\varpi )$ and $\stackrel{\text{̠}}{n}(Ϣ,ϳ)$, via equations 18, the reconstructed electron density is

$$n(r,\varpi )=\frac{1}{2{Ϟ}^2}\frac{\sqrt{1+\frac{1}{4}{\left(\frac{p}{r}-\frac{r}{p}\right)}^2}}{\frac{r}{2}\left(\frac{p}{r}-\frac{r}{p}\right)}{\cup }_0^2\mathrm{‎}\mathit{dϳ}{\cup }_0^{\mathrm{\propto }}\mathrm{‎}dϢ\frac{1}{\frac{1}{2Ϣ}\left(\frac{p}{r}-\frac{r}{p}\right)-\cos (\varpi -ϳ)}\frac{\complement }{\complement Ϣ}\left(\frac{Ϣ\stackrel{\text{̠}}{n}(Ϣ,ϳ)}{\sqrt{1+{Ϣ}^2}}\right)\mathrm{.}$$E22

4.2.2. External scanning

This scanning mode is illustrated by Fig. 10. Data acquisition is the same as for internal scanning, except for an object of compact support, the rotational motion may be replaced by a back and forth radar type sweeping motion.

Figure 10.

Second CST modality: external scanning

The working is similar to the internal scanning mode, except that one uses the scattering angle range $Ϟ/2<ϧ<Ϟ$ and the "external" arc of circle given by the equation

$$r=r\left(\cos \right(\varpi -ϳ\left)\right)=p\left(\sqrt{1+{Ϣ}^2{\cos }^2(\varpi -ϳ)}+Ϣ\cos (\varpi -ϳ)\right),$$E23

where $Ϣ=-\cot ϧ>0$ and $(ϳ-Ϟ/2)<\varpi <(ϳ+Ϟ/2)$. This arc of circle Radon transform has not yet been considered in [42]. Similarly the Radon transform of the electron density on this arc of circle is now expressed by the following Chebyshev integral transform for its angular components

$$\frac{Ϣ{\stackrel{\text{̠}}{n}}_l\left(Ϣ\right)}{\sqrt{1+{Ϣ}^2}}=2{\cup }_p^p\mathrm{‎}\frac{\mathit{dr}}{\sqrt{1-\frac{1}{4{Ϣ}^2}{\left(\frac{r}{p}-\frac{p}{r}\right)}^2}}{n}_l\left(r\right)\cos \left[l{\cos }^{-1}\left(\frac{1}{2Ϣ}\left(\frac{r}{p}-\frac{p}{r}\right)\right)\right]\mathrm{.}$$E24

This equation will take the form of equation 10, as equation 19, when the intermediate variable

$$g\mathrm{\text{'}}=\frac{1}{2}\left(\frac{r}{p}-\frac{p}{r}\right),$$E25

and the functions

$${\stackrel{\text{̠}}{N}}_l\left(Ϣ\right)=\frac{Ϣ{\stackrel{\text{̠}}{n}}_l\left(Ϣ\right)}{\sqrt{1+{Ϣ}^2}}\mathit{and}{N}_l\left(g\mathrm{\text{'}}\right)=\frac{p(\sqrt{1+{g\mathrm{\text{'}}}^2}+g\mathrm{\text{'}})}{\sqrt{1+{g\mathrm{\text{'}}}^2}}{n}_l\left(p\right(\sqrt{1+{g\mathrm{\text{'}}}^2}+g\mathrm{\text{'}}\left)\right),$$E26

are used. Finally, following the same steps as in the previous scanning mode, the reconstruction formula for $n(r,\varpi )$ reads

$$n(r,\varpi )=\frac{1}{2{Ϟ}^2}\frac{\sqrt{1+\frac{1}{4}{\left(\frac{r}{p}-\frac{p}{r}\right)}^2}}{\frac{r}{2}\left(\frac{r}{p}-\frac{p}{r}\right)}{\cup }_0^2\mathrm{‎}\mathit{dϳ}{\cup }_0^{\mathrm{\propto }}\mathrm{‎}dϢ\frac{1}{\frac{1}{2Ϣ}\left(\frac{r}{p}-\frac{p}{r}\right)-\cos (\varpi -ϳ)}\frac{\complement }{\complement Ϣ}\left(\frac{Ϣ\stackrel{\text{̠}}{n}(Ϣ,ϳ)}{\sqrt{1+{Ϣ}^2}}\right)\mathrm{.}$$E27

4.3. Third CST modality (Truong-Nguyen 2011)

In this subsection we discuss two scanning modes of a third CST modality proposed in [54]. This modality has originated from a search for curves in the plane such that the Radon transform on these curves may be reduced to a Chebyshev integral transform of the form of equation 10, for circles passing through a fixed point. It turns out that Radon transforms on arcs of circles orthogonal to a fixed circle do have this property. The pair source-detector still moves on a circle of radius $p$ (still an adjustable parameter) and centered at $\mathbf{O}$, the origin of a polar coordinate system. But their separation distance is no longer constant as in the second CST modality: it depends on the scattering angle $ϧ$. The positions of $\mathbf{S}$ and $\mathbf{D}$ are given by an opening angle ${\vartheta }_0$, measured from the symmetry axis of the circular arc.

4.3.1. Internal scanning

In this case, the Radon transform is defined on the following circular arc of equation

$$r=p\left(Ϣ\cos \right(\varpi -ϳ)-\sqrt{{Ϣ}^2{\cos }^2(\varpi -ϳ)-1}),$$E28

where $Ϣ=1/\cos {\vartheta }_0$ and $-{\vartheta }_0<(\varpi -ϳ)<{\vartheta }_0$. Inspection of Fig. 11 shows that ${\vartheta }_0=(Ϟ/2-ϧ)$ for $0<ϧ<Ϟ/2$. Hence $Ϣ=1/\sin ϧ>1$. So data acquisition works as follows. For fixed $ϳ$, the pair source - detector must be simultaneously displaced on the circle of radius $p$ so that the opening angle $\stackrel{\text{̠}}{\mathit{SOD}}=2{\vartheta }_0=(Ϟ-2ϧ)$, before registering at $\mathbf{D}$ a photon count. The number of isogonic lines is thus equal to the number of positions of the pair $(\mathbf{S},\mathbf{D})$.

Figure 11.

Third CST modality: internal scanning

The integration arc element being now

$$\mathit{ds}=\mathit{rd}\vartheta \sqrt{\frac{{Ϣ}^2-1}{{Ϣ}^2{\cos }^2\vartheta }-1}=\mathit{dr}\frac{\sqrt{{Ϣ}^2-1}}{Ϣ\sin \vartheta },$$E29

the Radon transform of the electron density $\stackrel{\text{̠}}{n}(Ϣ,ϳ)$ becomes, in terms of the angular components $n_l\left(r\right)$ and ${\stackrel{\text{̠}}{n}}_l\left(Ϣ\right)$, a Chebyshev transform similar to equation 17, i.e.

$$\frac{Ϣ{\stackrel{\text{̠}}{n}}_l\left(Ϣ\right)}{\sqrt{{Ϣ}^2-1}}=2{\cup }_p^{(}\mathrm{‎}\frac{\mathit{dr}}{\sqrt{1-\frac{1}{4{Ϣ}^2}{\left(\frac{p}{r}+\frac{r}{p}\right)}^2}}{n}_l\left(r\right)\cos \left[l{\cos }^{-1}\left(\frac{1}{2Ϣ}\left(\frac{p}{r}+\frac{r}{p}\right)\right)\right]\mathrm{.}$$E30

The structural similarity with equations 15 and 24 suggests an intermediate variable $g"$ of the form

$$g"=\frac{1}{2}\left(\frac{r}{p}+\frac{p}{r}\right),$$E31

which yields the following Chebyshev transform

$$\frac{Ϣ{\stackrel{\text{̠}}{n}}_l\left(Ϣ\right)}{\sqrt{{Ϣ}^2-1}}=2{\cup }_1^{Ϣ}\mathrm{‎}\mathit{dg}"\mathrm{}\frac{\cos l{\cos }^{-1}\left(\frac{g"}{Ϣ}\right)}{\sqrt{1-\frac{g{"}^2}{{Ϣ}^2}}}\frac{p(g"-\sqrt{g{"}^2-1})}{\sqrt{g{"}^2-1}}{n}_l\left(p\right(g"-\sqrt{g{"}^2-1}\left)\right)\mathrm{.}$$E32

Now redefining the functions by

$${\stackrel{\text{̠}}{N}}_l\left(Ϣ\right)=\frac{Ϣ{\stackrel{\text{̠}}{n}}_l\left(Ϣ\right)}{\sqrt{{Ϣ}^2-1}}\mathit{and}{N}_l(g")=\frac{p\left(\sqrt{g"-g{"}^2-1}\right)}{\sqrt{g{"}^2-1}}{n}_l\left(p\right(g"-\sqrt{g{"}^2-1}\left)\right),$$E33

we obtain the form of equation 10, nevertheless with a lower integration bound equal to $Ϣ=1$ and not zero. This does not spoil the inversion procedure based on the identity 11, as shown in [54]. Consequently we recover the electron density in closed form as

$$n(r,\varpi )=\frac{1}{2{Ϟ}^2}\frac{\sqrt{\frac{1}{4}{\left(\frac{p}{r}+\frac{r}{p}\right)}^2-1}}{\frac{r}{2}\left(\frac{p}{r}+\frac{r}{p}\right)}{\cup }_0^2\mathrm{‎}\mathit{dϳ}{\cup }_0^{\mathrm{\propto }}\mathrm{‎}dϢ\frac{1}{\frac{1}{2Ϣ}\left(\frac{p}{r}+\frac{r}{p}\right)-\cos (\varpi -ϳ)}\frac{\complement }{\complement Ϣ}\left(\frac{Ϣ\stackrel{\text{̠}}{n}(Ϣ,ϳ)}{\sqrt{{Ϣ}^2-1}}\right)\mathrm{.}$$E34

4.3.2. External scanning

This scanning mode is shown in Fig. 12. It has the same data acquisition procedure as for internal scanning. It is appropriate for large objects or objects situated far away from the scanner, e.g. buried objects underground or undersea, with a back and forth radar type of scanning motion.

Figure 12.

Third CST modality: external scanning

The Radon transform is now defined on the external circular arc with respect to the reference circle of radius $p$ of equation

$$r=p\left(Ϣ\cos \right(\varpi -ϳ)+\sqrt{{Ϣ}^2{\cos }^2(\varpi -ϳ)-1})\mathrm{.}$$E35

Here $Ϟ/2<ϧ<Ϟ$and ${\vartheta }_0=(ϧ-Ϟ/2)$

Note that the two arcs of equations 28 and 35 belong to the same circle.

. As for the case of internal scanning, we have also $Ϣ=1/\cos {\vartheta }_0$ and $-{\vartheta }_0<(\varpi -ϳ)<{\vartheta }_0$ with the same data acquisition procedure.

The integration arc element is the same as in the previous scanning mode, see equation 29. Then the Radon transform of the electron density $\stackrel{\text{̠}}{n}(Ϣ,ϳ)$ becomes, in terms of function angular components, the following Chebyshev transform

$$\frac{Ϣ{\stackrel{\text{̠}}{n}}_l\left(Ϣ\right)}{\sqrt{{Ϣ}^2-1}}=2{\cup }_p^p\mathrm{‎}\frac{\mathit{dr}}{\sqrt{1-\frac{1}{4{Ϣ}^2}{\left(\frac{p}{r}+\frac{r}{p}\right)}^2}}{n}_l\left(r\right)\cos \left[l{\cos }^{-1}\left(\frac{1}{2Ϣ}\left(\frac{p}{r}+\frac{r}{p}\right)\right)\right]\mathrm{.}$$E36

We see that the intermediate variable $g"$ of equation (31) can be used again so that

$$\frac{Ϣ{\stackrel{\text{̠}}{n}}_l\left(Ϣ\right)}{\sqrt{{Ϣ}^2-1}}=2{\cup }_1^{Ϣ}\mathrm{‎}\mathit{dg}"\mathrm{}\frac{\cos l{\cos }^{-1}\left(\frac{g"}{Ϣ}\right)}{\sqrt{1-\frac{g{"}^2}{{Ϣ}^2}}}\frac{p(g"+\sqrt{g{"}^2-1})}{\sqrt{g{"}^2-1}}{n}_l\left(p\right(g"+\sqrt{g{"}^2-1}\left)\right)\mathrm{.}$$E37

Now redefining the functions by

$${\stackrel{\text{̠}}{N}}_l\left(Ϣ\right)=\frac{Ϣ{\stackrel{\text{̠}}{n}}_l\left(Ϣ\right)}{\sqrt{{Ϣ}^2-1}}\mathit{and}{N}_l(g")=\frac{p\left(\sqrt{g"+g{"}^2-1}\right)}{\sqrt{g{"}^2-1}}{n}_l\left(p\right(g"+\sqrt{g{"}^2-1}\left)\right),$$E38

we obtain the form of equation 10, and surprisingly the same equation as for internal scanning. Consequently one end up with the same reconstruction formula (34), which is a remarkable advantage for this modality. This is due to the fact that the two arcs are on the same circle.

5. Numerical simulations

This section is devoted to showing that the analytic formulas of inversion do lead to robust computation algorithms. As we have shown, it is sufficient develop a single reconstruction algorithm for the Radon transform on circles passing through a fixed point, since all other cases can be reduced to that one by a change of variable and a change of functions. Since we use the formalism of angular Fourier components, we shall adapt a procedure set up long ago by C. H. Chapman & P. W. Carey [12] for the classical Radon transform.

5.1. Algorithm

The idea is to start with the result of A. M. Cormack [15] (see his equation $18$b) for the reconstructed angular Fourier component of $n(r,\varpi )$ given by

$$n_l\left(r\right)=\frac{1}{Ϟr}\left({\cup }_0^r\mathrm{‎}\mathit{dp}\frac{d{\stackrel{\text{̠}}{n}}_l\left(p\right)}{\mathit{dp}}\frac{{\left((r/p)-\sqrt{(r/p{)}^2-1}\right)}^l}{\sqrt{(r/p{)}^2-1}}-{\cup }_r^{\mathrm{\propto }}\mathrm{‎}\mathit{dp}\frac{d{\stackrel{\text{̠}}{n}}_l\left(p\right)}{\mathit{dp}}{U}_{l-1}(r/p)\right)\mathrm{.}$$E39

We compute first $n_l\left(r\right)$ from a discretized version of equation 39, which shall be set up now, and get $n(r,\varpi )$ back by its angular Fourier series using equation 9.

To this end, we change the form of equation 39 by making the following change of variables

$$p=\frac{r}{\cosh x}\mathit{andp}=\frac{r}{\cos x},$$E40

respectively in the first and in the second integral of (39). Hence

$$n_l\left(r\right)=\frac{1}{Ϟ}\left({\cup }_1^{\mathrm{\propto }}\mathrm{‎}\mathit{dx}\frac{d{\stackrel{\text{̠}}{n}}_l\left(p\right)}{\mathit{dp}}\frac{e^{-\mathit{lx}}}{{\cosh }^{2}x}-{\cup }_0^{Ϟ}\mathrm{‎}\mathit{dx}\frac{d{\stackrel{\text{̠}}{n}}_l\left(p\right)}{\mathit{dp}}\frac{\sin \mathit{lx}}{{\cos }^{2}x}\right)\mathrm{.}$$E41

Then these integrals will be approximated by discrete sums in which we take for simplicity the same step $\Upsilon$ for the variable $r$ as well as for the variable $p$. Hence $n_l\left(r\right)\eqsim n_l\left(j\Upsilon \right),$ with

$$n_l\left(j\Upsilon \right)=\frac{1}{Ϟ}\left(\coprod _k^{=}\mathrm{‎}\frac{d{\stackrel{\text{̠}}{n}}_l\left(p\right)}{\mathit{dp}}{\cup }_{x_j^k}^{x_j^{k+1}}\mathrm{‎}\mathit{dx}\frac{e^{-\mathit{lx}}}{{\cosh }^{2}x}-\coprod _k^{=}\mathrm{‎}\frac{d{\stackrel{\text{̠}}{n}}_l\left(p\right)}{\mathit{dp}}{\cup }_{x_j^k}^{x_j^{k+1}}\mathrm{‎}\mathit{dx}\frac{\sin \mathit{lx}}{{\cos }^{2}x}\right),$$E42

where $K$ is the maximal value taken by $k$. Now introducing the indefinite integrals

$$I_l\left(x\right)=\cup \mathrm{‎}\mathit{dz}\frac{e^{-\mathit{lz}}}{{\cosh }^{2}z}\mathit{and}{J}_l\left(x\right)=\cup \mathrm{‎}\mathit{dz}\frac{\sin \mathit{lz}}{{\cos }^{2}z},$$E43

we see that they obey the following recursion relations, see [50]

$$(l-2)I_l\left(x\right)=2\mathrm{}\left(\tan x\sin (l-2)x-\cos (l-2)x\right)-l{I}_{l-2}\left(x\right),$$E44

with $I_0\left(x\right)=1/\cos x$ and $I_1\left(x\right)=-2\ln \left(\cos x\right)$, and

$$l{J}_{l+1}\left(x\right)=-\frac{1}{x^l(1+x^2}-(l+2)J_{l-1}\left(x\right),$$E45

with

$$J_0\left(x\right)=\frac{x}{2(1+x^2)}+\frac{1}{2}{\tan }^{-1}\mathit{xand}{J}_1\left(x\right)=\ln \left(\frac{x}{\sqrt{1+x^2}}\right)+\frac{1}{2(1+x^2)}\mathrm{.}$$E46

The derivatives of $d{\stackrel{\text{̠}}{n}}_l\left(p\right)/\mathit{dp}$ are replaced by a simple linear interpolation ${\stackrel{\mathrm{\text{'}}}{n}}_l\left(k\right)$, i.e.

$$\frac{d{\stackrel{\text{̠}}{n}}_l\left(p\right)}{\mathit{dp}}\eqsim \frac{n_l(k+1)-n_l\left(k\right)}{\Upsilon }={\stackrel{\mathrm{\text{'}}}{n}}_l\left(k\right)\mathrm{.}$$E47

Then the discretized value of the reconstructed angular component of $n_l\left(r\right)$, with $r\mathrm{~}j\Upsilon$, is

$$n_l\left(j\Upsilon \right)=\frac{1}{Ϟ}\left(\coprod _k^{=}\mathrm{‎}{\stackrel{\mathrm{\text{'}}}{n}}_l\left(k\right)\left(I_l\left(\right(k+1\left)\Upsilon \right)-I_l\left(k\Upsilon \right)\right)-\coprod _k^{=}\mathrm{‎}{\stackrel{\mathrm{\text{'}}}{n}}_l\left(k\right)\left(J_l\left(\right(k+1\left)\Upsilon \right)-J_l\left(k\Upsilon \right)\right)\right)\mathrm{.}$$E48

Equation 48 shall be used in numerical simulations

For difficulties on the handling of the inversion formula of S. J. Norton, see [55]

.

5.2. Simulation results

In this subsection, we present numerical simulations on the first CST modality applied to the Shepp-Logan medical phantom, see [50]. The source $\mathbf{S}$ is placed below on the left of the image and the detector moves along the line $\mathit{SD}$. The relevant space is represented by $256\times 256$ (length unit)${}^2$. Let $N_{ϳ}$ be the number of angular steps and let $N_p$ be the number of radii for circles going through the origin $O$. Then the following sampling steps ${\Upsilon }_{ϳ}=2Ϟ/N_{ϳ}\mathit{and}\Upsilon =4\mathrm{å}256/N_p$ are taken. Since there is no rotation around the object, the $(ϳ,p)$-space is very large in comparison to the studied image. To have a good representation of the object, we need a large maximum value of $p$, (e.g. four times the image size).

To estimate the reconstruction quality, we use the normalized mean square error ($\mathit{NMSE}$) and the normalized mean absolute error ($\mathit{NMAE}$) (expressed as a percentage), defined by

$\mathit{NMSE}=\frac{100}{N^2}\frac{\underset{(i,j)\nabla [1,N{]}^2}{\coprod \mathrm{‎}}{\left|{\mathcal{I}}_r(i,j)-{\mathcal{I}}_o(i,j)\right|}^2}{\underset{(i,j)\nabla [1,N{]}^2}{\max }\left\{{\mathcal{I}}_o\right(i,j){\}}^2}\mathit{andNMAE}=\frac{100}{N^2}\frac{\underset{(i,j)\nabla [1,N{]}^2}{\coprod \mathrm{‎}}\left|{\mathcal{I}}_r(i,j)-{\mathcal{I}}_o(i,j)\right|}{\underset{(i,j)\nabla [1,N{]}^2}{\max }\left\{{\mathcal{I}}_o\right(i,j\left)\right\}},$

where ${\mathcal{I}}_r$ is the reconstructed image and ${\mathcal{I}}_o$ is the original image.

In order to reduce the artifacts generated around the object by the circular harmonic decomposition approach, we can use information given by the contours of the projection data. If the projection data are vanishing, this means that the corresponding circular arc does not intersect the object. Therefore if the object of interest is bounded in space, a null set in the projection space corresponds to a null set in the original space. This approach gives a very interesting image quality in the reconstruction of the Shepp-Logan phantom. Contours and small structures are nicely recovered. A maximum value $p_{\mathit{max}}$ is to be fixed in numerical computation. But a sharp cut-off of $p$ and the ensuing loss of data generate significant artifacts. Moreover an increase of $p_{\mathit{max}}$ leads to an increase of the detector length in the first CST modality.

The original Shepp-Logan phantom is given in Fig. 13.

Figure 13.

Shepp-Logan original phantom, reprinted from [51]

Fig. 14 illustrates the first CST modality data, obtained by application of the Radon transform on circles passing through the coordinate origin. Finally the reconstructed image by the first CST modality is given in Fig. 15.

Figure 14.

Data of the first CST modality (Norton 1994), reprinted from [51]

Figure 15.

Shepp-Logan phantom reconstructed by the first CST modality (Norton 1994), reprinted from [51]

Figure 16.

Second CST modality simulation results illustration

Now using the same algorithm but in different function spaces, one can achieve the reconstruction with the two other CST modalities. With the second CST modality for small objects, collected data and image reconstruction are given in Fig. 16, see [51].

We see that generally a reasonable image quality in the reconstruction of the medical phantom is obtained. Contours and small structures are correctly recovered. The numerical error measurements are very close to those of the ordinary Radon transform and are even better for a medical phantom. Artifacts arise if the studied object occupies the whole medium since the boundary parts are not well scanned and less information are available from these parts. As the set up works with radial parameters, an object of circular symmetry fits better to this modality: this is the case of a medical phantom.

As for the CST third modality for small objects, the simulations results are displayed in Fig. 17, see [43].

Figure 17.

Third CST modality simulation results illustration

The modified Chapman & Carey approach used here gives in general a reasonable image quality in the reconstruction of the medical phantom (Fig. 17) with (Normalized Mean Absolute Error =2.2% (NMAE) and Normalized Mean Square Error = 0.038 (NMSE). Contours and small structures are well recovered. The numerical error measurements obtained using the algorithm described above are of the same order of magnitude as in existing scanning devices. Nevertheless we observe some artifacts. They can be reduced by a better acquisition of the data. Indeed the $(ϳ,Ϣ)$-space is very large but the $(ϳ,ϧ)$-space is bounded. Hence with a good resolution, the loss of data is less important.

5.3. Open problems

Of course what is presented here does not claim neither completeness nor perfection. There are still many important problems to tackle.

1. Attenuation. This question, as mentioned before, has no analytic answer, even in the near future. One way to solve practically the problem is to produce an approximate iterative corrective algorithm, see e.g. [51].

2. Beam spreading. In the limit of vanishingly small detecting pixel, the inclusion of beam spreading factors proportional to the inverse distances (or inverse square distances) $S{M}^{-1,-2}$ (or $M{D}^{-1,-2}$) does not spoiled the inversion procedure. One may redefine an "effective" electron density and an "effective" Radon transform. This is due to the fact that the product $(\mathit{SM}\times \mathit{MD}{)}^{-1,-2}$ is of the form of a product of a function of $Ϣ$ and a function of $r$, see [42, 54].

3. Multiple scattering. This degrading factor is due to the use of wide angle collimator and is difficult to be removed by theoretical means. As in other gamma-ray imaging systems, such as Single Photon Emission Computed Tomography (SPECT) and Positron Emission Tomography (PET), this problem can be only approximatively fixed by corrective measures adapted to each particular scanner.

4. Noise. Radiation processes carry fluctuations as noise which could drastically affect numerical simulations in image reconstruction. So algorithm robustness should be tested against noise systematically. There is no universal receipt and an adapted de-noising procedure is to be set up for each scanner type. An example of such de-noising procedure is given in [51].

6. Conclusion and perspectives

Compton scatter tomography (CST) is an outstanding example of the use of Compton scattered radiation for non-invasive imaging. CST provides high resolution and high sensitivity imaging with virtually any material, including lightweight structures and organic matter, which normally pose problems in conventional x-ray computed tomography because of low contrast. However it also suffers from the usual degrading factors of radiation imaging, such as multiple scattering, attenuation, noise and fluctuations, etc., which have to be treated by appropriate (existing or to be developed) methods. Yet it has undergone a very impressive theoretical evolution. In this chapter, we have shown how the expression of the differential cross-section of the Compton effect has led to the idea of measuring the value of matter electron density point by point. The ensuing progress brought to light the successive line by line and plane by plane scanning procedures, however with limited success. A quantum leap in redefining CST was made in the 80’s by N. N. Kondic, who advocated the use or wide angle collimators and pointed out the occurrence of integral measurements along circular arcs corresponding to a definite measured photon energy. Then a concept of integral transform of Radon type has emerged from this type of integral measurement in [9]. Later in the 90’s, elaborating on a proposal of N. N. Kondic, T. H. Prettyman has constructed a scanner on this principle and called an energy dispersive projective scatterometer. He also realized that a relative motion between this apparatus and the object is necessary to generate the data needed for image reconstruction. This is the starting point for recent CST modalities which are now based on Radon transform on Kondic isogonic circular arcs. We have presented and discussed three existing CST modalities based on three "circular-arc" Radon transforms. Two of them are of recent origin and have arisen as new elements of integral geometry. The remarkable aspect is that they share a common inversion method based on the inversion of the Radon transform on circles intersecting a fixed point. As numerical simulations show their feasibility and viability as imaging systems, CST appears as a promising area of research and development which could bring interesting advances in the field of non-invasive imaging for medicine and industry. Work towards more efficient ways for data processing, e.g. setting faster back-projection types of reconstruction methods, is underway in this context.