Inversion of Weighted Divergent Beam and Cone Transforms

In this paper, we investigate the relations between the Radon and weighted divergent beam and cone transforms. Novel inversion formulas are derived for the latter two. The weighted cone transform arises, for instance, in image reconstruction from the data obtained by Compton cameras, which have promising applications in various fields, including biomedical and homeland security imaging and gamma ray astronomy. The inversion formulas are applicable for a wide variety of detector geometries in any dimension. The results of numerical implementation of some of the formulas in dimensions two and three are also provided.


Introduction
In this paper, we focus mainly on analytic and numerical inversion of an integral transform (cone or Compton transform) that maps a function to its integrals over conical surfaces with a weight equal to some power of the distance from the cone's vertex. It arises in various imaging techniques, most prominently, in modeling of the data provided by the so-called Compton camera, which has novel applications in various fields including medical and industrial imaging, homeland security, and gamma ray astronomy [1,2,5,7,29,35]. In Compton camera setting, the vertices of the cones correspond to the locations of the detection sites on the scattering detector. More information on the working principle of a Compton camera can be found, for example, in [1,5,7,29,32].
It has been mentioned in various papers, e.g. [5,23,31] that, depending upon the engineering of the detector, various power weights can appear in the surface integral. However, more work needs to be done to determine the weight factor that accurately represents the projections obtained from a Compton camera. Several works, e.g. [3-6, 8, 12, 17, 18, 22, 27, 31-33, and references therein] concentrated on the case of pure surface measure on the cone. Here, we consider a weight that is equal to some power of the distance to the vertex (detection site). An alternative inversion formula for such transform that assumes the vertices of the cones are located on a given straight line is provided in [24]. A reconstruction formula for such transform defined on the cones having vertices on a hyperplane and a central axis orthogonal to this hyperplane is derived in [14]. In comparison, the formulas we derive allow for a wide variety of cone vertex (a.k.a. detector or source) geometries, which do not allow for harmonic analysis, but satisfy what we call in this paper Tuy's condition (Definition 3.4).
A closely intertwined with the weighted cone transform is what is called weighted divergent beam transform, which integrates a function over rays with a weight equal to some power of the distance to the starting point (source) of the ray. We thus study it in some details, which leads eventually to the desired weighted cone transform inversions. When the weight factor is not present, this is the well studied and important for the 3D X-ray CT divergent (or cone) beam transform (see e.g. [11,13,15,19,20,25,30,34, and references therein]).
In order to avoid being distracted from the main purpose of this text, we assume that the functions in question belong to the Schwartz space S of smooth fast decaying functions. This allows us to skip discussions of applicability of various transforms. However, as in the case of Radon transform (see, e.g. [21,26]), the results have a much wider area of applicability, since the derived formulas can be extended by continuity (although we do not do this in the current text) to some wider functional spaces. This is confirmed, in particular, by our successful numerical implementations for discontinuous (piecewise continuous) phantoms. The issues of appropriate functional spaces will be addressed elsewhere.
We also adopt the standard abuse of notations, writing the action of a distribution T on a test function ϕ, T, ϕ , as T (x)ϕ(x)dx.
The paper is organized as follows. In section 2, we define the weighted divergent beam and cone transforms, and describe a simple relation between them. In section 3, we present a variety of inversion formulas for the weighted divergent beam transform (Theorems 3.6 and 3.8). We then derive another integral relation between the weighted divergent beam and cone transforms, which leads to new inversion formulas for the n-dimensional weighted cone transform (Theorem 3.15). In section 4, we investigate the relation between the Radon and weighted divergent beam and cone transforms. This enables us to derive other inversion formulas for the latter two (Theorem 4.4). Section 5 contains the results of numerical implementation of some of the inversion formulas for the weighted cone transform in dimensions two and three for two different vertex geometries, as well as examples of numerical inversion of two weighted divergent beam transforms in dimension three. Conclusions and remarks can be found in section 6, followed by the acknowledgments section.

Weighted Divergent Beam and Cone Transforms
In this section, we define the closely related weighted divergent beam and cone transforms.
where u ∈ R n is the source of the beam {u + ρσ}| ρ≥0 and σ ∈ S n−1 is the unit vector in the direction of the beam.
Consider now a circular cone 1 S in R n . Its surface can be parametrized by a triple (u, β, ψ), where u ∈ R n is the cone's vertex (apex) 2 , the unit vector β ∈ S n−1 is directed toward cone's interior along the cone's axis, and ψ ∈ (0, π) is the opening angle (see Fig. 1
where dS is the surface measure on the cone S. In other words, where dx is the Lebesgue measure on R n . Remark 2.3. We note that k = n − 2 corresponds to the case of pure surface measure on the cone. Remark 2.4. We will say just "weighted" cone or divergent beam transform, when no confusion about the value of k can arise.
Changing variables in (3) as x = u + ρσ for ρ ∈ [0, ∞) and σ ∈ S n−1 , and using the fact that δ is homogeneous of degree −1, we make the following simple observation: 1 The word "cone" in this paper always means a surface, rather than solid cone. 2 In the Compton camera imaging, cone's vertex corresponds to a detection location. 3 At this step, one can allow all real values k > −1, while later on in sections 3 and 4, k ∈ Z + := {0, 1, 2, . . . } will be important.
By letting t = cos ψ, we can rewrite (with an abuse of notation) C k f as

Inversion of Weighted Divergent Beam and Cone Transforms
In this section, we present a variety of inversion formulas for the weighted divergent beam transform. We then derive another integral relation between the weighted divergent beam and cone transforms, which enables us to develop new inversion formulas for the n-dimensional weighted cone transform.

Inversion of the Weighted Divergent Beam
Transform. If f ∈ S(R n ), for each u ∈ R n , D k u f (σ) can be uniquely extended to a smooth function on R n \{0} homogeneous of degree −(k + 1): This function is locally integrable with respect to x ∈ R n , provided k < n − 1, and has a well-defined Fourier transform as a tempered distribution (see e.g. [10]), i.e., for each ϕ ∈ S(R n ), In the following, we derive inversion formulas for the divergent beam transform that are analogs of the well known [25,34] Tuy's inversion formula, which addresses the case when k = 0 in dimension three, and the sources (detectors in the Compton camera case) move along a curve. Definition 3.1. In the rest of the paper, the shorthand notations ∂ uj and ∂ u will be used for the partial derivatives ∂/∂u j and gradient ∇ u with respect to the variables u.
where ∆ u := j ∂ 2 uj is the Laplace operator with respect to the variable u, and its power when k is not an odd integer is understood as the corresponding Riesz potential (see, e.g. [26]).
Proof. Let f ∈ S(R n ). For any θ ∈ S n−1 , the Fourier transform of D k u f satisfies Indeed, for any ϕ ∈ S(R n ), due to D k u f being homogeneous of degree −(k + 1) and changing to polar variables y = sω, we have Now, changing variables in s to ρ = s/r, and then letting y = u + rω, we get which implies (9). The following simple formula holds for any unit vector θ: Thus, applying (k+1)/2-th power of the Laplace operator with respect to u to (9), we obtain Now, recalling the Fourier inversion formula in polar coordinates and comparing with (11), we obtain the desired formula Remark 3.3. Considering formula (8), one realizes quickly that it is not very useful, since it requires "sources" u of the beams to be available throughout the whole space. In the Compton camera case, as well as in 3D CT, this would require detectors/sources to be placed throughout the object imaged, which is impossible. Moreover, in this case, one deals with just a deconvolution problem, and a severely overdetermined one at that (the dimension of the data used is 2n − 1 instead of n). Thus, there must exist formulas requiring much less data, in particular allowing the detectors u to be situated only outside the object being imaged (e.g. Tuy's formula only requires an arc of external sources). This is also related to the interesting question about "admissible" complexes of cones that provide enough data for stable reconstruction. We have already briefly addressed this issue in [22,32] and plan to have more detailed discussion elsewhere.
Here we show an example of how such deficiency can be alleviated for the weighted divergent beam transform.

Definition 3.4.
• Let M ⊂ R n be a smooth d-dimensional submanifold. We will say that it satisfies the Tuy's condition with respect to a subset V ⊂ R n , if any hyperplane intersecting V has a non-tangential intersection with M.
Equivalently: for any x ∈ V and unit vector θ ∈ S n−1 , there exists a point u ∈ M such that θ · x = θ · u, and θ is not normal to M at the point u.
• We denote by P u the orthogonal projection onto the tangent space to M at the point u ∈ M.
Remark 3.5. Notice that the above condition is a strengthened version of what was called admissibility condition in [22,32].
Theorem 3.6. Let k be an odd natural number and M ⊂ R n satisfies Tuy's condition with respect to a compact V . Then, for any homogeneous linear elliptic differential operator L(u, ∂ u ) of order k + 1 on M and any smooth function f supported in V , the following inversion formula holds: where u ∈ M is related to x and θ as in the Tuy's condition.
If M is one-dimensional, then k can be any natural number (in this case, when k = 0, one ends up with the standard Tuy's formula).
Remark 3.7. Notice that u ∈ M in (13) depends on both x and θ and that L(u, P u θ) does not vanish if the Tuy's condition is satisfied.
Proof. The proof follows exactly the one of Theorem 3.2, using at the end the formula for the symbol of a homogeneous differential operator of order k +1 (instead of a power of the Laplacian in (10)): (14) L(u, ∂ u )e iρθ·u = (iρ) k+1 L(u, P u θ)e iρθ·u and noticing that the factor L(u, P u θ) does not vanish, due to the ellipticity and homogeneity of the operator and the Tuy's condition.
A serious deficiency in Theorem 3.6 is that, unless M is one-dimensional, only odd values of k are allowed. This issue can be resolved, paying the price of having a more complex formula.
Indeed, consider the following first order linear differential operator acting tangentially to M, with coefficients depending upon u ∈ M and θ ∈ S n−1 : Let also L(u, ∂ u ) be an operator like in Theorem 3.6, but of order k. Applying the composition O • (aL) to the exponential e iρθ·u , where a := 1/L(u, P u θ), and using (14), we get Since the order of the composition is odd, this enables us to extended the inversion formula to the even values of k: Theorem 3.8. Let k be an even natural number and M ⊂ R n satisfies Tuy's condition with respect to a compact V . Then, for any homogeneous elliptic differential operator L(u, ∂ u ) of order k on M and any smooth function f supported in V , the following inversion formula holds: where u is related to x and θ as in the Tuy's condition.
Proof. The proof stays exactly the same, except instead of (14), we use (16).
Remark 3.9. Since we are dealing with a severely overdetermined transform, the variety of possible inversion formulas is large and is not exhausted by the ones above. For instance, instead of using the operator 1

3.2.
Inversion of the Weighted Cone Transform. The following result presents a relation between the weighted divergent beam and cone transforms, which will be instrumental in the inversion of the latter one.
(Notice that C k f (u, β, t) = 0 for |t| > 1, and is smooth for |t| < 1.) Proof. By the representation (5) of the weighted cone transform, we have We note that when h is a regular distribution near t = ±1, the identity (18) reads as Definition 3.11. Let k < n − 1. We introduce the following distribution: Proposition 3.12. The following identity holds: Proof. For each ϕ ∈ S(R n ), since D k u f is homogeneous of degree −(k + 1), we have which implies (21). Now, by combining (19) and (21), and using the inversion formula for the weighted divergent beam transform (8), we obtain an inversion formula for the weighted cone transform.
Theorem 3.13. Let f ∈ S(R n ), and k < n − 1. Then, Remark 3.14. The same deficiency applies here that was mentioned in remark 3.3: the formula requires the cones to be available with all vertices u throughout the space, which is unacceptable for many imaging applications (e.g. Compton ones). Fortunately, a similar remedy as for the divergent beam transform exists, which we address next.
Theorem 3.15. Let k ∈ Z + . Suppose that M ⊂ R n satisfies the Tuy's condition with respect to a compact V ⊂ R n , and f is a smooth function supported in V . Then, depending on the parity of k, the following inversion formulas hold: (i) if k is odd, then for any homogeneous linear elliptic differential operator L(u, ∂ u ) of order k + 1 on M (ii) if k is even, then for any homogeneous linear elliptic differential operator where O is given as in (15), and u ∈ M is related to x and θ as in the Tuy's condition.

Relations with the Radon Transform: Other Inversion Formulas
In this section, we present a relation between the Radon and weighted divergent beam and cone transforms, and from it we develop other analytical inversion formulas for the latter two in any dimension, when k ∈ Z + . We recall first that the n-dimensional Radon transform R maps a function f on R n into the set of its integrals over the hyperplanes of R n . Namely, if ω ∈ S n−1 and s ∈ R, (25) Rf (ω, s) = R ω f (s) := x·ω=s f (x)dx.
In this notation, the Radon transform of f is the integral of f over the hyperplane perpendicular to ω at the signed distance s from the origin. We will use the following well known inversion formula for the Radon transform (see e.g. [9,16,26]). For any f ∈ S(R n ), where H is the Hilbert transform in R defined as the principal value integral   (ω, s). We now present a relation between the Radon and the weighted divergent beam and cone transforms for any dimension n and any k ∈ Z + . Analogous relation for the usual divergent beam transform (k = 0) is obtained in [15] (see also [25,Chapter 2]).
Let k ∈ Z + , and h be the function on R defined by if k > n − 2 and k − n is odd, if k > n − 2 and k − n is even, We note that h is homogeneous of degree k − n + 1, and for k > n − 2, h (k−n+2) = δ(t) (see e.g. [10]).
Proof. The first equality is already obtained in Proposition 3.10. By definition of the weighted divergent beam transform, we have due to the homogeneity of h. Letting x = u + ρσ, we obtain Remark 4.2. In dimension three, for k = 1, the relation (29) gives the following (geometrically obvious) formula: which is used and in [5].

Remark 4.3.
We notice that Proposition 4.1 is valid for any choice of h as long as it is homogeneous of degree k − n + 1 and regular around ±1. Indeed, h(t) = t k−n+1 would work, too. In fact, the relation (29) is proven for such h and is used to derive an inversion formula for the cone transform for k = 0 and k = 1 in dimension three in [31], and for k = 0 in dimension two in [1,17]. For the usual divergent beam transform (k = 0) in dimension three, the applications of various functions h can be found in [25, Chapter 2, and references therein]. Now, using the differentiation property of the convolution and the inversion formula (26) for the Radon transform, we obtain the following formula, which can be used for inversion of both the weighted divergent beam and cone transforms: Suppose that for any u ∈ M and β ∈ S n−1 , s = u · β and Then, for any f ∈ S(R n ), where G (k+1) is the (k + 1)-st derivative of G with respect to s, h is given in (28), and H is the Hilbert transform (27) with respect to s.
Proof. Using h given in (28) for k ≥ n − 2, we have and for k < n − 2, we get Hence, the application of Radon transform inversion (26) gives the result.
Remark 4.5. The reconstruction formula (31) is independent of the geometry of the manifold M that u belongs to. Indeed, for the weighted cone transform, it is sufficient that for any s ∈ R and β ∈ S n−1 , there is a vertex u = sβ + y, where y⊥β, and the weighted cone data is available at all angles ψ for these u and β. In other words, the requirement for the reconstruction is that any hyperplane intersecting the domain of reconstruction contains the vertex of a cone with the axis normal to the plane and all opening angles. For the weighted divergent beam transform, the corresponding condition is that for any s ∈ R and β ∈ S n−1 , there is a source u = sβ + y, where y⊥β, and the weighted divergent beam data is available at all directions σ for this source u.

Numerical Implementation of Theorem 4.4
In this section, we present the results of numerical implementation of Theorem 4.4. In dimension two, the weighted cone transform and divergent beam transforms are similar, so we only provide examples for the weighted cone transform for k = 1 using two different vertex geometries. We then give the reconstruction results for the weighted cone transform in dimension three for k = 0 and k = 2 using a spherical vertex geometry. Examples showing the reaction of the algorithms to Gaussian white noise in the data are also provided. We also present an example of numerical inversion of the weighted divergent beam transform in dimension three for the cases k = 1 and k = 2 using a spherical source geometry.
All phantoms are placed off-center of the vertex curve/surface, to avoid unintended use of rotational invariance. Care was taken to avoid other possibilities of committing an inverse crime, by making the forward and inverse algorithms as unrelated as possible.
5.1. 2D Image Reconstruction from Cone Data for k=1. In dimension two, for k = 0, the relation (29) gives C 0 (u, β, π/2) = Rf (β, u·β), which is geometrically obvious (also see [5]). Thus, we focus here on the case k = 1 only. Since a cone in 2D is represented by two rays with a common vertex, the weighted cone transform for k = 1 is given by where β(φ) = (cos(φ), sin(φ)). The inversion formula (31) now reads as For the numerical implementation of (32), we considered the phantom where D 1 and D 2 are the concentric disks centered at (0, 0.4) with radii 0.25 and 0.5, respectively. Here, χ Di denotes the characteristic function of each disk (see Fig. 2). The cone projection data of the phantom f is simulated numerically using 256 counts for vertices u, 400 counts for central axis directions β and 90 counts for opening angles ψ.   As our inversion formula is valid for arbitrary geometry of vertices, we considered both a circular and a square geometry of vertices of the cones. Figure 3 shows the results of reconstruction from cone projections where the vertices cover the unit circle. The density plot and the y-axis profile of the reconstruction are provided in (a) and (b), respectively. The results with a 5% Gaussian white noise added to the cone data is shown in Figure 4.
In Figure 5, we provide the results of reconstruction from cone projections where the vertices are placed along the sides of a square with sides of length two. The density plot and the y-axis profile of the reconstruction are provided in Figure 5 (a) and (b), respectively. Figure 6 shows the results with a 5% Gaussian white noise added to the cone data.  In the case of the square geometry (but not in the circular one), some cornerrelated effects appear along the diagonals, as shown in Figure 7. They can be eliminated by using a finer discretization in β. 5.2. 3D Image Reconstruction from Weighted Cone Data for k=0 and k=2. In dimension three, the relation (29) is geometrically obvious for the case k=1 (see remark 4.2), and is numerically implemented in [5] using spherical harmonic expansions. Here, we provide examples of reconstruction from the weighted cone data for k = 0 and k = 2. Theorem 4.4 gives the following inversion formula for k = 0: and for k = 2, we have In our examples, the vertices of the cones cover the unit sphere S 2 in R 3 and the phantom is the characteristic function of the 3D ball of radius 0.5 units located strictly inside and off-center of this sphere.
The forward simulations of weighted cone projections were done numerically using 1800 counts for vertices u on the unit sphere, 1800 counts for unit vectors for the cone axis directions β and 200 counts for opening angles ψ. For the discretization of the sphere, we used a uniform mesh for both the azimuthal and the polar angles. Figure 8 shows the three cross sections of the spherical phantom and of its reconstructions from the cone data obtained via (33). The comparison of the phantom and the reconstruction given in Figure 8 in terms of their coordinate axis profiles is provided in Figure 9.   Figure 10 shows the three cross sections of the spherical phantom and of its reconstructions from the cone data obtained via (34). The comparison of the phantom and the reconstruction given in Figure 10 in terms of their coordinate axis profiles is provided in Figure 11.  Figure 11. Comparison of the x-axis (left), y-axis (center) and z-axis (right) profiles of the reconstruction and the phantom given in Fig. 10. Figures 12 and 13 show the results of reconstruction from weighted cone data for k = 2 contaminated with 5% Gaussian white noise.  5.3. 3D Image Reconstruction from Weighted Divergent Beam Data for k=1 and k=2. In dimension three, when k = 0, Theorem 4.4 reduces to Grangeat's formula [13]. Here, we provide examples of reconstruction from the weighted divergent beam data for k = 1 and k = 2 using a spherical source geometry. For k = 1, Theorem 4.4 gives the following inversion formula: and for k = 2, we have The forward simulations of weighted divergent beam projections were done numerically using 1800 counts for sources u on the unit sphere, 30K counts for unit directions σ. For the triangulation of the sphere, we used the algorithm given in [28]. Figure 14 shows the three cross sections of the spherical phantom and of its reconstructions from the weighted divergent beam data obtained via (35). The comparison of the phantom and the reconstruction given in Figure 14 in terms of their coordinate axis profiles is provided in Figure 15.   Figure 16 shows the three cross sections of the spherical phantom and of its reconstructions from the weighted divergent beam data obtained via (36). The comparison of the phantom and the reconstruction given in Figure 16 in terms of their coordinate axis profiles is provided in Figure 17.

Conclusions and Remarks
• In this paper, we present several novel inversion formulas for the weighted divergent beam and cone transforms. In most combinations of the weights and dimensions, such formulas have not been known before. Even in the non-weighted case, the formulas are different from developed previously, in particular in [22,32]. The main (but not the only one) trigger for such studies is applicability to Compton camera imaging. • One of the most important features, in the authors' view, is that the new formulas are adjustable to a wide variety of (detector) geometries. We introduce the class of such geometries satisfying what we call in the text the Tuy's condition (its weaker form was called admissibility in [22,32]). Most of the previously derived formulas required very symmetric geometries, allowing for harmonic analysis tools to be used. • The latter remark is related to the important issue of understanding the geometries that allow for (stable) reconstruction. They deserve a much more thorough study, which we plan to address in another publication.
• As it was mentioned in the introduction, to avoid being distracted from the main purpose of this text, we assume that the functions to be reconstructed belong to the Schwartz space S. However, as in the case of Radon transform (see, e.g. [21,26]), the results undoubtedly have a much wider area of applicability, since the derived formulas can be extended by continuity to some wider function spaces. Although we do not do this in the current text, this conclusion is confirmed, in particular, by our successful numerical implementations for discontinuous (piecewise continuous) phantoms. The issues of appropriate function spaces will be addressed elsewhere. • Practical soundness of the derived inversion techniques is shown by their numerical implementation in the most interesting dimensions two and three. One should also notice, that the new algorithm of the 3D cone transform inversion works much faster than some of the ones developed in [22]. The reason is that a much coarser mesh (1.8K nodes) on the sphere suffices, rather than 30K used in [22].