Exact cone beam reconstruction formulae for functions and their gradients for spherical and flat detectors

We derive unified inversion formulae for the cone beam transform similar to the Radon transform. Reinterpreting Grangeat’s formula we find a relation between the Radon transform of the gradient of the searched-for function and a quantity computable from cone beam data. This gives a uniqueness result for the cone beam transform of compactly supported functions under much weaker assumptions than the Tuy–Kirillov condition. Furthermore this relation leads to an exact formula for the direct calculation of derivatives of the density distribution; but here, similar to the classical Radon transform, complete Radon data are needed, hence the Tuy–Kirillov condition has to be imposed. Numerical experiments reported in Hahn B N et al (2013 Meas. Sci. Technol. 24 125601) indicate that these calculations are less corrupted by beam-hardening noise. Finally, we present flat detector versions for these results, which are mathematically less attractive but important for applications.


Introduction
In the following we consider the x-ray reconstruction problem in three dimensions when the data is measured by firing an x-ray tube emitting rays to a 2D detector. The movement of the combination source-detector determines the different scanning geometries. In many industrial applications the source is moved on a circle around the object. Then the Tuy-Kirillov condition, demanding complete data, is not fulfilled. For deriving a Grangeat-type relation this is not necessary. Also a uniqueness result can be formulated under much weaker conditions. For the closed form inversion formula, this condition is essential, similarly to the Radon inversion formulae. We present results for reconstructing both the function itself as well as the direct calculation for the gradient, which seems to be less affected by beam hardening as results from real data presented in [9] show.
Cone beam reconstruction formulae have a long history. The famous Feldkamp algorithm [4] is an approximation assuming that the source is sufficiently far away from the object, so that data in almost parallel planes are sampled. The reconstruction itself is then achieved by special weighting of two-dimensional algorithms. For a long time and still nowadays this algorithm is applied in nondestructive testing. For exact reconstructions the paper of Tuy, [23] plays an important role, not for the implementation, but for the condition stating that the cone beam data can be considered as complete three-dimensional Radon data. A relation between cone beam data and three-dimensional Radon data has been explicitly established by Grangeat [7]. His algorithm was numerically unstable. Finch [5] gave an inversion algorithm, which was very attractive, but inefficient. In the authorʼs group approximate inversion formulas were developed [3,17] using the approximate inverse as an inversion tool. Katsevich [11,12] developed an exact inversion algorithm based on Grangeatʼs formula and its derivative. Zhao et al [24] presented a general framework for deriving inversion formulas. Gel'fand and Graev [6] started a completely new direction for inversion formulae using the Hilbert transform and rebinning of the data. From the many papers in this direction we mention [2]. Due to the immense importance of the field and the many contributions the list cannot be complete.
In the present paper we derive new inversion formulae again based on the formula of Grangeat and inversion formulae for the Radon transform. The result, where first the data are differentiated, was presented in [15]. Using generalized Funk or Funk-Radon transforms this can now be presented in condensed form. Using the hemispherical transform, [21] we find a formula of rho-filtered layergram type, where first a backprojection and then a differentiation is performed. Motivated by feature reconstruction, [16] we derive formulas for the direct reconstruction of derivatives of the searched-for function. This is more efficient and more stable than first reconstructing the density and then applying formulas for numerical differentiation. From Grangeatʼs formula we easily find a relation between the gradient of a function and its cone beam transform.
Finally we derive those results also for flat detectors. In principle this is a mere change of coordinates from spherical to flat detectors. But the direct derivation of the analogue of Grangeatʼs formula allows for simpler representation of the results.
The paper is organized as follows. Section 2 summarizes some results for x-ray and cone beam transform and their relation. In section 3 the other integral transforms used in the paper are discussed, namely generalized Funk transforms, backprojections and the Radon transform. In section 4 we present a general approach for deriving relations between different integral transforms, resulting in a generalization of the Grangeat formula, relating the Radon transform of the gradient of a function with the Funk transform of its cone beam transform. This leads to a new uniqueness theorem for the cone beam transform. Finally we describe the change of coordinates from Radon to cone beam data. Section 5 contains the inversion formulae for the cone beam transform and section 6 its versions for the reconstruction of derivatives. Finally in section 7 we present the corresponding results for flat detectors, which are important for practical applications.

X-ray and cone beam transform
The x-ray transform of a function is defined via line integrals of the function. Classically we orient the lines in a parallel fashion leading to . Consequently the data space has four dimensions where the searched-for function is a three-dimensional object, the reconstruction problem is overdetermined.
In synchrotron applications the object is rotated around one axis, leading to q q Î´ŷ S , 0 2 ( ) , where S 0 2 is a great circle on the sphere S 2 . Hence the collected data are three-dimensional. The reconstruction problem is reduced to slice-by-slice reconstruction or more efficiently special inversion formulas can be designed in order to combine in one step all data for one direction, see [8].
In cone beam tomography the rays are differently parametrized. Here we start from a source position a and consider all rays emanating from there, leading to The source position is moved on a curve Γ outside of the support of f, hence the data set is given on = GT S. 3 2

( )
The dual operator of D as mapping between L 2 spaces is given as ⎛ This transform can be related to the above x-ray transform. Due to the fact that in the case of the x-ray transform integration is performed over the whole line, here equality only holds if the direction θ is pointing from the source position towards the object. Denoting with q E x the orthogonal projection of x to q^we observe where H denotes the Heaviside function.

Related integral transforms
In this section we summarize some transforms used in the following. Generalized Funk transforms play an important role for the filtering of the data on the detector as well as in the step after changing the coordinates from cone beam to Radon. Different weights appear in the backprojection leading to the generalized transforms. Finally, the Radon transform and inversion formulae are presented.

Generalized Funk transforms
With the Delta distribution δ and its derivatives we define for   -Î k 1 the following transforms. See also [10].
For k=0 we have the classical Funk or Funk-Radon transform Finally we note that with = - where H is the Heaviside function we get the hemispherical transform, see Rubin [21],

Generalized backprojection for the cone beam transform
In the following we introduce different weights as used for the different reconstruction formulae and the reconstruction of derivatives. ⎛ where * D is the dual transform of D as mapping from W L 2 ( ) to L T 2 ( ), see (4).

Radon transform
The Radon transform of a function is defined as integrals over hyperplanes. Hence This transform maps from a three-dimensional space to a three-dimensional space, hence the problem is in contrast to the x-ray inversion not overdetermined. Finally we mention inversion formulae, see [19] Here a -I is the Riesz potential defined via Fourier transform as

Formulae of Grangeat-type, uniqueness and the Tuy-Kirillov condition
For deriving inversion formulas for a transform it is helpful to relate this transform to one where inversion formulae are known. One possibility to derive a relation between different transformations T 1 and T 2 is to find functions or distributions h 1 and h 2 such that the application of the dual transforms results in The most prominent example in tomography are Fourier,  = T 1 and Radon transform for leads to the projection theorem. There is a simple relation between the 3D-Radon transform and two 2D x-ray transform relying on Fubini's theorem for expressing integrals over a plane, first used by Maass [18] to derive a singular value decomposition for the parallel beam x-ray transform. The relation between Radon and cone beam transform is not as simple. Essentially one has, going from the cone beam transform to Radon transform, polar coordinates which results in the measure w t t d d . In order to get the additional factor t in the measure we differentiate the cone beam transform with respect to θ.
Formally we use the above mentioned strategy with T 1 the Radon transform for fixed direction ω and T 2 the cone beam transform for fixed source position a. Consequently h 1 has to be a function of s where h 2 is a function of the direction of the x-ray θ.
and with Assuming h to be homogenous of degree −2 we see Choosing h to be the first derivative of the Delta distribution, known to be homogenous of degree −2, we find With the generalized Funk transform this takes on the form of the formula of Greangeat: where  1 ( ) acts on the second variable of f D . A further differentiation of this formula leads together with (13) to the inversion formula of Katsevich, [11]. Using ) we can get the following cone beam to Radon map for the gradient of the searched-for function.
For sufficiently smooth and compactly supported functions this relation directly leads to uniqueness results for cone beam tomography. We assume in the following theorem that the support of f is contained in a ball B ρ around 0 with radius ρ and that the searched-for function belongs to the Sobolev space  ) for all ω ä S 1 and all s ä I 1 .
Let the scanning curve G fulfill the following condition: The geometric interpretation is that all planes intersecting f supp( ) and in addition are orthogonal to all w Î S 1 , intersect Γ.
) . Now the uniqueness result for the Radon transform results in a uniquely determined  Î r f L B 2 3 ( ( )) . The homogeneous boundary conditions for the compactly supported f leads to the unique solution for the cone beam problem.
With the uniqueness results for L 2 functions from [22], and generalizations to other integration measures in [1], we can state the following uniqueness result for the cone beam transform. For uniqueness theorems for the Radon transform it is obviously not necessary to have complete data. As example in the two-dimensional case the limited angle problem may serve.
where ¢ a is the derivative of the scanning curve in a. The first part says that complete Radon data are available, and the transversality condition allows for transforming the integration variable from Radon data to those of cone beam data. In contrast to the above given condition (19), all planes intersecting the support of f have to intersect the scanning curve, not only those planes perpendicular to w Î S 1 . When  w = s a then we have to count how may source positions Î G a fulfill this condition. This number is called the Crofton symbol n :

Inversion formulae for the 3d cone beam transform
Starting from the formula of the Grangeat (17) relating cone beam and the Radon transform and using the two inversion formulae for the Radon transform (13,14) we can formulate the following result. Note that (27) was, not in this condensed form, first given in [15]. The starting point for the second formula is now again (17) and the second of the inversion formulae for the 3D Radon transform. Instead of the mostly used form, where first a derivative of the data is calculated and then the backprojection is performed we chose the version where the Laplacian is applied to the result of the backprojection. In two dimensions this is not appropriate because there the operation after the backprojection is neither local nor easily computable. Here we profit from the localness of the inversion formula of the Radon transform in odd dimensions. We observe that the number of differentiations in both formulas is the same, when we interpret  -1 ( ) as an inverse differentiation. In (27) the differentiation operator acts on the direction, which together with the discontinuity of the Crofton symbol needs differentiation in the distributional sense. In contrast to this, in (28) we integrate over the directions, and later differentiate with respect to the spatial variable.
As consequence of theroem 5.1 we note the following.
we calculate the second derivative of f R by applying again the derivative of the Delta distribution in order to use the first inversion formula for the Radon transform. Only then we change the coordinates to those of the cone beam transform avoiding that way the tedious differentiation similar to the formula of Katsevich, [11], as first performed in [15].
. With the assumption that the Tuy-Kirillov condition is fulfilled we can change the variables as  w = s a with the transformation of the integration measure as described in the last section (23) resulting in and perform the backprojection of the Radon transform to get Using again the homogeneity of d¢ of degree −2 we get after changing the order of integration With the above used change of variables we get where we used that H is homogeneous of degree 0, so that this integral over S 2 can be rewritten as   . Finally we apply the Laplacian.

Calculation of derivatives
For finding inversion formulae for calculating directly the derivatives of the searched-for function from the cone beam data we start from theorem 4.1.
Theorem 6.1. The gradient of the searched-for function is determined using the vector valued multiplication with the multiplication operator M Γ given in (25) as Proof. We apply the two inversion formulae (13,14) on the result of theorem 4.1 to get in the first case with and as in the proof of theorem 5.1 as

Applying the backprojection with
where in the final step we have used that d is homogeneous of degree −3. Similarly we get starting from (14) with where we used that δ is homogeneous of degree −1.

Flat detector versions
Mathematically equivalent are the results for flat detectors. Due to their importance for practical applications we derive them in the following.

Cone beam transform for flat detectors and its dual
As mathematical model we consider line integrals starting from a source position a and hitting a virtual detector through 0 perpendicular to a at the position h Îâ . In contrast to the classical notion we parametrize the rays now as leading to the transform, which we now call X instead of D to make the differences obvious.
Hence as the second variable we use instead of a direction q Î S 2 for a flat detector the points on the detector h Îâ .

Inverse Problems 32 (2016) 115005 A K Louis
To determine the dual operator we calculate the intersection of a ray emanating from a through x with the detector. The detector plane is spanned by  a a a a a a = Î = a S span , , , with 0. 35 : , , , 39 Now we are in the position to calculate the dual operator for fixed source position a.
which completes the proof.

Inverse Problems 32 (2016) 115005 A K Louis
We now consider also the dual operator of the cone beam transform for a source curve  G Ì 3 which is outside of the support of the searched-for function f. :

Flat detector version of the formula of Grangeat
To find a relation between the Radon transform and the flat detector cone beam transform we apply the strategy described in the beginning of section 4.
The argument of the function h above can be calculated as The main difference between the spherical and flat detector inversion is in the data filtering step. For spherical detectors we apply the Funk transform of order 1 on the sphere, for flat detectors it is a derivative of a line integral over the detector. The multiplication operators differ only by a factor   a . Formally the dual of the cone beam transform * D m appears, but it stems from the backprojection of the Radon transform.
For numerical implementation either the idea from the approximate inverse, [14,16], or direct regularization of the differential operators, [13] have to be applied. Faster implementations are more likely possible with the concept of approximate inverse, as future research may show.

Conclusion
For these closed form inversion formulae it is unavoidable that the Crofton symbol and especially the derivative of the scanning curve appear. For practical measurements at discrete source positions a the way from one to the next source position plays no role, the data are simply independent of the way. So, a further challenge is to derive inversion methods for those semi-discrete problems.