**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 18b) for the reconstructed angular Fourier component of *n*(*r*, *θ*) given by

$$m\_l(r) = \frac{1}{\pi r} \left( \int\_0^r dp \, \frac{d\widehat{n}\_l(p)}{dp} \, \frac{\left( (r/p) - \sqrt{(r/p)^2 - 1} \right)^l}{\sqrt{(r/p)^2 - 1}} - \int\_r^\infty dp \, \frac{d\widehat{n}\_l(p)}{dp} \, \mathcal{U}\_{l-1}(r/p) \right). \tag{39}$$

We compute first *nl*(*r*) from a discretized version of equation 39, which shall be set up now, and get *n*(*r*, *θ*) 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 \mathfrak{x}} \quad \text{and} \quad p = \frac{r}{\cos \mathfrak{x}'} \tag{40}$$

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

$$m\_l(r) = \frac{1}{\pi} \left( \int\_1^\infty dx \, \frac{d\hat{n}\_l(p)}{dp} \, \frac{e^{-lx}}{\cosh^2 x} - \int\_0^{\pi/2} dx \, \frac{d\hat{n}\_l(p)}{dp} \, \frac{\sin lx}{\cos^2 x} \right). \tag{41}$$

Then these integrals will be approximated by discrete sums in which we take for simplicity the same step *δ* for the variable *r* as well as for the variable *p*. Hence *nl*(*r*) � *nl*(*jδ*), with

$$n\_l(\mathbf{j}\delta) = \frac{1}{\pi} \left( \sum\_{k=1}^{j-1} \frac{d\widehat{n}\_l(p)}{dp} \int\_{\mathbf{x}\_j^k}^{\mathbf{x}\_j^{k+1}} d\mathbf{x} \, \frac{e^{-l\mathbf{x}}}{\cosh^2 \mathbf{x}} - \sum\_{k=j}^{K-1} \frac{d\widehat{n}\_l(p)}{dp} \int\_{\mathbf{x}\_j^k}^{\mathbf{x}\_j^{k+1}} d\mathbf{x} \, \frac{\sin l\mathbf{x}}{\cos^2 \mathbf{x}} \right), \tag{42}$$

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

$$I\_l(\mathbf{x}) = \int^\mathbf{x} dz \, \frac{e^{-lz}}{\cosh^2 z} \qquad \text{and} \qquad I\_l(\mathbf{x}) = \int^\mathbf{x} dz \, \frac{\sin lz}{\cos^2 z'} \tag{43}$$

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

$$(l-2)I\_l(\mathbf{x}) = 2\left(\tan\mathbf{x}\cdot\sin(l-2)\mathbf{x} - \cos(l-2)\mathbf{x}\right) - l\,I\_{l-2}(\mathbf{x}),\tag{44}$$

with *I*0(*x*) = 1/ cos *x* and *I*1(*x*) = −2 ln(cos *x*), and

$$\|I\|\_{l+1}(\mathbf{x}) = -\frac{1}{\mathbf{x}^{l}(1+\mathbf{x}^{2})} - (l+2)f\_{l-1}(\mathbf{x}),\tag{45}$$

with

18 Will-be-set-by-IN-TECH

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.

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

The idea is to start with the result of A. M. Cormack [15] (see his equation 18b) for the

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

We compute first *nl*(*r*) from a discretized version of equation 39, which shall be set up now,

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

*e*−*lx* cosh<sup>2</sup> *<sup>x</sup>* <sup>−</sup>

Then these integrals will be approximated by discrete sums in which we take for simplicity the same step *δ* for the variable *r* as well as for the variable *p*. Hence *nl*(*r*) � *nl*(*jδ*), with

> *dx <sup>e</sup>*−*lx* cosh<sup>2</sup> *<sup>x</sup>* <sup>−</sup>

�*l*

and *<sup>p</sup>* <sup>=</sup> *<sup>r</sup>*

� *π*/2 0

> *K*−1 ∑ *k*=*j*

� ∞ *r*

cos *x*

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

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

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

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

, (40)

�

*dx* sin *lx* cos2 *x*

sin *lx* cos2 *x*

� *<sup>x</sup>k*+<sup>1</sup> *j xk j*

⎞

. (41)

⎠ , (42)

⎞

⎟⎠ . (39)

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

*<sup>p</sup>*(*g*" <sup>+</sup> �*g*"<sup>2</sup> <sup>−</sup> <sup>1</sup>) �*g*"<sup>2</sup> <sup>−</sup> <sup>1</sup>

�*g*" <sup>+</sup> *<sup>g</sup>*"<sup>2</sup> <sup>−</sup> <sup>1</sup>) �*g*"<sup>2</sup> <sup>−</sup> <sup>1</sup>

*nl*(*p*(*g*" +

*nl*(*p*(*g*" +

�

�

*<sup>g</sup>*"<sup>2</sup> − <sup>1</sup>)). (37)

*<sup>g</sup>*"<sup>2</sup> − <sup>1</sup>)), (38)

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

� *g*" *τ* �

and *Nl*(*g*") = *<sup>p</sup>*(

by C. H. Chapman & P. W. Carey [12] for the classical Radon transform.

reconstructed angular Fourier component of *n*(*r*, *θ*) given by

�

and get *n*(*r*, *θ*) back by its angular Fourier series using equation 9.

*<sup>p</sup>* <sup>=</sup> *<sup>r</sup>* cosh *x*

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

> > � *<sup>x</sup>k*+<sup>1</sup> *j xk j*

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

�� <sup>∞</sup> 1

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

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

cos *l* cos−<sup>1</sup>

� <sup>1</sup> <sup>−</sup> *<sup>g</sup>*"2 *τ*2

*<sup>τ</sup> <sup>n</sup>*�*l*(*τ*) <sup>√</sup>*τ*<sup>2</sup> <sup>−</sup> <sup>1</sup> <sup>=</sup> <sup>2</sup> � *τ* 1 *dg*"

Now redefining the functions by

**5. Numerical simulations**

**5.1. Algorithm**

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

*π r*

⎛

⎜⎝ � *r* 0

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

*π*

⎛ ⎝ *j*−1 ∑ *k*=1

*nl*(*jδ*) = <sup>1</sup>

*π*

<sup>√</sup>*τ*<sup>2</sup> <sup>−</sup> <sup>1</sup>

*<sup>N</sup>*�*l*(*τ*) = *<sup>τ</sup> <sup>n</sup>*�*l*(*τ*)

$$f\_0(\mathbf{x}) = \frac{\mathbf{x}}{2(1+\mathbf{x}^2)} + \frac{1}{2}\tan^{-1}\mathbf{x} \quad \text{and} \quad f\_1(\mathbf{x}) = \ln\left(\frac{\mathbf{x}}{\sqrt{1+\mathbf{x}^2}}\right) + \frac{1}{2(1+\mathbf{x}^2)}.\tag{46}$$

The derivatives of *dn*�*l*(*p*)/*dp* are replaced by a simple linear interpolation *<sup>n</sup>*�� *<sup>l</sup>*(*k*), *i.e.*

$$\frac{d\widehat{n}\_l(p)}{dp} \simeq \frac{n\_l(k+1) - n\_l(k)}{\delta} = \widehat{n}'\_l(k). \tag{47}$$

Then the discretized value of the reconstructed angular component of *nl*(*r*), with *r* ∼ *jδ*, is

$$n\_l(j\delta) = \frac{1}{\pi} \left( \sum\_{k=1}^{j-1} \widehat{n}\_l^\flat(k) \left( I\_l((k+1)\delta) - I\_l(k\delta) \right) - \sum\_{k=j}^{K-1} \widehat{n}\_l^\flat(k) \left( I\_l((k+1)\delta) - I\_l(k\delta) \right) \right). \tag{48}$$

Equation 48 shall be used in numerical simulations3.

#### **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 **S** is placed below on the left of the image and the detector moves along the line *SD*. The relevant space is represented by 256 × 256 (length unit)2. Let *N<sup>φ</sup>* be the number of angular steps and let *Np* be the number of radii for circles going through the origin *O*. Then the following sampling steps *δφ* = 2*π*/*N<sup>φ</sup>* and *δ* = 4 256/*Np* 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 (*NMSE*) and the normalized mean absolute error (*NMAE*) (expressed as a percentage), defined by

$$NMSE = \frac{100}{N^2} \frac{\sum\_{(i,j) \in [1,N]^2} |\mathcal{Z}\_\mathbf{r}(i,j) - \mathcal{Z}\_\mathbf{o}(i,j)|^2}{\max\_{(i,j) \in [1,N]^2} \{\mathcal{Z}\_\mathbf{o}(i,j)\}^2} \quad \text{and} \quad NMAE = \frac{100}{N^2} \frac{|\mathcal{Z}\_\mathbf{r}(i,j) - \mathcal{Z}\_\mathbf{o}(i,j)|}{\max\_{(i,j) \in [1,N]^2} \{\mathcal{Z}\_\mathbf{o}(i,j)\}},$$

<sup>3</sup> For difficulties on the handling of the inversion formula of S. J. Norton, see [55]

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

where I*<sup>r</sup>* is the reconstructed image and I*<sup>o</sup>* 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 *pmax* 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 *pmax* leads to an increase of the detector length in the first CST modality.

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

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

[51]

17, see [43].

**Figure 15.** Shepp-Logan phantom reconstructed by the first CST modality (Norton 1994), reprinted from

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

As the set up works with radial parameters, an object of circular symmetry fits better to this

As for the CST third modality for small objects, the simulations results are displayed in Fig.

modality: this is the case of a medical phantom.

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

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.

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

20 Will-be-set-by-IN-TECH

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 *pmax* 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 *pmax* leads to an increase of the detector length

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

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

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.

objects, collected data and image reconstruction are given in Fig. 16, see [51].

where I*<sup>r</sup>* is the reconstructed image and I*<sup>o</sup>* is the original image.

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

**Figure 13.** Shepp-Logan original phantom, reprinted from [51]

CST modality is given in Fig. 15.

in the first CST modality.

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

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

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

(a) Data of the second CST modality, reprinted from [50]

(a) Data of the third CST modality, reprinted from [43]

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

(b) Image reconstruction by the third CST modality reprinted from [43]

Of course what is presented here does not claim neither completeness nor perfection. There

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

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) *SM*−1,−<sup>2</sup>

**Figure 17.** Third CST modality simulation results illustration

are still many important problems to tackle.

corrective algorithm, see *e.g.* [51].

**5.3. Open problems**

**Figure 16.** Second CST modality simulation results illustration

[50]

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.

(b) Image reconstruction by the third CST modality reprinted from [43]

**Figure 17.** Third CST modality simulation results illustration

#### **5.3. Open problems**

22 Will-be-set-by-IN-TECH

(a) Data of the second CST modality, reprinted from [50]

(b) Image reconstruction by the second CST modality, reprinted from

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

[50]

**Figure 16.** Second CST modality simulation results illustration

a good resolution, the loss of data is less important.

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) *SM*−1,−<sup>2</sup>

(or *MD*−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 (*SM* <sup>×</sup> *MD*)−1,−<sup>2</sup> is of the form of a product of a function of *<sup>τ</sup>* and a function of *<sup>r</sup>*, see [42, 54]. **Author details**

*F-95302 Cergy-Pontoise Cedex, France*

361-375. ISSN 0029-5450

173, pp. 64-67, ISSN 0766-5210

1209-1213. ISSN 0025-5327

0019-0578

502-507. ISSN 0167-5087, 0168-9002

Vol. 22, No. 2, pp. 229-244. ISSN 0031-9155

*research B*, Vol. 196, pp. 161-168. ISSN 0168-583X

*Problems*, Vol. 2, pp. 23-49. ISSN 0266-5611

*Laboratoire de Physique Théorique et Modélisation, University of Cergy-Pontoise/CNRS UMR 8089,*

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

*Laboratoire Équipes Traitement de l'Information et Systèmes, ETIS-ENSEA/University of*

[1] Adejumo O. O., Balogun F. A.& Egbedokun G. G. O., (2011). Developing a Compton Scattering Tomography System for Soil Studies: Theory, *Journal of Sustainable Development and Environmental Protection*, Vol. 1, No. 3, pp. 73-81. ISSN 2251-0605 [2] Anghaie S., Humphries L. L. & Diaz N. J., (1990). Material characterization and flaw detection, sizing, and location by the differential gamma scattering spectroscopy technique. Part1: Development of theoretical basis, *Nuclear Technology*, Vol. 91, pp.

[3] Arendtsz N. V. & Hussein E. M. A., (1993). Electron density tomography with Compton scattered radiation, *SPIE Mathematical Methods in Medical Imaging II* , Vol. 2035, Joseph

[4] Babot D., Le Floch C. & Peix G., (1992). Contrôle et caractérisation des matériaux composites par tomodensitométrie Compton, *Revue Pratique du Contrôle Industriel*, Vol.

[5] Balogun F. A. & Cruvinel P. E., (2003). Compton scattering tomography in soil compaction study, *Nuclear Instruments and Methods in Physics Research A*, Vol. 505, pp.

[6] Barrett H. H., (1984). The Radon Transform and its Applications, in *Progress in Optics*,

[7] Battista J. J., Santon L. W. & Bronskill M. J., (1977). Compton scatter imaging of transverse sections: corrections for multiple scatter and attenuation, *Phys. Med. Biol.*,

[8] Berodias M. G. & Peix M. G., (1988). Nondestructive measurement of density and effective atomic number by photon scattering, *Materials Evaluation*. Vol. 46, No. 8, pp.

[9] Bodette D. E. & Jacobs A. M., (1985). Tomographic two-phase flow distribution measurement using gamma-ray scattering, *ISA 1985*, Paper 85-0362, pp. 23-34. ISSN

[10] Brateman L., A. M. Jacobs & L. T. Fitzgerald, (1984). Compton scatter axial tomography

[11] Brunetti A., Cesareo R., Golosio B., Luciano P. & Ruggero A., (2002). Cork quality estimation by using Compton tomography, *Nuclear instruments and methods in Physics*

[12] Chapman C. H. & Carey P. W., (1986). The circular harmonic Radon transform, *Inverse*

with x-rays: SCAT-CAT, *Phys. Med. Biol*, Vol. 29, pp. 1353, ISSN 0031-9155

*Cergy-Pontoise/CNRS UMR 8051, F-95014 Cergy-Pontoise Cedex, France*

N. Wilson; David C. Wilson; Eds., pp. 230-41. ISSN 0924-9907

vol. 21, pp. 219-286, Ed. E. Wolf, North Holland. ISSN 0079-6638

Truong T. T.

Nguyen M. K.

**7. References**

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