1 Introduction
Compton Scattering Tomography (CST) is an imaging technique whose objective is to exploit wisely Compton scattered photons by the object to scan in order to reconstruct its electron density map. Since the early proposition of this type of imaging by Lale [16], Clarke [7] and Farmer [10], studies about CST systems proved already promising results in medical imaging, for instance for the identification of lung tumours [22, 15], but also for earthquake engineering [11], cultural heritage imaging [20, 13], landmine detection [14] and agricultural measurements [9]. In fact, CST has made it possible to widen the fields of application for tomography to the imaging of onesided large objects, because, in some configurations, sources and detectors can be placed at the same side of the object.
The proposition of CST systems is strongly related to the study of the associated integral transform, which models data acquisition [31, 18]. These integral transforms are generalizations of the classical Radon transform on lines studied by Radon [21] and Cormack [8]. While, in twodimensional CST, manifolds are on circle arcs or double circle arcs, in three dimensions, data is acquired on toric surfaces. These circular geometries for the considered manifolds originate from the Compton effect. This physical phenomenon occurs when a photon emitted by the source with energy collides with an electron as it passes through matter. This photon is scattered, and deviated by an angle from its original direction. The photon looses also a part of its energy. The Compton formula gives us the onetoone correspondence between the energy of scattered photons and the related scattering angle
(1.1) 
is the electron mass and the speed of light. This relation ensures also that scattering sites of photons with identical energy are located on a circle arc labelled by the scattering angle . Figure 1 illustrates the general functioning principle of a CST system.
Several issues around the study of these generalized Radon transforms are then of interest, such as existence and uniqueness of the solution, its stability or the range conditions. The main important problem remains the reconstruction of the image, mostly coming from the proposition of analytical inverse formulas.
In that way, the first studied twodimensional CST system was made of a fixed line of detectors containing a fixed source [19]. See Figure (a)a. The associated scanning circle arcs consist of a family of semicircles with a common fixed end point (the point source) and another extremity (the considered detector) on the line. Then, circular geometries for CST have been proposed, the first one considers a pair source  detector diametrically opposed in rotation around the object [17, 24, 25] and data acquisition is modelled by a Radon transform on circle arcs having a fixed sourcedetector chord length. See Figure (b)b. The second consists also of a pair source  detector moving on a circle, but the distance between the two is no longer constant and depends on the scattering angle [30] (Fig. (c)c). With this modality, data measurement is performed on a family of circles orthogonal to the fixed circular path of the pair source  detector. In a third modality (see Fig. (d)d), on the contrary, it was considered a set of detectors placed on a ring containing the source, thus obtaining a completely fixed CST modality [28, 29, 5, 26]. The corresponding Radon transform measures the contribution of the photons scattered at the points located on circles arcs having a common extremity (the point source) and another (the considered detector) on the detector ring. Some configurations, mentioned above, employ collimators to split up photons coming from different circle arcs. Geometries without collimation have also been studied. Consequently, for a given position of detector, the acquisition is performed on a family of double circle arcs and the amount of registered data is potentially doubled. This feature may also reduce the acquisition time and finally radiation exposition in comparison with a similar geometry with collimation. Among the proposed CST systems without collimation, one supposes a fixed source and a detector rotating around this source [27]. See Figure (e)e. The other one, proposed in [34], is made of a source and a detector which translate along a line (Fig. (f)f).
The direct reconstruction of volumes is also of interest with the proposition of threedimensional CST systems. These modalities use uncollimated sources and detectors, and data acquisition is performed on toric surfaces. In many cases, these 3D modalities correspond to extensions of 2D systems. Circular geometries become thus spherical or cylindrical [36, 23, 6] and linear geometries become planar [34].
Here, we are interested in the twodimensional modality proposed in [34]
. The purpose of this article is to develop an alternative formulation suitable for the development of a faster and efficient reconstruction algorithm. The associated reconstruction algorithm will use only classical tools such as Fast Fourier Transform (FFT) algorithm.
The paper is outlined as follows. Section 2 recalls the general setup of the system and the model for data acquisition for the modality. Section 3 introduces the main result of the paper, that is, the alternative formulation for the inverse Radon transform on double circle arcs. Section 4 will give the discrete formulation of the forward operator, as well as the proposed strategy to reconstruct the object under study. Section 5 discusses the obtained simulation results with a study of the influence of some parameters on reconstruction quality.
2 Setup and measurement model of the CST system under study
2.1 Setup
The system under study is made of a source, assumed to be monochromatic, and a detector separated by a fixed distance from each other. The source and the detector move respectively on a horizontal line of equation and . The horizontal position of the pair sourcedetector is labelled by (see Figure 3). Alternatively, this system may be sketched with fixed lines of sources and detectors that will be used in pair. The object, placed below the detector path, is scanned transversely. No collimation is used at the detector, hence the acquisition is performed on a family of double circle arcs (called toric sections in the original publication [34])^{1}^{1}1We position ourselves in the same frame of study as the original article, with the same working assumptions. Thus, first order Compton scattering is the only source of attenuation for radiation and data acquisition is performed with a pair source detector, assumed to be pointlike. These conditions are common in the literature [17, 36, 35, 4, 33, 28, 27] and have been already discussed in [19, 30, 36, 27].. We parameterize these circle arcs and define the corresponding Radon transform in the next paragraph.
2.2 Modelling of data acquisition using the CST system
Given a scattering angle , data acquisition is performed on a family of double circle arcs of radius where (or, equivalently, ). For parameterization, these double circle arcs are obtained with the union of four half arcs denoted , of respective equation
and . See Figure 3.
The Radon transform which mathematically models data measurement with this CST system is then defined as follows :
Definition 1.
Let be an unknown function, nonnegative, continuous and compactly supported in the half plane . The Radon transform on double circle arcs maps into the set of its integrals over the family of double circle arcs as
(2.1) 
where refers to the elementary arc length measure on the considered double circle arc. Then, after computation of the arc length measure, we have the explicit reformulation for [34, Proposition 3.1] :
(2.2) 
where .
In the original study of this modality, the invertibility of the corresponding Radon transform as well as its analytical inversion formula has been established. The invertibility was proven using the theory of integral equations and resulted in a Volterra integral equation with a weakly singular kernel in the Fourier domain. This study also leads to a formulation for inversion formula as an integral transformation with a kernel computed iteratively. The numerical calculation of this kind of kernel may require high computational time and/or memory. Furthermore, as mentioned in the Remark 3.4 of the original paper, the proposed approach by Webber is severely illposed, particularly in terms of stability. Implementing such a method can lead to large instabilities, even when these are due to small changes in the data.
3 An alternative formulation for the inversion formula of the Radon transform on double circle arcs
In this section, we state the main result of the paper, a different formulation for the associated inversion formula, that will be easier to implement numerically. Let us introduce before some notations that will be used in the proofs.
3.1 Notations
It is useful to define the following transform pairs.
Definition 2 (Fourier transform).
Let be a compactly supported function in . The ndimensional Fourier transform of , denoted , is given by
(3.1) 
with . The inverse Fourier transform is
(3.2) 
Definition 3 (Fourier cosine transform [12]).
Let be a compactly supported function in . The Fourier cosine transform of , denoted , is given by
(3.3) 
with . The inverse Fourier cosine transform is
(3.4) 
We define also the Hankel transform.
Definition 4 (Hankel transform [3]).
Let be a compactly supported function in . The zeroorder Hankel transform of is defined as
(3.5) 
where stands for the Bessel function of the first kind of order .
Finally, we recall the integral representation of the Bessel function :
(3.6) 
3.2 Inversion formula
Proposition 1.
Denoting the operator whose Fourier transform according to the first variable is
(3.7) 
if and when , the unknown function is completely recovered from as follows
(3.8) 
Proof.
With the change of variables in (2.2) and taking the Fourier transform of respectively to variable , one gets
(3.9) 
With the second change of variables , one gets
(3.10) 
where stands for the onedimensional Fourier transform relatively to the variable .
Using relation (3.7), multiplying both sides of (3.10) by with and integrating with respect to variable , for , one recognizes the Hankel transform of , denoted
(3.11) 
Then, with the double substitution , one gets
(3.12) 
The result of the integral is given in the table [2, p. 55, eq. (35)]. Finally, one gets for
(3.13) 
and if .
Let . Then with the fact for ,
(3.14) 
The lefthand side is the Fourier cosine transform of according the variable . We can then extract , applying the inverse cosine transform to (3.14)
(3.15) 
where . The final equation is obtained going back to variable , and applying the inverse Fourier transform. ∎
Remark 1.
The projections (3.7) contain zeros in the denominator, since the cosine function vanishes when , . From (3.11) to the end of the demonstration, it was supposed that is different from . Furthermore, this may be a source of instability in the simulations. The addition of a regularization parameter for simulations is discussed in Section 4.2 to prevent this.
Remark 2.
Another reconstruction algorithm is also possible from the projections of . This process is achieved performing geometric inversion. Geometric inversion is a mapping converting a point into a point such that , where is a constant value. The mapped point has the same direction as the original point but a distance of to the origin of the considered coordinate system. As an example, geometric inversion converts circles passing through the origin into straight lines. In the present case, the Radon transform on double circle arcs is converted into a Radon transform on an apparent family of circle arcs of similar geometry as the one studied in [30, 32]. Although the inverse problem can be alternatively solved using geometric inversion, the approach we employ here is more straightforward.
4 Numerical formulations for the forward and inverse transform
4.1 Image formation
Let be the number of positions for the pair source  detector and the number of double scanning circle arcs per sensor position. We denote , and , the discrete variables corresponding respectively to and . The matrix of projection data is then computed, writing (2.2) under a discrete form, with the change of variables
(4.1) 
where is the sampling angular distance of . The above Cartesian parameterization allows having a constant distance between running points of the considered scanning circle arcs during simulations.
4.2 Image reconstruction
For image reconstruction, we need to compute the projections in the Fourier domain according to (3.7). This expression contains zeros in the denominator. This may induce instabilities on reconstruction. For simulations, we add a small regularization parameter denoted
(4.2) 
In terms of computational cost, the most demanding step in the implementation of (3.8) is the calculation of . The idea is to establish a relation between the above operator with the Fourier transform of , defined as follows
(4.3) 
We have now the following proposition.
Proposition 2.
Let . is related with the twodimensional Fourier transform of the operator as
(4.4) 
Consequently, from the inversion formula (3.8), it follows in the Fourier domain
(4.5) 
where is the twodimensional Fourier transform of .
Proof.
From the definitions of Fourier and Hankel transforms and with the integral representation of the Bessel function, one gets
(4.6) 
There is an angle which corresponds to the angular coordinate of the point , such that and . Using the property of periodicity of trigonometric functions, it follows that
(4.7) 
Changing variables and ,
(4.8) 
This leads to the reconstruction algorithm summed up in Algorithm 1.
5 Simulations results
The original object used for simulations is Derenzo phantom, an object made of multiple circles of different sizes. The circles in the object also allow the study of the performance of the algorithm in front of different contrasts and spatial resolution, as well as its ability to reconstruct features locally tangent to lines of any slope. The unit length used here is the pixel. We suppose thus that the distance between the source and the detector paths of the modality is two pixels. The size of the object is pixels in all simulations. Furthermore, given the linear geometry of this modality, there is no loss of generality to consider the object centred relatively to the axis.
5.1 Data acquisition
Data measurement is calculated according to (4.1). Figure 4 shows an example of data obtained for Derenzo phantom.
5.2 Influence of some parameters on reconstruction quality
In the following paragraphs, we study the influence of the other general parameters of the system such as the position of the object, the number of required positions for sensors or the number of scanning circle arcs. In addition to a visual comparison of the reconstruction quality, we propose here to measure quantitatively the error rate between the original object and the reconstruction with the Normalized Mean Squared Error (NMSE) where refers to the norm.
The regularization parameter , which has to be small, was arbitrarily set to .
5.2.1 Position of the object relative to the detector path
We analysed here the influence of the position of the object on reconstruction quality. Firstly, the number of positions for the pair source  detector and the number of scanning circles per position was chosen to largely satisfy the wellknown condition [3] and are set arbitrarily to and . A convenient choice for these parameters will be discussed later. We performed various acquisition, modifying the gap between the detector path and the upper part of the object (see Figure 3). Consequently, the object is in the square of Cartesian coordinates . Figure 5 shows the reconstruction results for and pixels which correspond respectively a position for the object in the respective domains , and along the zaxis. The difference between the three reconstructions is on the top of the object. If the object is close to the line of movement of the detector, then this part is less well reconstructed. Indeed, this distance between the object and the detector path allows having arcs of circle tangent horizontally to this part of the object. An offset of seems to be a good tradeoff between quality of reconstruction and the applicability of such a measure in practical use. For the rest of the simulations, is set to .
5.2.2 Number of necessary positions for the pair sourcedetector
We studied then the number of different positions required for a good quality of reconstruction. The influence of two running parameters is analysed, first, the farthest position from the object for the sourcedetector pair (that is an array of length for the source and detector paths) and the distance between two adjacent positions of the pair.
Figures (a)a, (b)b and (c)c show the reconstruction results when the farthest position from the object to the pair source detector is respectively and with a common set to . Reconstruction from a domain appears to be blurred with strong artefacts in the upper parts of the image. For and , reconstruction quality seems to be visually equivalent, even if the NMSE for is higher. This may be due to numerical approximations. For the rest of the simulation, is set to .
The influence of the distance between two adjacent positions of the pair is now evaluated. Figures (d)d, (e)e and (f)f show the result for and pixels. This represents a respective amount of and detectors per unit length. In Figure (d)d , we can see streaks suggesting a lack of data for reconstructing the object. On the contrary, the doubling of the number of detectors between Fig. (e)e and Fig. (f)f does not bring a better quality of reconstruction. Consequently, the use of one detector per unit length seems to be a good tradeoff, and the value will remain constant in the rest of the paper.
5.2.3 Number of scanning circles per position of the pair sourcedetector
The number of scanning circles necessary for reconstruction is now under study. In the same way, two parameters are of interest, that is, the maximum radius of the scanning double circle arcs to be taken and the discretization step that have to be chosen. We first evaluate the consequences of the value of with three examples on Fig. (a)a, (b)b and (c)c where is set respectively to and and . For , the reconstruction suggests a lack of data in front of the obtained results for and . Notice the higher NMSE for
, probably due to numerical approximations.
5.2.4 Discussions
The above simulations results show some interesting issues for the CST modality and the reconstruction quality that can be expected with such a system. First, this system is able to scan objects whose depth is largely oversized relative to the sourcedetector distance. However, this seems to be counterbalanced by the need for a consequent length for the source and detector linear paths, since sufficient reconstruction results appear for source and detector paths of length size six times greater than the object size. The necessary number of scanning circle arcs is also very important since, if we relate the values of with the scattering angles, the reconstruction quality is largely improved when is high whereas the angular distance between two values of is very small. Considering a larger distance between the source and the detector paths for the system may reduce partially the required amount of projection data.
Moreover, the geometry offers a sufficient reconstruction quality for every tangent with arbitrary slopes. However, one can notice that vertical slopes are slightly less well reconstructed than the other ones. This can be seen if we pay carefully attention to the left and right sides of the reconstruction. This may be problematic if the object to scan is essentially made of vertical features. One way to avoid this is to perform additional scans by rotating the object, if possible.
Some artefacts that look like shadows around the circles which compose the object remains clearly visible. The issue of artefacts in CST has already been addressed in different manners, for instance in [38], where microlocal analysis was employed to alleviate artefacts in reconstructions from limited data. Moreover, in [37]
, a penalized iterative algorithm was developed for a mixed modality. Regarding our approach, it can be combined in a pipeline with postprocessing stages based on machine learning, as we did in
[1] for limited data issues in classical computed tomography. Some work is on the way.6 Concluding remarks
In this article, we proposed a new reconstruction algorithm for a recently proposed CST modality with translational geometry. The algorithm can be numerically implemented efficiently and reconstructions exhibit good quality. A quantitative study of the required data has also been carried out. This study proved the ability for such a CST system to reconstruct larger objects than the system itself. This advantage is moderated by the fact that it is necessary to have an important amount of sourcedetector positions, which can make the acquisition time longer. An interesting issue of this work concerns a reconstruction algorithm for the threedimensional extension of this system with planar paths for the source and the detector.
7 Acknowledgements
We would like to thank Prof. T. T. Truong for stimulating discussions.
C. Tarpau research work is supported by grants from Région ÎledeFrance (in Mathematics and Innovation) 20182021 and LabEx MMEDII (Modèles Mathématiques et Économiques de la Dynamique, de l’Incertitude et des Interactions) (No. ANR11LBX002301).
J. Cebeiro research work is supported by a postdoctoral grant from the University of San Martín. He is also partially supported by SOARDAFOSR (grant number FA95501810523).
References

[1]
(202107)
Deep morphological networkbased artifact suppression for limitedangle tomography.
In
Proceedings of the 25th International Conference on Image Processing, Computer Vision and Pattern Recognition (IPCV’21)
, Las Vegas, United States. Cited by: §5.2.4.  [2] (1954) Tables of integral transforms. Vol. 1, McGrawHill Book Compagny, New York. External Links: ISBN 070195498 Cited by: §3.2.
 [3] (1990) Numerical transforms. Science 248 (II May), pp. 697–704. Cited by: §5.2.1, Definition 4.
 [4] (2017) New “improved” Compton scatter tomography modality for investigative imaging of onesided large objects. Inverse Problems in Science and Engineering 25 (11), pp. 1676–1696. Cited by: footnote 1.
 [5] (201811) 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.
 [6] (2021) On a three dimensional Compton scattering tomography system with fixed source. Inverse Problems 37 (5), pp. 054001. Cited by: §1.
 [7] (1969) Comptonscattered gamma rays in diagnostic radiography. In Medical Radioisotope Scintigraphy. VI Proceedings of a Symposium on Medical Radioisotope Scintigraphy, Cited by: §1.
 [8] (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.
 [9] (2006) Compton scattering tomography for agricultural measurements. Engenharia Agricola 26 (1), pp. 151–160. Cited by: §1.
 [10] (1971) A new approach to the determination of anatomical crosssections of the body by compton scattering of gammarays. Physics in Medicine & Biology 16 (4), pp. 577. Cited by: §1.
 [11] (1983) Compton interaction tomography I. Feasibility studies for applications in earthquake engineering. IEEE Transactions on Nuclear Science 30 (2), pp. 1680–1684. Cited by: §1.
 [12] (201409) Table of integrals, series, and products; 8th ed.. Academic Press, Amsterdam. External Links: Link, Document Cited by: Definition 3.
 [13] (2010) Compton scatter imaging: a tool for historical exploration. Applied Radiation and Isotopes 68 (6), pp. 993–1005. Cited by: §1.
 [14] (2005) On the use of radiation scattering for the detection of landmines. Radiation Physics and Chemistry 73 (1), pp. 7–19. Cited by: §1.
 [15] (2018) Characterization of Comptonscatter imaging with an analytical simulation method. Physics in Medicine & Biology 63 (2), pp. 025016. Cited by: §1.
 [16] (1959) The examination of internal tissues, using gammaray scatter with a possible extension to megavoltage radiography. Physics in Medicine & Biology 4 (2), pp. 159. Cited by: §1.
 [17] (2010) Inversion of a new circulararc Radon transform for Compton scattering tomography. Inverse Problems 26 (6), pp. 065005. Cited by: §1, footnote 1.
 [18] (2006) Imagerie par rayonnement gamma diffusé. Hermès Science. Cited by: §1.
 [19] (1994) Compton scattering tomography. Journal of applied physics 76 (4), pp. 2007–2015. Cited by: §1, footnote 1.
 [20] (2017) Threedimensional imaging of flat natural and cultural heritage objects by a Compton scattering modality. Journal of Electronic Imaging 26 (1), pp. 011026. Cited by: §1.
 [21] (1917) Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten. Akad. Wiss. 69, pp. 262–277. Cited by: §1.
 [22] (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.
 [23] (2018) 3D Compton scattering imaging and contour reconstruction for a class of Radon transforms. Inverse Problems 34 (7), pp. 075004. Cited by: §1.
 [24] (2012) Novel numerical inversions of two circulararc Radon transforms in Compton scattering tomography. Inverse Problems in Science and Engineering 20 (6), pp. 809–839. Cited by: §1.
 [25] (2013) Combined modalities of Compton scattering tomography. IEEE Transactions on Nuclear Science 60 (3), pp. 1570–1577. Cited by: §1.
 [26] (2017) Compton Scattering Tomography: Feature Reconstruction and RotationFree Modality. SIAM Journal on Imaging Sciences 10 (4), pp. 2217–2249. External Links: Document, Link, https://doi.org/10.1137/17M1120105 Cited by: §1.
 [27] (2020) Analytic inversion of a Radon transform on double circular arcs with applications in Compton Scattering Tomography. IEEE Transactions on Computational Imaging 6 (), pp. 958–967. External Links: Document Cited by: §1, footnote 1.
 [28] (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, footnote 1.
 [29] (2020) Compton scattering imaging system with two scanning configurations. Journal of Electronic Imaging 29 (1), pp. 013005. Cited by: §1.
 [30] (2011) Radon transforms on generalized Cormack’s curves and a new Compton scatter tomography modality. Inverse Problems 27 (12), pp. 125001. Cited by: §1, Remark 2, footnote 1.
 [31] (2012) Recent developments on Compton scatter tomography: theory and numerical simulations. In Numerical SimulationFrom Theory to Industry, Cited by: §1.
 [32] (2020) Function reconstruction from reflection symmetric radon data. Symmetry 12 (6). External Links: Link, ISSN 20738994, Document Cited by: Remark 2.
 [33] (2019) Compton scatter tomography in annular domains. Inverse Problems 35 (5), pp. 054005. Cited by: footnote 1.
 [34] (2020) Compton scattering tomography in translational geometries. Inverse Problems 36 (2), pp. 025007. Cited by: §1, §1, §1, Figure 3, §2.1, Definition 1.
 [35] (2019) Microlocal analysis of a spindle transform. AIMS Inverse Problems and Imaging 13 (2), pp. 231–261. Note: External Links: ISSN 19308337, Document, Link Cited by: footnote 1.
 [36] (2018) Three dimensional Compton scattering tomography. Inverse Problems 34 (8), pp. 084001. Cited by: §1, footnote 1.
 [37] (202007) A joint reconstruction and lambda tomography regularization technique for energyresolved xray imaging. 36 (7), pp. 074002. External Links: Document, Link Cited by: §5.2.4.
 [38] (2020) Microlocal analysis of a compton tomography problem. SIAM Journal on Imaging Sciences 13 (2), pp. 746–774. Cited by: §5.2.4.
Comments
There are no comments yet.