Analytic inversion of a Radon transform on double circular arcs with applications in Compton Scattering Tomography

12/14/2019 ∙ by Cécilia Tarpau, et al. ∙ 1

In this work we introduce a new Radon transform which arises from a new modality of Compton Scattering Tomography (CST) which is made of a single detector rotating around a fixed source. Unlike some previous CST no collimator is used at the detector, such a system allows us to collect scattered photons coming from two opposite sides of the source-detector segment, hence the manifold of the new Radon transform is a family of double circular arcs. As first main theoretical result, an analytic inverse formula is established for this new Radon transform. This is achieved through the formulation of the transform in terms of circular harmonic expansion satisfying the consistency condition in Cormack's sense. Moreover, a fast and efficient numerical implementation via an alternative formulation based on Hilbert transform is carried out. Simulation results illustrate the theoretical feasibility of the new system. From a practical point of view, a collimator-free system considerably increases the amount of data, which is particularly significant in a scatter imaging system.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Since the seminal works of Radon [20] and Cormack [6] on the classical Radon transform on straight lines, several studies attempted to generalize this integral transform along various two and three dimensional manifolds. Among these generalizations, many concern circular paths, with integrals either on complete circles with Circular Radon Transforms (CRTs) or on parts of circles with Circular Arc Radon Transforms (CARTs). One can refer for example to the work of Cormack on the family of circles passing through the origin [7, 8], the work of Ambartsoumian on circles centered on a circle [1] for thermoacoustic and photoacoustic tomography or the work of Redding [21] about a circular Radon transform with applications in synthetic aperture radar imaging. About Radon transform on circular arcs, Nguyen and Truong [16, 30] proposed the inversion Radon transforms on different families of circular arcs modelling data acquisition and image reconstruction of new modalities of Compton Scattering Tomography (CST). Syed [26] proposed a numerical inversion for circular arcs with a fixed angular span with applications in photoacoustic tomography.

1.1 Basics of Compton Scattering Tomography (CST)

In Compton scattering tomography, the objective is to use radiation scattered by Compton effect to reconstruct the object. This imaging technique is an advance relatively to conventional Computed Tomography (CT), which only uses primary (non-deviated) radiation and considers scattered radiation as moise. Hence, the use of CST systems solves the radiation scattering problem and offers a better contrast in reconstructions [11, 22].

In fact, if a photon emitted with energy meets an electron inside matter, this photon is scattered off an angle with its original direction. Its energy after collision is given by the Compton formula

(1.1)

with the electron mass and the speed of light.

This one-to-one correspondence between energy and angle means that a detector measuring at is actually collecting the photons scattered on a circular arc. This circular arc passes through the source and the detector and subtends an angle . Consequently, assuming Compton effect is the only source of energy attenuation for emitted radiation and considering only first order scattering, modalities of Compton scatter tomography lead to CARTs on different families of circular arcs according to the chosen geometry of the system. The inversion of these Radon transforms represents the theoretical challenge raised by CST modalities, required for image reconstruction.

The first CST system was proposed by Norton [17]. This modality is made of a fixed source and a line of detectors passing through the source. Image acquisition is performed on half circular arcs having a fixed common end-point in the source and the other ends on the straight line of detectors. The reader can also refer to [23] where Rigaud gave another numerical inversion method for the Radon transform associated to Norton’s modality.

The second CST modality, proposed by Nguyen and Truong [16, 24, 23], is composed of a pair source - detector diametrically opposed on a circle. This system has lead to a CART based on circular arcs having a chord of fixed length.

A third modality has been proposed recently by Nguyen, Tarpau and Cebeiro [27, 29, 28, 3]. This modality is composed of a fixed source and a detector ring passing through the source. This one has led to a new Radon transform on circular arcs having a fixed common end-point in the source and the other ends on the detector ring. A similar system has been studied by Rigaud [25] and leads to the proposition of another method of reconstruction based on an approximate inversion and the use of mollifiers.

1.2 Motivation of this work: introduction of a new CST modality

The context of this study is the proposition of a new CST system. This modality is based on a fixed source and a detector moving on a circle of radius around the source (see Fig. 1), and localized by its angular position . Hence, Cartesian coordinates of are

(1.2)

The object to scan is placed outside of this circle.

Figure 1: Setup of the new CST modality
(In black) dotted circle: detector path
(In red) continuous curve: scanning double arcs, dotted curves: portion of circles not used for acquisition

A remarkable feature of this new simple system (with only one fixed source and one detector) consists in the absence of collimation on detectors. This characteristic enables an increase of the amount of acquired data for a given position of the detector and thus a reduction in the acquisition time. Thus, for an angular position of detector and given a scattering angle , two circular arcs correspond. This is why modelling of data acquisition with this collimation-free system leads to a Radon transform on double circular arcs (DCART).

1.3 Objectives and outline of the paper

Image reconstruction from projections obtained with this modality requires inversion of the Radon transform on double circular arcs (DCART). In this paper, we derive the analytic inversion of this integral transform.

This paper is outlined as follows. Section II introduces the modelling of image formulation corresponding to the forward DCART. Section III presents an analytic inversion formula through a formulation of the problem in terms of circular harmonic expansions of a function and its Radon transform on double arcs. Section IV deals with the establishment of consistency conditions in order to obtain a reconstruction formula based on Hilbert transform suitable for an efficient numerical implementation in Section V. Numerical simulations illustrating the theoretical feasibility of the new proposed system are presented in Section V. Section VI gives perspectives of this work.

2 Modelling of the collimation-free CST modality and the corresponding forward DCART

Figure 2: Parametrisation of double scanning circular arcs and detector path

Acquisitions using this sytem are made on a double family of circular arcs and of centre and whose union is denoted in the rest of the paper. This family of double circular arcs is then defined relatively to the diameter of scanning circles and the angular position of the detector (see Fig. 2). In fact, using or in order to define the family of double circular arcs is equivalent, according to the relation

(2.1)

Then, one can obtain

(2.2)

where polar equations of and are respectively

(2.3)

and

(2.4)

where .

This leads to a first formulation for the Radon transform on double circular arcs (DCART):

Definition 2.1 (Generalized Radon transform on double circular arcs associated to the collimation-free modality).

Denoting as an unknown function, non negative, continuous and compactly supported outside of the circle of radius , data acquisition with this new modality leads to the generalized Radon transform on the family of double circular arcs :

(2.5)

with a parameter associated to the considered double arcs defined in (2.2).

3 Analytic inversion of the Radon transform on double circular arcs

In this section, we give the detailed procedure for inverting the DCART. The procedure uses circular harmonic expansion to establish the forward transform in circular components (3.17), and then, with a similar approach as Cormack’s one, one obtains (3.23), which is the inverse formula in circular harmonic expansion.

3.1 Circular harmonic expansion of and

Functions and can be decomposed in terms of Fourier series, where and are respectively their circular harmonic expansion components :

(3.1)
(3.2)

where :

(3.3)
(3.4)

3.2 Connection between and

Data projection on double circular arcs can be decomposed on two Radon transforms on the families of circular arcs and . Consequently, we denote and the Radon transforms on respectively and and defined as follows:

(3.5)
(3.6)

with and the parameters respectively associated to and .

Hence, is the sum of and

(3.7)

One can decompose and in Fourier series to obtain and similarly as for .

Consequently, linearity of circular expansion yields

(3.8)

First step of this demonstration is to make explicit the circular expansion of and . Given that the reasoning is equivalent for both and , the result is only established for .

We use the expansion of in (3.1) to write as

(3.9)

Noticing the symmetrical property of about (see Fig. 3), (3.9) is rewritten as follows

(3.10)

with denotes the half part of circular arc , where .

Figure 3: Symmetry of about

Observing that

(3.11)

and plugging (3.11) in (3.10), one obtains circular harmonic expansion of

(3.12)

Straightforward computations show that

(3.13)

and

(3.14)

Equation (3.12) becomes

(3.15)

Similarly, we get

(3.16)

Hence, from the addition of (3.15) and (3.16), the connection between circular components of and can be written

(3.17)

Equation (3.17) is the equation of image formation in the circular harmonic expansion.

3.3 On inversion of

Denoting

(3.18)

and multipliying both sides of (3.17) by

one gets

(3.19)

Rearranging (3.19) gives

(3.20)

In [7], it has been proven that

(3.21)

Hence, substituting (3.21) in (3.20) and differentiating with respect to the variable , one gets

(3.22)

Going back to the variable and with the change of variables , one finds

(3.23)

Equation (3.23) is the equation of image reconstruction in the circular harmonic expansion. Furthermore, one can observe that (3.23) illustrates the Cormack’s hole theorem: in order to determine by its circular harmonic expansion , the knowledge of projections of in the annular domain is sufficient.

4 Consistency conditions

We introduce consistency conditions [8, 30] in terms of Cormack sense in order to deduce a closed formulation of (3.23), more suitable for numerical computation.

Equation (3.23) can be rewritten using the n-th order Tchebychev polynomial of first kind as follows :

(4.1)

Using the following relationship between Tchebychev polynomials of first kind and second kind [8]

(4.2)

and according to (3.18), (4.1) becomes :

(4.3)

Furthermore, from (3.17), one can obtain

for and

(4.4)

The order of integration can be changed in (4.4) to obtain:

(4.5)

With the change of variables , the previous -integral becomes

(4.6)

We arrive to an integral whose result is given in [13]

(4.7)

where refers to the gamma function.

For and denoting , one can remark that are poles of . Since is pair, then for are also poles of .

Consequently, one can conclude that this integral vanishes [30]:

For and ,

(4.8)

Furthermore, since Tchebychev polynomials of second kind is a linear sum of polynomials (see [13]), (4.8) leads to :

For ,

(4.9)

5 A closed-form inverse formula

This section uses consistency conditions established in (4.9) to obtain the final closed-form inverse formula (5.9).

From (4.9), one can deduce

(5.1)

Consequently, using (5.1) and the relation between Tchebychev polynomials of first and second kind (4.2), one can obtain

(5.2)

Then, we denote the projections of circular harmonic components

(5.3)
(5.4)

Using (3.1) and (5.4) in (5.2), one can obtain

(5.5)

We need the two following lemmas to simplify (5.5).

Lemma 5.1 ([30], Lemma 3.8).

For , one has

(5.6)
Lemma 5.2 ([30], Lemma 3.9).

For , we have

(5.7)

Hence, one can obtain :

(5.8)

where refers to the Cauchy principal value.

Equation (5.8) is a first closed formulation for image reconstruction with in polar coordinates. Moreover, an alternative formulation is obtained shoving the derivative inside the -integral. This goal is achieved using substitution and the fact that for . Then, we go back to the original variables to obtain:

(5.9)

Equation (5.9) is the final equation of image reconstruction in polar coordinates.

The intermediate result (5.2) is also an inversion formula based on circular harmonic expansion. However, the proposed formulation (5.9) is more attractive for numerical reconstruction than (5.2). In fact, (5.2) can be numerically implemented using a method similar to Chapman and Cary numerical approach [5], see for instance [23]. Nevertheless, the technique requires the evaluation of Tchebychev functions in (5.2) in terms of primitives integrals that are evaluated recursively. This implies either higher computational time or memory requirements. On the other hand, (5.9) can be rewritten using the Hilbert transform and implemented in a more efficient way using standard tools of discrete Fourier analysis, as we are going to show in the next section.

6 Reconstruction algorithm and numerical results

6.1 Parameter choices

The size of the object (Shepp-Logan phantom, see Fig. 4) is pixels. The detector is moving on a ring of radius pixels with an angular step of degree. Then, the detector will have an amount of different positions to collect data. For the number of double scanning circular arcs for each detector position, we choose . In fact, the parameter choices should satisfy the condition according to [2]. This condition is the lower limitation for a data acquisition in an ideal case, that is with an uniform sampling and the absence of perturbations such as noise or missing data.

In our case, we have to overestimate because of the non uniform sampling due to the acquisition on the family of double arcs. Furthermore, we face also to a missing data problem because in the ideal case, scanning double arcs with a diameter have to be considered to have the whole acquired data. On simulation, we choose a maximal diameter , in a real setting this parameter is linked to the energy resolution of the detector.

Figure 4: Original object: Shepp-Logan phantom

6.2 Image formation

6.2.1 General algorithm

A parametrisation in Cartesian coordinates (instead of equations in polar coordinates) is preferable to perform numerical simulations in order to have the same distance between adjacent running points on the considered scanning circle.

Figure 5: Cartesian parametrisation of and

Hence, a scanning arc can be seen as the shift of a circle centred at the origin of identical radius, to its center with a restriction of the domain of the variable :

(6.1)

Consequently, Cartesian parametrisation of and are respectively

(6.2)

where

and

(6.3)

where

Then, image acquisition on double arcs is modelled by the equation

(6.4)

6.2.2 Simulation results

Figure 6 shows the result of image formation for the Shepp-Logan phantom.

Figure 6: Image formation of the original object (Fig. 4)

6.3 Image reconstruction

6.3.1 General algorithm

In order to reconstruct the function, we use the Hilbert transform which is related to the Cauchy principal value as follows:

(6.5)

In fact, the Hilbert transform is easily computed in the Fourier domain, where we have :

(6.6)

where

denotes the one-dimensional Fourier transform.

Equation (5.9) of image reconstruction becomes consequently

(6.7)

Finally, using the correspondence between polar coordinates and Cartesian coordinates

(6.8)

image reconstruction equation used for simulation in Cartesian coordinates is

(6.9)

The projections are computed via the circular harmonic components of with (5.4). However, zeros in the denominator may be source of instability and regularization may be required. According to [15], we add a regularization parameter in (3.18) (equal to in the proposed simulation) in order to compute the circular harmonic components of :

(6.10)

Algorithm 1 summarizes the different steps for reconstructing the object from (6.9).

Data: , projections on double circular arcs of function
Result:
1 Compute circular harmonic expansion of to compute with (6.10) and recompose ;
2 Compute discrete derivation of relative to variable and multiply the result by ;
3 Write the Hilbert transform as a filtering operation in Fourier domain using (6.6);
4 For each

, interpolate the data on the considered scanning circles

of (6.9);
5 Weight the result using the factor ;
6 Sum the weighted interpolations on all direction ;
7 Weight the result by ;
Algorithm 1 Reconstruction of object

6.3.2 Evaluation of quality of reconstruction

In order to evaluate the quality of the reconstructions, we use the Normalized Mean Square Error (NMSE) defined by:

(6.11)

where and are respectively the original and the reconstructed images.

6.3.3 Simulation results and remarks

Fig. 7 shows the result of image reconstruction. It can be noted that the object is well reconstructed as a whole. In particular, small details of the object are well visible at reconstruction.

Figure 7: Image reconstruction - NMSE = 0.01

7 Additional considerations and perspectives

Some interesting issues arise when studying the new Compton scattering modality.

7.1 Consideration of attenuation in the model

In a practical setting, attenuation may be considered in the forward model. This leads to a major challenge from the mathematical point of view for which no exact solution is known at present. Nevertheless, a number of ways to deal with the problem of attenuation in Compton Scattering Tomography have been studied like the generalized Chang Correction (GCC) [4], an iterative algorithm that allows the a posteriori correction of the reconstruction or the Iterative Pre Correction algorithm [14] which corrects data. This later method, which has been adapted for the previous Nguyen and Truong’s modality in [10], requires having the analytic inverse Radon transform found in the non attenuated case, i.e. the problem that we have efficiently solved here.

7.2 On artifacts and image quality

Image reconstruction (see Fig. 7) has artifacts, located on the lower left edge and upper right edge of the reconstruction. These artifacts are due to missing data. In fact, projections used for reconstruction are limited by circular arcs having a finite diameter . Some interesting works about the consequences of data restriction have been proposed by Nosmas [18], Frikel and Quinto [9, 19, 12] for the classical Radon transform on straight lines. Some filtering strategies arising in these works may be suitable for image improvement in our DCART.

7.3 A proposed extension of this modality in three dimensions

In addition, an extension of this modality to three dimensions can be also proposed: a fixed source and a detector moving on a sphere centered at the source. This modality leads to a new toric Radon transform, resulting of the rotation of the circular arcs around the source-detector axis. The theoretical challenge raised by this extension consists in the inversion of this new toric Radon transform. Some work on this direction is on the way.

8 Concluding remarks

A new modality of Compton Scattering Tomography has been studied in this paper. The design has attractive features such as compactness and simplicity. Furthermore, this CST modality is a collimation-free modality and this allows to increase data and simultaneously reduce acquisition time. This modality leads to a new Radon transform on double circular arcs, for which a closed analytic inversion formula is established in this paper.

In addition, an efficient reconstruction algorithm based on the Hilbert transform has been developed.

Simulation results illustrate the feasibility of this new system and the good performance of the reconstrution algorithm.

An efficient method based on Hilbert transform may enable reconstruction in a 3D version of the design where the amount of data turns critical. However, more research work must be done in order to develop a 3D extension.

9 Funding and acknowledgements

C. Tarpau research work is supported by grants from Région Ile-de-France (in Mathematics and Innovation) 2018-2021 and LabEx MME-DII (Modèles Mathématiques et Economiques de la Dynamique, de l’Incertitude et des Interactions) (No. ANR-11-LBX-0023-01).

J. Cebeiro and M. A. Morvidone are suported by SOARD-AFOSR (grant number FA9550-18-1-0523). J. Cebeiro is supported a CONICET postdoctoral grant ( 171800).

The authors would like to thank also the Institute for Advanced Studies (IAS) of the University of Cergy Pontoise, France for the financial support given to this research project.

References

  • [1] G. Ambartsoumian, R. Gouia-Zarrad, and M. A. Lewis (2010) Inversion of the circular Radon transform on an annulus. Inverse Problems 26 (10), pp. 105015. Cited by: §1.
  • [2] R. N. Bracewell (1990) Numerical transforms. Science 248 (II May), pp. 697–704. Cited by: §6.1.
  • [3] J. Cebeiro, M. K. Nguyen, M. Morvidone, and C. Tarpau (2018-11) An interior Compton Scatter Tomography. In 25th IEEE Nuclear Science Symposium and Medical Imaging Conference 2018 (IEEE NSS/MIC’18), Sydney, Australia. Cited by: §1.1.
  • [4] L. Chang (1978) A method for attenuation correction in radionuclide computed tomography. IEEE Transactions on Nuclear Science 25 (1), pp. 638–643. Cited by: §7.1.
  • [5] C. Chapman and P. Cary (1986) The circular harmonic Radon transform. Inverse Problems 2 (1), pp. 23. Cited by: §5.
  • [6] A. M. Cormack (1963) Representation of a function by its line integrals, with some radiological applications. Journal of Applied Physics 34 (9), pp. 2722–2727. Cited by: §1.
  • [7] A. M. Cormack (1981) The Radon transform on a family of curves in the plane. Proceedings of the American Mathematical Society 83 (2), pp. 325–330. Cited by: §1, §3.3.
  • [8] A. M. Cormack (1984) Radon’s problem–old and new. In SIAM-AMS Proceedings, Vol. 14, pp. 33–39. Cited by: §1, §4, §4.
  • [9] J. Frikel and E. T. Quinto (2015) Artifacts in incomplete data tomography with applications to photoacoustic tomography and sonar. SIAM Journal on Applied Mathematics 75 (2), pp. 703–725. Cited by: §7.2.
  • [10] O. A. O. Guerrero, G. Rigaud, R. Régnier, and M. K. Nguyen (2013) Attenuation correction in a new modality of Compton Scattering Tomography. Cited by: §7.1.
  • [11] K. C. Jones, G. Redler, A. Templeton, D. Bernard, J. V. Turian, and J. C. Chu (2018) Characterization of Compton-scatter imaging with an analytical simulation method. Physics in Medicine & Biology 63 (2), pp. 025016. Cited by: §1.1.
  • [12] V. P. Krishnan and E. T. Quinto (2014) Microlocal analysis in tomography. Handbook of mathematical methods in imaging, pp. 1–50. Cited by: §7.2.
  • [13] W. Magnus, F. Oberhettinger, and R. P. Soni (2013) Formulas and theorems for the special functions of mathematical physics. Vol. 52, Springer Science & Business Media. Cited by: §4, §4.
  • [14] A. Maze, R. Collorec, P. Bourguet, and M. Pierpittte (1992) Iterative reconstruction methods for non uniform attenuation distribution in SPECT. In 1992 14th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, Vol. 5, pp. 1821–1822. Cited by: §7.1.
  • [15] S. Moon and M. Haltmeier (2017) Analytic inversion of a conical Radon transform arising in application of Compton cameras on the cylinder. SIAM Journal on imaging sciences 10 (2), pp. 535–557. Cited by: §6.3.1.
  • [16] M. K. Nguyen and T. T. Truong (2010) Inversion of a new circular-arc Radon transform for Compton scattering tomography. Inverse Problems 26 (6), pp. 065005. Cited by: §1.1, §1.
  • [17] S. J. Norton (1994) Compton scattering tomography. Journal of applied physics 76 (4), pp. 2007–2015. Cited by: §1.1.
  • [18] J. Nosmas (1996) Analyse Microlocale et Tomographie Géométrique. In ICAOS’96, pp. 338–344. Cited by: §7.2.
  • [19] E. T. Quinto (2017) Artifacts and visible singularities in limited data X-ray tomography. Sensing and Imaging 18 (1), pp. 9. Cited by: §7.2.
  • [20] J. Radon (1917) Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten. Akad. Wiss. 69, pp. 262–277. Cited by: §1.
  • [21] N. J. Redding and G. N. Newsam (2001) Inverting the circular Radon transform. Technical report Defense Science and Technology Organisation Salibury (Australia). Cited by: §1.
  • [22] G. Redler, K. C. Jones, A. Templeton, D. Bernard, J. Turian, and J. C. Chu (2018) Compton scatter imaging: a promising modality for image guidance in lung stereotactic body radiation therapy. Medical physics 45 (3), pp. 1233–1240. Cited by: §1.1.
  • [23] G. Rigaud, M. K. Nguyen, and A. K. Louis (2012) Novel numerical inversions of two circular-arc Radon transforms in Compton scattering tomography. Inverse Problems in Science and Engineering 20 (6), pp. 809–839. Cited by: §1.1, §1.1, §5.
  • [24] G. Rigaud, R. Régnier, M. K. Nguyen, and H. Zaidi (2013) Combined modalities of Compton scattering tomography. IEEE Transactions on Nuclear Science 60 (3), pp. 1570–1577. Cited by: §1.1.
  • [25] Gaël. Rigaud (2017) Compton Scattering Tomography: Feature Reconstruction and Rotation-Free Modality. SIAM Journal on Imaging Sciences 10 (4), pp. 2217–2249. External Links: Document, Link, https://doi.org/10.1137/17M1120105 Cited by: §1.1.
  • [26] T. A. Syed, V. P. Krishnan, and J. Sivaswamy (2016) Numerical Inversion of Circular arc Radon Transform. IEEE Transactions on Computational Imaging 2 (4), pp. 540–549. Cited by: §1.
  • [27] C. Tarpau, J. Cebeiro, M. Morvidone, and M. K. Nguyen (2019) A new concept of Compton Scattering tomography and the development of the corresponding circular Radon transform. IEEE Transactions on Radiation and Plasma Medical Sciences (accepted for publication). Note: [10.1109/TRPMS.2019.2943555] Cited by: §1.1.
  • [28] C. Tarpau and M. K. Nguyen (2018-07) A novel modality of Compton Scattering Tomography, Image formation and Reconstruction. In

    Proceedings of the 22nd International Conference on Image Processing, Computer Vision, and Pattern Recognition (IPCV’18)

    ,
    Las Vegas, United States, pp. 60–65. Cited by: §1.1.
  • [29] C. Tarpau and M. K. Nguyen (2019) Scattering imaging system with dual configuration. In 14th International Conference on Quality Control by Artificial Vision, Vol. 11172, pp. 111720Y. Cited by: §1.1.
  • [30] T. T. Truong and M. K. Nguyen (2011) Radon transforms on generalized Cormack’s curves and a new Compton scatter tomography modality. Inverse Problems 27 (12), pp. 125001. Cited by: §1, §4, §4, Lemma 5.1, Lemma 5.2.