**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

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

Recent Developments on Compton Scatter Tomography: Theory and Numerical Simulations 107

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

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

framework for CST, came into play.

whereby attenuation and other effects are set aside.

density from which several collective properties can be determined.

**4. Recent CST modalities**

**Figure 4.** Philips Research Laboratory Compton Scatter Scanner (COMSCAN)

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 1cm 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**

6 Will-be-set-by-IN-TECH

**Figure 4.** Philips Research Laboratory Compton Scatter Scanner (COMSCAN)

technique seemed to suffer from lack of sensitivity.

**Figure 5.** CST with plane by plane scanning

scanner.

at fixed scattering angle *ω* using an movable collimated gamma camera. However this

By now, Compton scatter tomography as set up by these workers seemed to be well established with spatial resolutions of better than 1cm 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

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

#### 8 Will-be-set-by-IN-TECH 108 Numerical Simulation – From Theory to Industry Recent Developments on Compton Scatter Tomography: Theory and Numerical Simulations <sup>9</sup>

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*(*ω*) 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*(**M**) is a function of two coordinates. The necessary scatterometer data to recover *n*(**M**) should depend also on two variables. One of them obviously would be *ω* the scattering angle (or alternatively the scattering energy *E*(*ω*), 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.

**Figure 6.** Sketch of Kondic's 1978 proposal

**Figure 7.** Sketch of Prettyman's 1993 Scatterometer

*<sup>n</sup>*(*ω*, *<sup>φ</sup>*) =

**R**<sup>2</sup>

*<sup>d</sup>***<sup>M</sup>** *<sup>I</sup>*<sup>0</sup> 4*π*

summing equation 3 over the scattering sites **<sup>M</sup>**. The result *<sup>n</sup>*(*ω*, *<sup>φ</sup>*) is called the attenuated

where *δ*(Circ.Arc) 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

2

Recent Developments on Compton Scatter Tomography: Theory and Numerical Simulations 109

*<sup>e</sup> P*(*ω*) A*out*(*MD*) *mσ*�(*MD*) *δ*(Circ.Arc),

(6)

Radon transform of *n*(**M**) on a circular arc. Then, for fixed **S** and **D**, we have

*mσ*(*SM*) A*in*(*SM*) *n*(**M**) *π r*

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.

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*(**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

**Figure 6.** Sketch of Kondic's 1978 proposal

8 Will-be-set-by-IN-TECH

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*(*ω*) 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

In a sectional slice, the electron density *n*(**M**) is a function of two coordinates. The necessary scatterometer data to recover *n*(**M**) should depend also on two variables. One of them obviously would be *ω* the scattering angle (or alternatively the scattering energy *E*(*ω*), 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

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

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

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*(**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

inversion problem of an attenuated Radon transforms on circular arcs in the plane.

missing in the formulation of this forward problem.

imaging algorithm, which should be also resistant to noise.

integrals.

**Figure 7.** Sketch of Prettyman's 1993 Scatterometer

summing equation 3 over the scattering sites **<sup>M</sup>**. The result *<sup>n</sup>*(*ω*, *<sup>φ</sup>*) is called the attenuated Radon transform of *n*(**M**) on a circular arc. Then, for fixed **S** and **D**, we have

$$
\widehat{m}(\omega,\phi) = \int\_{\mathbb{R}^2} d\mathbf{M} \, \frac{I\_0}{4\pi} \, m\_\sigma(SM) \, \mathcal{A}\_{\text{int}}(SM) \, n(\mathbf{M}) \, \pi \, r\_\varepsilon^2 P(\omega) \, \mathcal{A}\_{\text{out}}(MD) \, m\_{\sigma'}(MD) \, \delta(\text{Circ.Arc}), \tag{6}
$$

where *δ*(Circ.Arc) 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

#### 10 Will-be-set-by-IN-TECH 110 Numerical Simulation – From Theory to Industry Recent Developments on Compton Scatter Tomography: Theory and Numerical Simulations <sup>11</sup>

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.

arc element is *ds* <sup>=</sup> *p dγ*, we have *<sup>n</sup>*�(**D**) = *<sup>n</sup>*�(*p*, *<sup>φ</sup>*) with

arc *SD*

*ds n*(*r*, *γ* + *φ*) = *p*

where, for ease of notations, we have absorbed in the definition of *<sup>n</sup>*�(**D**) the Compton differential cross-section and the emitted flux density of the source and assumed no attenuation and beam spreading factors on the paths *SM* and *MD*. Equation 7 is not really the circular Radon transform of *n*(**M**) 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*(*ω*). Note that the lower arc *SD*, is related to the scattering angle (*π* − *ω*) and the scattering energy *E*(*π* − *ω*), which is according to equation 1 not equal to *E*(*ω*). Yet Norton argued that for an object situated above the line *SD*, which means that the support of *n*(*r*, *θ*) 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

> *dp <sup>∂</sup>n*�(*p*, *<sup>φ</sup>*) *∂p*

*il<sup>θ</sup>* , with *<sup>n</sup>*�(*p*, *<sup>φ</sup>*) = ∑

cos−1(*r*/*p*)

�*l*

�<sup>1</sup> <sup>−</sup> (*r*/*p*)<sup>2</sup>

�

�(*r*/*p*)<sup>2</sup> <sup>−</sup> <sup>1</sup> <sup>−</sup>

where *Ul*−1(cos *<sup>x</sup>*) = sin *lx*/ sin *<sup>x</sup>*. Then *<sup>n</sup>*(*r*, *<sup>θ</sup>*) in equation 8 is obtained by resumming the

In his work of 1994 [44], Norton has also derived an alternative inversion formula *via* the

We now give some details on how equation 8 is derived. Assuming that both *n*(*r*, *θ*) and

1 *r*/*p* − cos(*θ* − *φ*)

*l*

�

cos *l*(cos−1(*t*/*x*)) �<sup>1</sup> <sup>−</sup> (*t*/*x*)<sup>2</sup> <sup>=</sup> *<sup>π</sup>*

> � ∞ *r*

*dp dn*�*l*(*p*)

*<sup>n</sup>*�*l*(*p*)*<sup>e</sup>*

� *π*/2 −*φ*

Recent Developments on Compton Scatter Tomography: Theory and Numerical Simulations 111

*dγ n*(*p* cos *γ*, *γ* + *φ*), (7)

. (8)

*ilφ*, (9)

<sup>2</sup> . (11)

⎞

⎟⎠ , (12)

*nl*(*r*), (10)

*dp Ul*−1(*r*/*p*)

*<sup>n</sup>*�(*p*, *<sup>φ</sup>*) = �

*<sup>n</sup>*(*r*, *<sup>θ</sup>*) = <sup>1</sup>

radial Fourier transform of *n*(*r*, *θ*).

2*π*2*r*

*<sup>n</sup>*�(*p*, *<sup>φ</sup>*) can be represented by their angular Fourier series

*l*

then equation 7 takes the form of a Chebyshev <sup>1</sup> transform, *i.e.*

*<sup>n</sup>*�*l*(*p*) = <sup>2</sup>

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

cosh �

� (*s*/*x*)

�

*n*(*r*, *θ*) = ∑

� *t s dx x*

*nl*(*r*) = <sup>1</sup>

series in equation 9.

*π r*

⎛

⎜⎝ � *r* 0

The inverse formula for *nl*(*r*) is given in [17] as

*dp dn*�*l*(*p*) *dp*

� 2*π* 0

*nl*(*r*)*e*

� *p* 0

*l* cosh−<sup>1</sup> (*s*/*x*)

<sup>2</sup> <sup>−</sup> <sup>1</sup>

(*r*/*p*) <sup>−</sup> �(*r*/*p*)<sup>2</sup> <sup>−</sup> <sup>1</sup>

<sup>1</sup> This name comes from the fact that the Chebyshev polynomial of the first kind is *Tl*(cos *x*) = cos *lx*.

*dr* cos *<sup>l</sup>* �

*dφ* � ∞ 0
