Fast Laplace Transforms for the Exponential Radon Transform

The Fourier slice theorem for the standard Radon transform generalizes to a Laplace counterpart when considering the exponential Radon transform. We show how to use this fact in combination with algorithms for the unequally spaced fast Laplace transform to construct fast and accurate methods for computing both the forward exponential Radon transform and the corresponding back-projection operator.


Introduction
One of the most popular ways of computing the standard two-dimensional Radon transform for parallel beam geometry relies on Fourier transforms via the Fourier slice theorem. A standard discretization with equally spaced samples in the spatial variables yields a relation between the sampled spatial data and its Fourier transform on a polar lattice. This can be rapidly evaluated by using a combination of standard FFT and algorithms for unequally spaced FFT. In this way the time complexity is essentially reduced by one order, from O(N 3 ) to O(N 2 log N ). The main purpose of this paper is to show that a similar reduction in computational cost can be obtained for the exponential Radon transform, by relying on algorithms for unequally spaced fast Laplace transforms instead of the counterpart for Fourier transforms. The same holds for the adjoint operator, sometimes called the exponential back-projection operator, which in particular enables us to implement a fast algorithm for the inverse exponential Radon transform, based on a formula by Tretiak and Metz [21]. The need for fast transforms is particularly important when employing iterative schemes for the solution of (ill-posed) inverse problems, arising e.g. from incomplete angular coverage in the measurements, or for employing noise reduction techniques.
Upon discretization, the Radon transform becomes an operator on finite dimensional space, having a matrix representation. The same goes for the inverse and adjoint Radon transform, and these discretizations are associated with approximation errors which depend on the choice of specific parameters etc. On top of that, there are errors that arise when replacing the discrete operator with its fast counterpart, which is the topic of this paper. We show that the errors that come from the latter part can be controlled to arbitrary precision level at a computational cost that is in practice proportional to − log( ). In this way, reconstruction issues concerning for instance noise sensitivity can be isolated to the (standard) sampling setup, with only minor influence coming from the proposed fast computational algorithms.
The inverse problem of reconstructing data modeled by the exponential Radon transform is ill-posed, meaning that additional regularization techniques often need to be applied in order to obtain suitable reconstructions of data. There are several different regularization approaches in the literature, and the particular type of regularization will typically depend on the structure of the sample measured on. We therefor restrict our attention to how to construct fast algorithms for computing the exponential Radon transform and its adjoint, and only briefly mention some regularization techniques that could be useful. In particular, we include a description on how to include a discretized version of the filter operator present in the filtered back-projection type of reconstruction formula for the inversion of Radon data, as it is the most commonly used reconstruction method used for standard Radon data.
Let s ∈ R and θ ∈ S 1 , where S 1 denotes the unit circle. The standard twodimensional Radon transform is defined as and the parameters s and θ are used to parameterize the set of lines in R 2 . The attenuated Radon transform in R 2 describes the mapping from functions f and μ to line integrals of the form Integrals of this type appear for instance in medicine, where f describes the intensity distribution of isotopes or radiopharmaceuticals inside a tissue, and where the function μ(x) ≥ 0 describes the attenuation of the tissue. For the case where μ is constant within a known convex body, then the data obtained from the attenuated Radon transform can be modeled by the simpler exponential Radon transform [20], The adjoint operator associated with the exponential Radon transform (1) can be written as cf. [13,21]. It is a weighted integral of g over lines passing through the point x, and with μ = 0 this operator is typically referred to as the back-projection operator.
In [13] a generalization of the Fourier slice theorem is derived for the exponential Radon transform along with an integral equation for the reconstruction problem of obtaining f given R μ f . A similar integral equation is also derived in [5]. There also exist integral equations for the attenuated Radon transform, where the function μ has certain constraints. An explicit inversion formula is given by Novikov [15], and a generalization is given in [7] for a somewhat larger class of functions. In [14] a result closely related to Novikov's formula is derived, using a different approach. Yet another approach for inversion for arbitrary realistic attenuation was presented by Arbuzov et. al [4].
In [21] properties of the exponential Radon transform are derived, along with an inversion formula of filtered back-projection type. In Sect. 2 we recast this formula in a particularly simple format, see Theorem 2.1. Since then, the exponential Radon transform and the associated inversion problem have been discussed by several authors. For instance, in [1] range conditions of operator R μ are shown using Paley-Wiener type theorems. Cormack-type inversion formulas are based on the circular harmonic expansion of the transform and solving integral equations of special type (generalized Cormack equations) and they are discussed in [16] with generalization for the attenuated Radon transform in [17]. An explicit integral formula for 180-degree data reconstruction is obtained in [18] along with numerical tests for the approximation of the integrals. Another approach to reconstruct data from 180-degree measurements is proposed in [11], the approach is analytical except for the pre-calculated inverse kernel, where the least-squares method is utilized. In [2] the exponential Radon transform is extended to various spaces of distributions. In [22,23] it was shown that the exponential Radon transform inversion can be expressed as a convolution on Euclidean motion group (the rigid motion on R 2 ). Fourier-based approaches are proposed in [10,12]. Here, the Fourier series expansion in the angular variable allow to increase the stability of the inverse transform. In this paper we address a method for the exponential Radon transform, i.e., for the case where the attenuation parameter μ is constant.

The Fourier-Laplace Slice Theorem and Inversion in the Continuous Case
The Fourier slice theorem which relates the one-dimensional Fourier transforms of Radon transformed data with the two-dimensional Fourier transform of data is fairly straightforward to generalize to the exponential Radon transform. To see this, we introduce the transformR μ as the one-dimensional Fourier transform of the exponential Radon transform with respect to the first variable s of R μ f , i.e.
where the dot denotes scalar product in R 2 . Since f by assumption has support in the bounded domain , its Fourier transform extends to an analytic function in C 2 , and hence we can express the outcome of the above calculation aŝ which is a generalization of the Fourier slice theorem, sometimes referred to as the Fourier-Laplace slice theorem. For more details, see [13,21]. Note thatR μ maps compactly supported smooth functions to rapidly decaying smooth functions in L 2 (Z ) where Z = (R \ {0}) × S 1 is a polar-type coordinate system, except for the fact that we allow negative radii. We now wish to express this as a transform acting into a space of functions on Cartesian-type coordinates. We therefore let R 2 denote the two copies of the punctured plane on top of each other. Concretely, we set R 2 = (R 2 \ {0}) × {1, −1} and define a bijection ι : We will use the notationξ = (ξ, α) for the independent variable in R 2 , where ξ is the Cartesian part ofξ and α simply the sign of σ in the identityξ = ι(σ, θ). For a function g on Z or R 2 we will employ the standard convention of suppressing ι in the notation, writing g(σ, θ ) when the former (polar type) coordinate system is used, and g(ξ) when the latter (Cartesian type) coordinate system is used. Given any non-zero x = (x 1 , x 2 ) ∈ R 2 we introduce notation x ⊥ = (−x 2 ,x 1 ) |x| and moreover we defineξ Note that with this notation we always haveξ ⊥ = θ ⊥ wheneverξ = ι(σ, θ). Finally, we let H be the measure on R 2 that coincides with the Lebesgue measure on each of the two planes. For μ ≥ 0 set ThenR Using this notation, we shall show that the inversion formula by Tretiak and Metz takes a simple form given in the following theorem: Theorem 2.1 Let f be a compactly supported smooth function. Then for μ ≥ 0, we have or, expressed differently, Proof The inversion formula by Tretiak and Metz (see Corollary 1 of [21]) is written as an integral formula involving a family of kernels depending on two parameters σ, , and a limit as these approach zero. The reason for the parameters σ, is that the corresponding integrals for σ = = 0 may be divergent. For compactly supported smooth functions it is easy to see that this is not the case and the formula takes the following simpler form where W is a convolution (in s) with the inverse Fourier transform of the function w(σ ) = |σ/2|1 {|σ |> μ 2π } , and 1 {|σ |> μ 2π } denotes the characteristic function of {|σ | > μ 2π }, (the normalization of the Fourier transform is different in [21], so the kernel looks slightly different there). This leads to where M w is multiplication by w. Note that, for a smooth rapidly decaying function g on Z we have The formula (5) now follows by (4), (7) and the above computation. The subsequent formula is an immediate consequence since H equals Lebesgue measure on each copy of the punctured plane, andξ ⊥ = (ξ, α) ⊥ has opposite directions depending on whether α equals 1 or −1, so the evaluation of (5) reduces to By the formulas (4) and (5) it is clear that both R μ and the inversion formula can be evaluated by using Laplace transforms. To use the continuous formulas derived in this section, we will need to make use of discretization schemes, and this will be discussed in the next section.
To illustrate how the fast Laplace transform work, let us consider a one-dimensional discrete counterpart of the Laplace integrals above, and specifically the rapid evaluation of sums of the form Let ϕ be a Gaussian with a fixed width, and define ϕ a through the Fourier relation ϕ a (x) = e ax ϕ(x). It then holds for a single exponential that which can be expressed as By applying this relationship to the discrete Laplace sum above, we obtain that Due to the presence of the (modulated) Gaussians, the integrand above will be smooth and it can be approximated (to arbitrary precision) with the trapezoidal rule. This implies that it can be evaluated rapidly at equally spaced points in x by using FFT. To make this work in practice, the Gaussian ϕ must be comparatively wide in order to avoid loss of accuracy (it cannot be too close to zero). This means that the ϕ (and ϕ a j ) will have comparatively short support. An additional practical requirement is to use some amount of oversampling compared to the set of equally spaced point (x) at which the discrete Laplace transform is evaluated. We will make use of an oversampling factor ν to describe this level of oversampling. A typical value is ν = 2. Moreover, the inner sum in the expression above contains discrete convolution between the data values f j and the modulated Gaussians ϕ a j . As these Gaussians decay rapidly, only a few of terms will contribute (within given precision), and we may therefor truncate the Gaussians. We will use the parameter M to denote where we can truncate the Gaussian. Specifically, this means that the discrete convolution above is performed over at 2M + 1 terms, and the contribution from the other terms can be discarded. The width parameter M (and the choice of the Guassian ϕ) is related to the desired level of accuracy as well as the oversampling parameter ν. For more details, see Appendix.

Discretization
We will now consider discretized versions of the exponential Radon tranform and its adjoint. The discretization will be fairly standard, using an equally spaced sampling in the spatial coordinates as well as in the (s, θ)-coordinates. We assume that f has compact support. Without loss of generality we assume that it has support on the disc {x : |x| < 1/2}, and we represent it by sampling on the Cartesian lattice for N even, and let III X denote the operator that samples functions on X . We also discretize R μ f (s, θ) at equally spaced points in θ that cover the whole interval [0, 2π). Specifically, let us assume that we have samples at a grid with N number of points in the s-direction and N θ points in the θ -direction and define Via the discrete Fourier transform, the nodes s m are associated with frequency nodes as well as a corresponding polar grid, We let III S and III denote the sampling operators on the two sets. As us customary in this field, we identify θ l with the point (cos(θ l ), sin(θ l )) on S 1 , wherever convenient.

Fast Computation of the Discretized R µ
We now compute a discretized version of the exponential Radon transform (1) in two steps using two-dimensional Laplace transforms [3] and FFT. A discrete counterpart If μ = 0 then sums of the above type can be approximately evaluated (with precision ) in O(N 2 log N + N N θ log(1/ )) by using unequally spaced fast Fourier transforms (USFFT): [6,8,9]. To account for the case where μ = 0, we will instead use unequally spaced fast Laplace transforms [3] (USFLT), at the same computational cost and the same precision dependence as for the USFFT. We will denote the corresponding operator L USFLT d . See Appendix for more information on this topic in two dimensions.
We denote the one-dimensional discrete inverse Fourier transform acting on the first variable by F * d : C → C S , i.e.
Summing up, we can approximately compute III S R μ f by F * d L USFLT d III X f . An analysis of the errors involved in this approximation is found in Appendix.

Error Analysis
The total error of our discretization/approximation is which by the triangle inequality and the fact that DFT is isometric can be estimated by the sum of the error caused in the discretization of the one-dimensional Fourier integrals, the discretization error for the Laplace transform and the error due to the fast evaluation of the discrete Laplace transform The first two errors are due to the discretization, and the third one is due to the proposed fast computational algorithm. As mention above, the error in (13) will be bounded by the choice of precision parameter , see Appendix. We will now show that the errors in (11) and (12) bear similarities to standard sampling results about essentially bandlimited functions.
To this end we introduce as a measurement of the frequency content of f outside [−N /2, N /2] 2 . Note that for bandlimited functions whose frequency support is contained in the central square with side length N it holds that B 0,N ( f ) = 0.
Proof By (3) it follows that Set g(σ ) = f (σ θ + μθ ⊥ ) and note that its frequency support lies within (−1/2, 1/2). The Poisson summation formula yields where the summation from −N /2 to N /2 − 1 corresponds to F * d . It is now easy to see that the difference (14) is bounded by B μ,N ( f ), yielding an ∞ estimate of (11). The corresponding 2 estimate follows by the equivalence of finite dimensional norms.
The estimate for (12) is similar. By the two-dimensional Poisson summation formula (applied to f (x)e −μθ ⊥ l ·x ) we have By the assumption on the support of f , the first double sum can be replaced by which yields a pointwise estimate of the difference in (12). The desired conclusion again follows by the equivalence of finite dimensional norms.

Fast Computation of the Discretized Adjoint
The back-projection operator R * μ in (2) can be discretized in a similar way. Note that which follows from (2) and the fact that g(s, θ )).
In the discretized setting, the inner integral (i.e. F s→σ ) corresponds to F d whereas the two outer integrals (i.e.R * μ ) correspond to (L USFLT d ) * . We remark that the latter do not correspond to a two-dimensional Fourier-Laplace integral (since the factor |σ | is missing), but it can still be evaluated fast, i.e., in O(N 2 log N ) time. In the terminology of the Appendix, the USFLT operation changes from a Z 2 → C 2 operation (for L USFLT d ) to a C 2 → Z 2 operation (for (L USFLT d ) * ), designed to rapidly evaluate sums of the type for all values n 1 N , n 2 N in X , where F k,l can be arbitrary numbers. For the evaluation of III S R μ * g via (L USFLT d ) * F d III S g, we obviously have F = F d III S g, but we shall see in the next section that other choices also will be relevant. An error analysis of this approximation can be done in a similar way as in Sect. 3.2, but we omit this.

Inversion from Complete Measurements
To invert measurements of the exponential Radon transform on the polar grid S, we use the formula (7); A discrete approximation g ofR μ f on is easily computed by applying F d to the available data III S R μ ( f ), and the operationR * −μ corresponds to (L USFLT d ) * . Since M w is multiplication by the weight w, the inversion could in principle be done via but the discontinuities of w will introduce large errors for low-frequency components of f . To compensate for this, we will evaluateR μ f on a polar grid˜ with a more dense radial sampling near the singularities σ = ± μ 2π . The operator (L USFLT d ) * computes sums of the type (16) over unequally spaced grids, and in particular there is no need to have equally spaced sampling in σ . However, if we increase the sampling near the singularity, the corresponding sum (16) will not be a good approximation of the corresponding integral inR * −μ , and we therefore need to compensate for this by modifying the weight w, using composite rules for quadratures. We briefly outline this in more detail below.
In practice, the bottleneck of the reconstruction algorithm lies in evaluating (L USFLT d ) * , and hence it is desirable to keep the total number of grid points low. We therefore propose to use a dense sampling in σ near the singularities ± μ 2π . We remark that values of μ 2π typically are less than 2, and that the original grid is sampled with radial sampling distance of unit length. An increased sampling near μ 2π can thus be achieved by increasing the amount of grid points in near the origin. Recall that f is assumed to have support on the unit circle, which implies that R μ ( f )(s, θ) is supported on [−1/2, 1/2] in the s-variable. By zero-padding of R μ ( f )(s, θ) in the s-direction and application of FFT, an oversampled representation of g(σ, θ ) can be obtained (on an equally spaced grid in σ , but not all values will be used).
To derive appropriate weights, it is sufficient to consider accurate quadratures for where we fix θ and set h(σ ) = 1 2 g(σ, θ )e −μx·θ ⊥ +2πiσ x·θ . Concerning the region around the singularity σ = μ 2π , we construct a composite quadrature by considering In this way, we obtain quadrature weights w k for the equally spaced approximation of (18). The weights can be combined to form composite quadrature rules that take the singularity at σ = μ 2π into account. The same procedure can be performed when dealing with the singularity at σ = − μ 2π . The left panels of Fig. 1 demonstrate the effect of the correction with the weights w k near both the singularities.
In order for the approximation step in (19) to be accurate, h has to be sufficiently oversampled. However, this is only needed in a vicinity of the singularity, and a regular sampling is sufficient away from the singularity. A direct usage of a trapezoidal rule would, however, use a sampling that is determined by the needed sampling rate around the singularity. We therefore split the integral as Now, the derived quadrature rule with a dense sampling can be used for the left integral in the right-hand side above, while a reduced sampling density can be used on the right integral. This procedure can also be used several times to gradually decrease the sampling rate. In the numerical simulations presented in Sect. 4, we use an oversampling factor of four to account for the singularity at μ 2π , and make sampling reduction in two steps, i.e., by a quadrature rule of the form where K 3 ≤ K 2 << K 1 . The setup is illustrated in the right panels of Fig. 1.

Numerical Simulations
In this section reconstruction results and computational speed tests are presented. C++ routines were written for the fast exponential Radon transform and the fast Laplace transform along with MEX interfaces to MATLAB. For accelerations we used the tools from Intel Composer XE 2016, including the Intel Math Kernel Library (MKL), Intel Performance Primitives (IPP), and OpenMP API. For the simulation we used a standard desktop with an Intel i7-3820K processor and eight OpenMP threads, and the computations were performed in single precision.
To begin with, accuracy tests were performed on the Shepp-Logan [19] phantom, where we use the modified version available in MATLAB. The phantom consists of linear combinations of characteristic functions of ellipses, and its support is inscribed in a circle, see the left panel of Fig. 2. The red line is depicted with correspondence to a specific angle θ and polar sample s. Values through this line multiplied by factor e μt are plotted in the left panel of Fig. 2. The upper part of Fig. 3 show the exponential Radon transform for different values of the exponential parameter μ.
In order to make accuracy comparisons, we do not carry out comparisons with the original function, but use instead a slightly smoothed out version of the phantom. This is to account for the fact that the characteristic function of the ellipses does not decay rapidly in the frequency domain, and hence high-frequency errors would be dominant in the error plots. A reconstruction method using only a finite frequency range should have as aim to recover that frequency range accurately, and it is therefore natural to measure the accuracy in this way. The minor smoothing effect can be noted in the plots of Fig. 2.
The lower part of Fig. 3 illustrates reconstruction errors for different values of μ chosen to be the same as for the corresponding projection from the upper part of the same figure. To show how the choice of μ affects the reconstruction we compute reconstructions using the fast Laplace based filtered back-projection. In the filtered back-projection method the high-frequency components are boosted, and boosting of high frequent noise is also caused by the exponential factor in the back-projection Fig. 2 Weighting with factor e μt for a line through the phantom image. Exponential weight y(t) = e μt is plotted by the dashed line, ratios between maximal and minimal weights equal to 1, 10, 100, 1000, respectively  operator. In particular for high attenuation parameters, this reconstruction technique becomes sensitive to noise. This is illustrated in Fig. 4, which shows filtered-backprojection reconstruction results for varying values of the attenuation parameter μ. The Radon data for these different attenuation values are contaminated with white Gaussian noise with a signal-to-noise ratio (SNR) of 20 dB.
Let us now turn our attention to the computational speed of the reconstruction method. Tables 1 and 2 show the total times of the simulations and the times per operation (to verify the time complexity) of applying the exponential Radon transform and backprojection, respectively, for different image sizes.

Conclusions
We have shown how to construct fast algorithms for the computation of the exponential Radon transform and the associated (adjoint) back-projection operator by using algorithms for fast Laplace transforms. In addition, we have included estimates for the discretization errors that arise, and shown how to separate them into parts where one bears similarities to standard sampling arguments of (almost) bandlimited functions, and where the other part is due to the approximation errors arising in the fast Laplace transform. The latter part can be controlled at arbitrary numerical precision. and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Computation of L USFLT d
For the sake of completeness, we present algorithms for fast evaluation of twodimensional discrete Laplace transforms. These algorithms are analogous to the one-dimensional ones presented in [3], but in this appendix they are adapted for computing the exponential Radon transforms. Let n 1 , n 2 andξ j be as in Sect. 3. We describe how to fast evaluate the operations where f n 1 ,n 2 and F j represent any complex numbers. Note that with f n 1 , , the operation (22) corresponds to the evaluation of (9) in Sect. 3, whereas (21) is needed in (16). Methods for unequally spaced fast Laplace transforms (USFLT) utilize convolution style operations with FFT to achieve approximations of the sums (21,22) with arbitrary, but fixed, accuracy. We make a heuristic description of how USFLT works, for a more precise formulation along with error estimates, see [3].
Let F denote the Fourier transform and introduce the Gaussian ϕ(x) = e −α x 2 1 +x 2 2 and modulated Gaussian ϕξ⊥ where α is a parameter which we discuss how to choose below. Note that in the frequency domain it holds that i.e. ϕξ⊥ j is a frequency modulated Gaussian. Next, by choosing g = f ϕ one can derive the relation The integral in (23) can be approximated by the trapezoidal rule because of the fact that ϕξ⊥ j are smooth functions that decay rapidly to zero (the amount of oscillation is limited by the parameter μ). It turns out that accurate results can be obtained by (v 1 , v 2 ) on an equally spaced lattice with some oversampling parameter ν. Crucial for the construction of fast algorithms is that the function ϕξ ⊥ has small (numerical) support, since this will keep down the computational times for convolutions with ϕξ ⊥ . It turns out that the (numerical) width of ϕξ ⊥ is proportional to − log , where is the desired precision. More precisely, values of ϕξ⊥ j that are smaller than a certain threshold level can be omitted. In practice we thus replace ϕξ⊥ Now, if we assume that the discrete Fourier transform (realized by FFT) provides an accurate approximation ofĝ at this oversampled grid, the evaluation of (22) can be approximated in three steps: 1. Division in the space domain 2. FFT 3. Convolution-type operation in the frequency domain We refer to [3] for proofs and more details, and remark that μ = log A in the notation of that paper. The approximations to prescribed precision can be achieved by a suitable choice of α in combination with proper oversampling ν. It then holds that the difference between the sum of (21) and that described by the above algorithm can be bounded for instance by where C is a constant that is independent of f . A similar result can also be obtained for the approximation of (22). Following [3], the parameter α should be chosen as where ν is an oversampling factor for Fourier transform (usually ν = 2). In practice this means that the integral (23) is replaced by a summation over k 1 , k 2 ∈ Z using ϕ ξ ⊥ j k 1 ν − ξ 1 j , k 2 ν − ξ 2 j . Note that for each j = 1 . . . J , these values will be nonzero only for values of k 1 , k 2 such that |k 1 − νξ 1 j | ≤ νT and |k 2 − νξ 2 j | ≤ νT . By introducing M = νT we find that the contributing values of k 1 , k 2 satisfy Finally, we have an approximation formula for (22) with accuracy : where the sums over k 1 , k 2 variables consist each of (2M + 1) non-zero terms, respectively. The constant M is typically equal to 8-10 for single precision of computations, depending on the choice of μ and ν. Algorithm L Z 2 →C 2 ( f ) describe these computational steps, with the assumption |ξ j | are inside a circle with radius slightly smaller than N 2 . The first operation (division) has O(N 2 ) computation cost, the second Algorithm L C 2 →Z 2 .
3. f n 1 ,n 2 = h n 1 ,n 2 ϕ 0 ( operation is evaluated by FFT in O(ν 2 N 2 log(N )) time, and the third operation cost O((2M + 1) 2 N N θ ) operations, since for each j there are only (2M + 1) 2 non-zero contributions. The total time complexity is thus O(ν 2 N 2 log(N ) + (2M + 1) 2 N N θ ). In terms of time complexity, this is dominated by FFT part, while in practice most time is typically spent on the third step. The algorithm for L C 2 →Z 2 ( f ) (21) can be also constructed by applying the above operations in reverse order, since L C 2 →Z 2 ( f ) is adjoint to L Z 2 →C 2 ( f ). These steps are gathered in Algorithm L Z 2 →C 2 below.