The Mathematical Foundations of 3D Compton Scatter Emission Imaging

The mathematical principles of tomographic imaging using detected (unscattered) X- or gamma-rays are based on the two-dimensional Radon transform and many of its variants. In this paper, we show that two new generalizations, called conical Radon transforms, are related to three-dimensional imaging processes based on detected Compton scattered radiation. The first class of conical Radon transform has been introduced recently to support imaging principles of collimated detector systems. The second class is new and is closely related to the Compton camera imaging principles and invertible under special conditions. As they are poised to play a major role in future designs of biomedical imaging systems, we present an account of their most important properties which may be relevant for active researchers in the field.


INTRODUCTION
During the last fifty years, progress in imaging systems using penetrating radiation for biomedical purposes has brought about new topics in mathematics and fueled intense research activities with far reaching results. The mathematics of imaging science has evolved to a full fledged discipline [1]. Transmission computer assisted tomography (CAT-scanning) is based on an integral transform in two dimensions discovered in the sixties by Cormack [2][3][4], who did not realize that J. Radon had already introduced and studied it in his seminal paper [5] in 1917. Subsequently, in an effort to reconstruct directly in three dimensions an object without having to assemble its two-dimensional sections, one is led to consider the so-called X-ray transform or cone beam transform. This transform is an off-spring of the Radon transform and maps a function in R 3 to its (straight) line integrals in R 3 . One way to reconstruct an object is to convert the line data into planar data in R 3 using the Grangeat technique [6]. Then application of the inversion formula of the three-dimensional Radon transform [7] yields the answer. A further generalization of the Radon transform in emission imaging, called the attenuated X-ray transform, accounts for the radiation loss in the traversed medium. The problem of its inversion has been a mathematical challenge for decades and has been solved only in 2001 [8], thanks to complex analysis methods applied to the stationary photon transport problem [9]. The success of this branch of mathematics, coined by Gelfand as integral geometry, goes far beyond the imaging science scope as it has brought significant progress in group representation theory, partial differential equations, boundary value problems, and so forth [7].
The standard Radon transform and its variants are related to a process of data collection in an interaction free propagation of radiation through matter, possibly affected by attenuation. Thus, radiation energy is not altered from emission to detection. However, a sizable part of traveling radiation does suffer an interaction with matter through Compton scattering [10]. Generally as the data quality is lowered, Compton scattering effects have been treated as noise and must be eliminated. So far, in most image processing methods, the aim is to deal away with it, for example, by filtering or by geometric rejection using an antiscatter grid of miniature lead septa. In conventional projection imaging, most of the radiation have been scattered in the patient body, the scatter-to-primary ratio can be as high as 10 [10]. But throwing out all the scattered radiation may not be a smart move since this also means a loss of sensitivity and certainly a loss of valuable information.
This is why we have recently advocated the use of scattered radiation to improve quality in image processing [11] as well as to construct new imaging principles [12]. This proposal generalizes the concept of a Compton camera, proposed by several workers some thirty years ago [13][14][15], in which the idea of electronic collimation is implemented. When radiation is detected at a lower energy than the originally emitted energy, there must be at least one Compton scattering occurring along its propagation path. But the vast majority of scattered radiation is due to mainly single scattering events [16]. To an emerging ray detected at a definite energy there corresponds an ensemble of incoming rays distributed on a circular cone of definite opening angle. The measured data consisting of collected emerging radiation at some detection point is viewed as the integral contribution of a source function on specific circular cones; it may be also called a conical projection of the source. The name of conical Radon transform is attributed to such integral transform and is added to the list of already known Radon transforms on geometric manifolds in R 3 : paraboloids [17], spheres [18], special surfaces [19], second-order surfaces [20], and so forth. Other proposals for use of scattered radiation in X-ray imaging which do not rely on conical projections have been made independently by [21].
Thus, the main topic addressed here is the mathematical framework defining the working principles of a new imaging technique. To keep the discussion transparent, we have adopted an ideal working assumption whereby attenuation is not taken into account. A well-known difficulty is that the nonuniform attenuation along the radiation path leads to tremendous mathematical complication. For a given projection, attenuation correction factors in SPECT rarely exceed 10 in virtually all clinical imaging situations [22]. An exact solution is beyond the scope of the present study. It should be pointed out that, in the case of the attenuated X-ray transform in emission imaging, a comprehensive analytic solution has been attained only recently [9]. Therefore, in practical situations standard corrections for photon attenuation should be envisaged.
In conventional emission imaging, the data collected for object reconstruction is formed by linear projections (cone beam projections) of unscattered radiation filtered by an energy acquisition window which has already discarded from 70% to 80% of the incident photon flux. Moreover, among the retained some have suffered scattering in the collimator and hence must be discarded. It has been reported that collimator scatter increases as the energy of the photopeak of interest increases from a low of 1.9% for Tc-99 m (141 keV) to a high of 29.4% for I-131 (364 keV) with the usual highenergy collimator. The penetration percentage also goes up with energy. Therefore, correction for photons that penetrate through, or scatter in, collimator septa is hardly important at all for Tc-99 m tracers [23].
In the proposed imaging technique, data for reconstruction are provided by the so-called conical projections, which gather radiation from point sources lying on large surfaces inside the emitting object. Although contamination of the scatter component by multiple scatter events may be as high as 30% of the total scatter according to Monte Carlo investigations [16,24], we believe that the signal of single scattered radiation is largely sufficient to make the new imaging principles work, in particular when advanced semiconductorbased detectors with better energy and spatial resolution and sensitivity will become available. These interesting issues will be explored in future work.
In this paper, we present a unified treatment of conical Radon transforms relevant in emission imaging by scattered radiation. In fact, we will be concerned with two classes of conical Radon transforms originating from image formation by Compton scattered radiation on a gamma camera with and without collimator. The first class of conical Radon transform uses circular cone sheets with fixed axis direction (C 1 -cones) whereas the second class deals with circular cone sheets with axis swinging around a point (C 2 -cones). Note that if this point goes to infinity in a given spatial direction, then the C 2 -cone goes over the C 1 -cone. In each case, we will start by showing how the image formation process leads to an integral transform and how this transform is related to the conical Radon transforms. Each conical transform will be introduced and its relevant properties for imaging purposes, in particular their invertibility, discussed. Conclusions and perspectives are given in the last section.

NOTATIONS
Let f be a real nonnegative integrable, smooth, with compact support function on R 3 . The same function in cylindrical (spherical) coordinates is noted f(f).
The definitions of various transforms of f (x, y, z) are as follows: Special functions [25] are as follows: (i) Y (x): Heaviside unit step function; (ii) j l (kr): spherical Bessel function of order l and variable (kr); (iii) P l (cos θ): Legendre polynomial of order l and of variable cos θ; solid angle in the direction of the unit vector k.

WORK SETTING
In this article, we consider the emission imaging problem, that is, the problem of reconstructing in R 3 a gamma-ray radiating object from its Compton scattered radiation data. This object is described by its activity volume density function f . Detection of scattered radiation is performed by a gamma camera in two instances: with or without collimator. The recorded data consists of the coordinates of the detection site, the surface flux density of photons at this site (pixel), and the value of their energy (list mode). Between the radiating T. T. Truong et al. object and the detector stands a scattering medium: it may be a volume or a layer as illustrated by Figure 1. Note that the object itself may also be a scattering medium, and for photon energies above 25 keV, over 50% of the interactions in biological tissues are scatterings [26]. Higher-order scattering events (of much lower probability of occurrence) will be the object in future studies.

Compton scattering
As Compton scattering plays a key role, we will recall some of its properties. The Compton effect discovered in 1923 [27] had served to confirm the particle (photon) nature of radiation, as proposed by A Einstein. Thus, energetic radiation under the form of X-or gamma-rays behave like particles and scatter with electrically charged particles in matter. In biomedical domains, X-or gamma-photons scatter electrons in the biological media they traverse. This scattering process has cylindrical symmetry around the direction of the incoming photon and the energy of the outgoing photon is given by the Compton formula [28] where ω is the scattering angle as measured from the incident photon direction, E 0 , the photon initial energy, ε = E 0 /mc 2 , and mc 2 the rest energy of the electron. Equation (1) shows that single-scattered photons have a continuous energy spectrum in the range Thus an emerging photon with an energy E in a direction of unit vector n may originate from an incoming photon of energy E 0 emitted from a site on a sheet of a circular cone of axis direction n and an opening angle ω.

Compton differential cross-section
At a scattering site M, the number of particles d 2 N sc scattered in a solid angle dΩ sc along a direction making an angle ω with the incident direction follows from the definition of the differential scattering cross-section (see Figure 2) when the following quantities are given: (a) φ in , the incident photon flux density, (b) n e (M)dM, the number of scatterers (electrons) around the scattering site M with volume dM.
We have The Compton differential cross-section has been computed in 1929 by Klein and Nishina [29]. It appears as the product of the area of a disk of radius r e , the classical radius of the electron, and a probability factor P(ω), that is, where r e = 2.82 × 10 −15 m and Thus, we see that, for a given incident energy, the angular distribution of scattered photons is no longer isotropic [28]. Hence, the final form of this number of scattered photons in the direction given by the angle ω is this number is basic to the image formation by scattered radiation.

Conical projections
In standard computer assisted tomography (CAT), the data is gathered under the form of line projections or integrals of a function (attenuation or activity) along straight lines. Here following the scattering mechanism, the data would appear as conical projections or integrals of a function f on a circular cone sheets. The integration is carried out with the Lebesgue measure of the cone in a chosen coordinate system. The result is a function of (i) the coordinates of the cone vertex, (ii) the parameters of the unit vector of the cone axis, (iii) the opening angle of the cone, that is, ω the scattering angle. Figure 3 displays the representations of line projection and conical projection. In the text we will have two types of conical projections corresponding to the two cases of image formation mentioned earlier. The question is now how to use the conical projections to reconstruct the source function f . We will treat the two cases separately in the coming sections.

THE C 1 -CONICAL RADON TRANSFORM
In this section, we consider the possibility of imaging a threedimensional object by collecting data on its scattered radiation on a gamma camera equipped with a collimator and show how the C 1 -conical Radon transform arises. Figure 1 shows the experimental arrangement with the location of the radiating object, the scattering medium, and the collimated gamma camera.

Image formation in gamma imaging by scattered radiation
To concentrate on the scattered imaging principle, we make some simplifying assumptions [12,30]: (1) absence of attenuation for the propagating radiation, (2) constant density of the electrons n e in the scattering medium, (3) isotropic emission from original radioisotopes in object.
To compute the photon flux density at a detection site (pixel) D on the collimated camera we start from (6). Radiation is emitted at point source S, will propagate to scattering site M and reach detection site D.
The incoming photon flux density φ in on scattering site M is now computed from the emission data from point source. Let f (S)dS be the number of gamma photons emitted per unit time by a volume dS in the object around site S. The emission being isotropic, the number of photons emitted in the direction Therefore, the incoming photon flux density at scattering site M is this means that the detected photon flux density at site D is Now, all the contributing point sources S, for given scattering center M, lie on a circular cone sheet of axis identified with − − → MD and opening angle ω, thus we must integrate with the measure δ(cone)dS first. Next, we must take into account all the scattering sites in the scattering medium situated on the line parallel to the collimator axis at site D. Hence, we must perform a second integration with the measure δ(line)dM. To sum up the detected photon flux density at D is 1 SM 2 n e dMπr 2 e P(ω) In the cylindrical coordinate system of Figure 4, the integration measure on the cone is r sin ω dr dφ and the measure T. T. Truong et al.
In practice as f and n e are volume densities, to keep the physical dimensions right, we should think of the cone sheet as having a small thickness e and the line on which M moves as having a tiny section s. These are constants and do not affect the mathematical reasoning, they will be dropped for the sake of expression simplicity, but should always be kept in mind.
Clearly this imaging equation is a compounded integral transform. Assuming that integration interchange is valid, if we define the first integral transform as h x D + r sin ω cos φ, y D + r sin ω sin φ, r cos ω we see that the imaging equation (12) is just a conical projection of the function h on the C 1 -cones with vertex on the plane xOy.
Let us define an interaction factor which also includes e and s: then (12) or in terms of our notations (see Section 2)

Properties of the C 1 -conical transform
The C 1 -conical Radon transform has been discussed in [31], where some of its properties have been studied. However, the problem of inversion will be handled here with a new approach. We will not go through the method of decomposition of functions in circular components but will show that there exists a variant of the central slice theorem [28] which provides the grounds to invert the transform as it is done in the standard Radon transform [28]. First, we observe that by definition in the chosen cylindrical coordinate system the C 1 -conical Radon transform of f , the cones having vertex on the xOy plane, is This can be rewritten under the form of a Fredholm integral equation of the first kind with a delta function kernel concentrated on the sheet of a circular cone [31] with We use now the Fourier representation of the delta function and the two-dimensional Fourier transform of f (x, y, z), f (x, y, z) = dp dqe 2iπ(px+qy) F 2 (p, q, z).
We can perform the integration over z, which restores a three-dimensional Fourier transform F 3 (p, q, r) of f and 6 International Journal of Biomedical Imaging yields a new form of f (x S , y S , ω), Let F 2 (p, q, ω) be Fourier component with respect to the coordinates x S and y S of f (x S , y S , ω), then The last integral of (23) can be computed using polar coordinates in (x, y) and (p, q) spaces, that is, (x−x S ) = ρ cos β, (y − y S ) = ρ sin β, dx dy = ρ dρ dβ, and p = k cos α, p = k sin α, k dk dα. Therefore, Integration on β yields a Bessel function of order zero in the last integral of (24): ρ dρ dβe 2iπ[νρ cos ω+kρ cos (β−α)] = ∞ 0 ρ dρe 2iπρν cos ω 2πJ 0 (2πkρ). (25) Thus, in two-dimensional Fourier space, the C 1 -conical Radon transform appears as Now, we perform the integration on ν which brings back the three-dimensional Fourier transform of f , This step reduces the C 1 -conical Radon transform to a Hankel of order zero: sin ω F 2 (p, q, ω) = ∞ 0 ρ dρ2πJ 0 2πρ p 2 + q 2 F 2 (p, q, ρ cot ω). (28) To perform the inversion of this Hankel transform, we should switch to an appropriate variable, that is, ζ = ρ cot ω. But care must be exercised as far as the range of ω is concerned. We will distinguish two cases.
(a) 0 < ω < π/2, then ζ > 0 as well as tan ω > 0 and (28) can be rewritten as (b) π/2 < ω < π, then ζ < 0 as well as tan ω < 0 and (28) can be rewritten as Application of the Hankel identity [25] 1 leads to the inversion of the Hankel transforms and yields the three-dimensional Fourier components of the object f . As the dual variable to ζ is τ p 2 + q 2 , with τ = tan ω, we get Hence, f can be recovered by inverse three-dimensional Fourier transform.

Remark 2.
In order to reconstruct the object by inverting the compound integral transform given by (12), we need to invert (13). This is quite easy since it can be viewed as a convolution in the variable z D between the two functions 1/z 2 D and f (. . . , . . . , z D + r cos ω). As this operation is not important to the topic of this paper, we refer the reader to [12,30].

Compton camera
That Compton effect which has been proposed as mechanism for imaging is known ever since the fifties. However, there are many ways to do the experimental setups. Most proposals are systems with collimated point source and pointlike detector, see, for example, [21]. In fact the Compton effect is used to probe the electron density of matter and applied often to nondestructive material control. Here we are T. T. Truong et al. interested in data collected by a gamma camera. In the previous section we have examined the case of a gamma camera with a lead collimator which has the disadvantage of rejecting many of the scattered photons. So to improve drastically detection sensitivity, the idea of a Compton camera has been proposed as early as 1974 by many workers [13][14][15]32]. The concept of a Compton camera is analogous to the scheme of Section 4 except that the scattering medium is now a thin scattering layer parallel to the face of a gamma camera without collimator. The data consists of C 2 -conical projections, the cone sheet axis converging to the detection site D; see Figure 5. Note that when D → ∞ in a given direction we recover the C 1 cones.
Following the image formation process in a Compton camera, we will show how a new conical Radon transform, the C 2 -Radon transform, comes up and sketch a proof of its invertibility under specific conditions. The true conical Radon transform of a Compton camera is not yet an analytical inversion formula.

Image formation by scattered radiation in Compton camera
The radiating object stands above the first scattering layer. Its primary rays hit the scattering layer and will be absorbed by the planar camera (see Figure 5). If only photons of energy E below the energy E 0 of the primary photons are recorded, then each detection site collects all possible conical projections coming from all directions in half space, delimited by the photon absorbing detector. This is the principle of electronic collimation which has been designed to improve sensitivity of gamma cameras. To reconstruct an object described by a source function f (x, y, z), we need a set of data consisting of conical projections depending also on three variables. Ideally one could select one detection site D, and consider all the projections along circular cones of opening angle ω and axis swinging around D but with cone vertex constrained to be on a plane. With these conditions a conical projection will depend on three parameters: the scattering angle ω and the two coordinates of the cone vertex on the scattering plane. Thus, we obtain a mapping of f onto a function of three variables. The inverse mapping, when explicitly worked out, would yield a correct imaging procedure by a Compton camera.
Following the assumptions of Section 4.1, the photon flux density at detection site D is evaluated in the same manner as for the case with collimator: × n e dMπr 2 e P(ω)sδ(line) Now in the chosen coordinate system (see Figure 5), it has the expression l being the distance OD, the explicit integration on the cone sheet δ(cone)dS will be given later in Section 5.3 since it is related to the C 2 -conical Radon transform we will be examining.
Up to now there exists only a few attempts to exactly solve this inversion problem. Cree and Bones [33] were the first to consider conical projections on a Compton camera which still has a collimator. Later on Basko et al. [34] as well as Parra [35] have designed inversion techniques based on properties of spherical harmonics as they consider conical projections as made up of cone beam projections, but have not touch really upon the problem of converting cone beam data into conical projection data, a problem similar to the one solved by Grangeat [6] for planar projections in R 3 . There are numerous approximate reconstruction methods, most of them using some back projection techniques (or numerical algorithms) and search for point sources as intersections of cone sheets reconstructed from coincidence measurements on the Compton camera. Lastly let us also cite some other original approaches based on statistical physics [36] as well as algebraic methods [37][38][39].
In this situation, it is of interest to find an analytic inversion of the special C 2 -conical Radon transform responsible for the imaging process in Compton camera. Before tacking this problem, we discuss here a more general problem in which the vertex of the scattering cone is not constraint to lie on a plane (or a surface) as in the Compton camera case. This would give us some freedom to find an inversion formula for another class of C 2 -conical Radon transform. We will come back to the Compton camera in another work. Thus, the mathematics of the Compton camera imaging process have led to a new class of conical Radon transform.

Properties of the C 2 -conical Radon transform
The C 2 -cone is generated by the rotation of a straight line making an angle ω around an axis direction n and meeting which may be regarded as the cone equation in a meridian section within polar coordinates (r, γ) and polar axis n. At ω = π/2 we recover the equation of a plane as a degenerate cone [28].
When evaluating a C 2 -conical projection, we integrate a nonnegative function f on one sheet of the cone. This is equivalent to saying that, giving a "mass" density, we compute the "mass" of a piece of cone surface limited by the intersection curve of the support of f with the cone. This "mass" may be calculated in any convenient coordinate system.
We will use the coordinate system which is expressed in (35). The area element of the cone is the product of the arc element of the line by the element of circle in a plane perpendicular to the cone axis.