Saturation Rates of Filtered Back Projection Approximations

This paper concerns the approximation of bivariate functions by using the wellknown filtered back projection (FBP) formula from computerized tomography. We prove error estimates and convergence rates for the FBP approximation of target functions from Sobolev spaces Hα(R2) of fractional order α > 0, where we bound the FBP approximation error, which is incurred by the application of a low-pass filter, with respect to the weaker norms of the rougher Sobolev spaces Hσ (R2), for 0≤ σ ≤ α . In particular, we generalize our previous results to non band-limited filter functions and show that the decay rate of the error saturates at fractional order depending on smoothness properties of the filter’s window function at the origin. The theoretical results are supported by numerical simulations.


Introduction
The method of filtered back projection (FBP) is a commonly used reconstruction technique in computerized tomography (CT), which deals with recovering the interior structure of an unknown object from given X-ray scans, where the measured X-ray data is interpreted as a set of line integrals for a bivariate function, termed attenuation function. To formulate the classical CT reconstruction problem mathematically, we regard for f ≡ f (x, y) its Radon transform R f ≡ R f (t, θ ), given by With this definition, the basic CT reconstruction problem can be formulated as follows.
Problem 1 (Reconstruction problem) On domain Ω ⊆ R 2 , reconstruct a bivariate function f ∈ L 1 (Ω ) from given Radon data Hence, the CT reconstruction problem seeks for the inversion of the Radon transform R. We remark that this problem has a very long history, dating back to Johann Radon, whose seminal work [17] provided an analytical inversion of R already in 1917. This has later led to the filtered back projection (FBP) formula (cf. [6,13]), where F is the univariate Fourier transform acting on the radial variable, i.e., for g ≡ g(t, θ ) h(x cos(θ ) + y sin(θ ), θ ) dθ for (x, y) ∈ R 2 .
In the following of this paper, we explain on (1) and its ingredients in more detail. For the moment, we only wish to remark that the FBP formula (1) is highly sensitive with respect to noise, due to the filter |S|. In standard stabilization methods for the FBP formula, |S| in (1) is replaced by a compactly supported low-pass filter A L : R −→ R, with bandwidth L > 0, of the form where W ∈ L ∞ (R) is an even window function of compact support supp(W ) ⊆ [−1, 1], cf. [14]. This modification leads us to an approximate FBP reconstruction formula Examples for window functions W of commonly used low-pass filters are shown in Table 1. However, the limitation to band-limited low-pass filters is not required, see e.g. [23]. Although the FBP method has been one of the standard reconstruction algorithms in CT for decades, its error analysis and convergence behaviour are not completely settled so far. In [15], Popov showed pointwise convergence restricted to a small class of piecewise smooth functions. The approach of Rieder and Schuster [19] leads to L 2 -convergence with suboptimal rates for compactly supported Sobolev functions. In contrast to this, in [18,20] Rieder et al. prove optimal L 2 -convergence rates for sufficiently smooth Sobolev functions.  However, the authors verify their assumptions only for a restricted class of filters based on B-splines. More recently, Qu [16] showed convergence without rates in the L 2 -norm for compactly supported L ∞ -functions and in points of continuity under additional assumptions. Note that [16] deals with the continuous problem, while [18,19,20] discuss discrete settings.
In this paper, we also focus on the continuous setting and analyse the inherent FBP reconstruction error e L = f − f L of the FBP approximation f L that is incurred by the application of the low-pass filter A L .
To this end, we prove novel error estimates in Sobolev spaces of fractional order and provide quantitative criteria to evaluate the performance of the utilized filter by means of its window function, which is not required to have compact support. Further, we prove convergence for the approximate FBP reconstruction in the treated Sobolev norms along with asymptotic convergence rates as the filter's bandwidth goes to infinity, where the smoothness order of the target function f is only required to be positive. Most notably, our results allow us to predict saturation of the order of convergence at fractional rates depending on smoothness properties of the filter's window function, which can be easily evaluated. The results presented in this paper generalize our previous findings in [1,2,3] and confirm the key observation that the flatness of the filter's window function at the origin determines the convergence behaviour of the approximate FBP reconstruction. The outline of this paper is as follows. In §2 we explain how low-pass filters can be used to stabilize the FBP formula (1) yielding an approximate inversion formula for bivariate functions. We analyse the approximation error of the FBP reconstruction in §3 and derive Sobolev error estimates for target functions from fractional Sobolev spaces. Assuming certain regularity of the filter's window, we prove asymptotic convergence rates as the bandwidth goes to infinity in §4, where the presented results generalize our previous work [1,2,3] to more general filters and show that the rate of convergence saturates at fractional order depending on smoothness properties of the window. We finally provide numerical experiments in §5 to support our theoretical results and conclude with summarizing remarks in §6.

Approximate FBP reconstruction formula
In this section we rigorously define the approximate FBP reconstruction f L in (3), i.e., where we more generally assume that the even window W : R −→ R satisfies W ∈ L ∞ (R) and | · |W (·) ∈ L 2 (R).
(W) In particular, we show that f L is defined almost everywhere on R 2 , for f ∈ L 1 (R 2 ), and can be rewritten as where * denotes the convolution product of bivariate functions in polar coordinates, given by and where we define the convolving function q L : Moreover, for f ∈ L 1 (R 2 ) ∩ L 2 (R 2 ) we show f L ∈ L 2 (R 2 ) and express f L directly in terms of the target function f via where the convolution kernel K L : R 2 −→ R is defined as Before we proceed with the derivation of the representations (4) and (6), for the reader's convenience we briefly collect a few relevant facts concerning the Radon transform R that we need in our subsequent analysis. Since the following results are well-known, we omit the proofs and refer to the literature instead, e.g. [7,13,22]. .
Moreover, if f is compactly supported, then R f has compact support as well.
Next we recall that the L 2 -norm of R f is bounded, provided that the function f belongs to L 2 c (R 2 ), i.e., f is square integrable and has compact support.
Lemma 2 indicates that R can be viewed as a densely defined unbounded linear operator from L 2 (R 2 ) to L 2 (R × [0, π)) with domain L 2 c (R 2 ). We now turn to its adjoint operator R # .
According to Lemma 3 the back projection B is, up to constant 1 π , the adjoint of R. In particular, for g ∈ L 2 (R × [0, π)) the function Bg is defined almost everywhere on R 2 and satisfies Bg ∈ L 2 loc (R 2 ).
We are now prepared to show that the approximate FBP reconstruction f L in (3) is defined almost everywhere on R 2 , for f ∈ L 1 (R 2 ), and can be represented as in (4).
Note that the convolving function q L = F −1 A L in (5) satisfies q L ∈ L 2 (R × [0, π)), since the filter A L in (2) is in L 2 (R) due to Assumption (W). Thus, the convolution kernel K L in (7) is defined almost everywhere on R 2 and satisfies K L ∈ L 2 loc (R 2 ) according to Lemma 3. In the following, we prove K L ∈ L 2 (R 2 ) and determine the bivariate Fourier transform of K L , as needed in the upcoming analysis of the FBP approximation error. To this end, we extend the window function W L = W ( · /L) to R 2 by its radialization W L : R 2 −→ R, i.e., where we let r(x, y) = x 2 + y 2 for (x, y) ∈ R 2 .
Proposition 2 Let W ∈ L ∞ (R) be even with | · |W (·) ∈ L 2 (R). Then, for any L > 0 the convolution kernel K L in (7) satisfies K L ∈ L 2 (R 2 ) and its Fourier transform is given by Proof Since W ∈ L ∞ (R) satisfies | · |W (·) ∈ L 2 (R), the bivariate window W L also satisfies W L ∈ L 2 (R 2 ) ∩ L ∞ (R 2 ) and F −1 W L ∈ L 2 (R 2 ) due to the Rayleigh-Plancherel theorem. Moreover, in L 2 -sense and, thus, for almost all (x, y) ∈ R 2 follows that A L (S) e iS(x cos(θ )+y sin(θ )) dS dθ by transforming (X,Y ) = (S cos(θ ), S sin(θ )) from Cartesian to polar coordinates. Recall that we have q L ∈ L 2 (R × [0, π)) and, for all θ ∈ [0, π), Thus, with the definition of the back projection operator B follows that Consequently, we have K L ∈ L 2 (R 2 ) and applying the Fourier inversion formula shows that holds in L 2 -sense and, in particular, almost everywhere on R 2 .
Before we proceed, we wish to add one remark concerning the convolution kernel K L . If the even window function W ∈ L ∞ (R) in addition to Assumption (W) satisfies | · |W (·) ∈ L 1 (R), the bivariate window W L ∈ L ∞ (R 2 ) is also in L 1 (R 2 ) and we have K L ∈ C 0 (R 2 ) ∩ L 2 (R 2 ) with We can now prove the desired representation in (6) in the L 2 -sense. In particular, we show that f L ∈ L 2 (R 2 ) and determine its Fourier transform as (8).
which in turn implies that Recall that the identity holds for almost all S ∈ R and fixed θ ∈ [0, π). Hence, we obtain Since (q L * R f )(·, θ ) ∈ L 2 (R) for all θ ∈ [0, π), the Fourier inversion formula holds again in L 2 -sense and, in particular, almost everywhere on R. Thus, for fixed θ ∈ [0, π) we get and the definition of the back projection B in combination with Proposition 1 yields for almost all (x, y) ∈ R 2 . This finally implies f L = f * K L ∈ L 2 (R 2 ).
Combining Propositions 2 and 3 allows us to determine the Fourier transform F f L of the approximate FBP reconstruction f L .
Then, for all L > 0 the Fourier transform F f L of the approximate FBP reconstruction f L is given by in L 2 -sense and, in particular, almost everywhere on R 2 .
We wish to add the following remark on the regularity of the FBP approximation f L . Since the target function f and the kernel K L are both in L 2 (R 2 ), the representation (6) implies that f L is continuous and vanishing at infinity, i.e., we have f L ∈ L 2 (R 2 ) ∩ C 0 (R 2 ).

Recall that we typically deal with even window functions
is automatically satisfied and the assumption f ∈ L 2 (R 2 ) can be omitted. Moreover, Corollary 1 shows that in this case the approximate FBP reconstruction formula (4) provides a band-limited approximation f L to the target function f . In particular, the FBP approximation f L is arbitrarily smooth. (7). Moreover, its Fourier transform is given by To close this section, we remark that the approach taken in this work is standard in the mathematics of computerized tomography and especially the representations (4) and (6) of the approximate FBP reconstruction f L are known from the literature, see e.g. [8,23]. However, we prove the statements under relatively mild conditions directly posed on the filter's window W and not on the corresponding convolution kernel K, unlike in [8,23].

Analysis of the inherent FBP reconstruction error
In the previous section we have seen that the application of a low-pass filter A L of the form with finite bandwidth L > 0 and an even window W ∈ L ∞ (R) satisfying | · |W (·) ∈ L 2 (R) leads to an approximate FPB reconstruction f L that can be expressed as For the sake of brevity, we call the application of the approximate FBP formula (9) an FBP method. Therefore, each FBP method provides one approximation f L to f , f L ≈ f , whose quality depends on the choice of the low-pass filter A L . We remark that in practical situations the target function f belongs to L 2 (R 2 ) and has compact support such that its Radon transform satisfies R f ∈ L 2 (R × [0, π)). However, the Radon data g = R f is usually not known precisely, but only up to an error bound δ > 0, and we have to reconstruct f from given data Applying the FBP method (9) to the noisy measurements g δ now yields the reconstruction and, using standard concepts from the theory of inverse problems and regularization (cf. [4]), we observe that the total reconstruction error can be split into an approximation error term and a data error term: For an analysis of the data error in the L 2 -norm in the case of band-limited low-pass filters we refer to [2,Section 8]. In this paper, we focus on the approximation error, i.e., we analyse the inherent error of the FBP method which is incurred by the application of the low-pass filter A L . More precisely, we wish to analyse the reconstruction error with respect to the filter's window function W and bandwidth L. We remark that pointwise and L ∞ -error estimates on e L in (10) are discussed by Munshi et al. in [9,10]. Their results are further supported by numerical experiments in [11]. Error bounds on the L p -norm of e L , in terms of L p -moduli of continuity of the target function f , are proven by Madych in [8]. But our approach is essentially different from previous approaches, see also [2, Section 3].
We prove Sobolev error estimates on e L for target functions f from Sobolev spaces and S (R 2 ) denotes the space of tempered distributions on R 2 . In relevant applications of (medical) image processing, Sobolev spaces of compactly supported functions, on an open and bounded domain Ω ⊂ R 2 , and of fractional order α > 0 play an important role (cf. [12]). In fact, the density function of an image in Ω ⊂ R 2 has usually jumps along smooth curves, but is otherwise smooth off these curve singularities. Such functions belong to any Sobolev space H α 0 (Ω ) with α < 1 2 and we can consider the density of an image as a function in a Sobolev space H α 0 (Ω ) whose order α is close to 1 2 . Consequently, throughout this section, we assume that we deal with target functions f ∈ L 1 (R 2 ) ∩ H α (R 2 ) for some α > 0 and that we are given some bandwidth L > 0 and an even window W ∈ L ∞ (R) satisfying | · |W (·) ∈ L 2 (R). Under these assumptions, we now prove Sobolev error estimates for the FBP reconstruction error e L in (10) with respect to the H σ -norm, for all 0 ≤ σ ≤ α. This in particular gives L 2 -error estimates, when σ = 0. Moreover, this generalizes our previous results in [1,2,3], where we more restrictively deal with band-limited low-pass filters, i.e., compactly supported window functions.
We start with showing that the FBP approximation f L belongs to the Sobolev space H σ (R 2 ) for 0 ≤ σ ≤ α. To this end, recall that we have f L ∈ L 2 (R 2 ) with F f L = W L · F f . Thus, the H σ -norm of f L can be bounded above by Let us now analyse the inherent FBP reconstruction error e L = f − f L with respect to the H σ -norm. For γ ≥ 0, we define so that, for any ε > 0, the H σ -norm of e L can be expressed as where we let For γ ≥ 0, we define the function so that we can bound the integral I ε 1 in (11) from above by In addition, for 0 ≤ σ ≤ α, we can bound the integral I ε 2 in (12) by Combining the estimates for I ε 1 and I ε 2 , we finally obtain so that we can summarize the discussion of this subsection as follows.
be even with | · |W (·) ∈ L 2 (R). Then, for 0 ≤ σ ≤ α and any ε > 0, the H σ -norm of the FBP reconstruction error e L = f − f L is bounded above by where We now prove that, under suitable assumptions on the window W , the approximate FBP reconstruction f L converges to the target function f as the bandwidth L goes to infinity.
In [1, Corollary 1] we have seen that for band-limited low-pass filters and functions f ∈ L 1 (R 2 ) ∩ H α (R 2 ) with α > 0 the H σ -norm of the FBP reconstruction error e L = f − f L converges to 0 for L → ∞, as long as 0 ≤ σ < α. Additionally, we had to assume that the even window function W is continuous on the interval [−1, 1] and satisfies W (0) = 1. In the following, we relax these assumptions and proof convergence for the H σ -norm of e L for all 0 ≤ σ ≤ α, where especially the case σ = α is included. However, the proof technique is not suitable for determining the rate of convergence.
with α ≥ 0 and let W ∈ L ∞ (R) be even with | · |W (·) ∈ L 2 (R). Further, let W be continuous at 0 with W (0) = 1. Then, for 0 ≤ σ ≤ α, the H σ -norm of the FBP reconstruction error e L = f − f L converges to 0 as L goes to ∞, i.e., we have seen that the FBP approximation f L belongs to the Sobolev space H σ (R 2 ) for all 0 ≤ σ ≤ α and that the H σ -norm of the FBP reconstruction error e L = f − f L is given by For all L > 0 and almost all (x, y) ∈ R 2 holds that Furthermore, for all (x, y) ∈ R 2 holds that since the window W is continuous at 0 and satisfies W (0) = 1. Therefore, we can apply Lebesgue's theorem on dominated convergence and obtain i.e., the FBP approximation f L converges to f in the H σ -norm for all 0 ≤ σ ≤ α.
As a corollary we obtain the convergence of the FBP reconstruction f L in L 2 (R 2 ), which is also stated in [16, Theorem 2.2] under stronger assumptions.

Convergence rates for approximate FBP reconstructions
In this section, we analyse the rate of convergence of the reconstruction error e L = f − f L in the H σ -norm for 0 ≤ σ ≤ α as the bandwidth L goes to ∞. To this end, let S * γ,W,ε (L) ∈ [0, ε], for γ ≥ 0 and ε > 0, denote the smallest maximizer in [0, ε] of the even function i.e., To determine the rate of convergence for e L σ , we assume that S * α−σ ,W,ε (L) is uniformly bounded away from 0 for 0 ≤ σ ≤ α, i.e., there exists a constant c α−σ ,W,ε > 0 such that Then, the error term Φ α−σ ,W,ε (L) is bounded above by Consequently, under Assumption (A) we obtain In summary, this yields the following result.
Theorem 3 (Rate of convergence in H σ ) Let f ∈ L 1 (R 2 ) ∩ H α (R 2 ) for α > 0 and let W ∈ L ∞ (R) be even with | · |W (·) ∈ L 2 (R). Further, let Assumption (A) be satisfied. Then, for 0 ≤ σ ≤ α and any ε > 0, the H σ -norm of the FBP reconstruction error e L = f − f L is bounded above by In particular, Note that the decay rate in (14) is determined by the difference between the smoothness α of the target function and the order σ of the Sobolev norm in which the error is measured.
Moreover, we remark that Assumption (A) is fulfilled for a large class of windows. For instance, let W ∈ L ∞ (R) satisfy W (S) = 1 for all S ∈ (−η, η) with η ∈ (0, ε]. Then, Assumption (A) is fulfilled with c α−σ ,W,ε = η for all 0 ≤ σ ≤ α and the H σ -error estimate (14) reads For ε = η, this reduces to In particular, Assumption (A) is satisfied for the classical Ram-Lak filter, see Table 1, in which case we obtain f − f L σ ≤ L −(α−σ ) f α and the decay rate of the FBP reconstruction error is determined by the difference between the smoothness α of the target function and the order σ of the considered Sobolev norm. However, Assumption (A) is not always satisfied for the other commonly used low-pass filters listed in Table 1 and, as we will see, the decay rate of the error can saturate depending on smoothness properties of the filter's window function at the origin.
Error estimates for A C -windows with L p -derivatives In the following, we consider the special case of even window functions W ∈ L ∞ (R) with | · |W (·) ∈ L 2 (R) that are absolutely continuous on

ε]) and the Fundamental Theorem of Calculus for A C -functions yields
Further, for some k ∈ N we assume that W satisfies and W (k) ∈ L p ([−ε, ε]) for some 1 < p ≤ ∞. Under these assumptions, we will prove that, for 0 ≤ σ ≤ α, the H σ -norm of the FBP reconstruction error e L = f − f L now behaves like According to the H σ -error estimate (13) from Theorem 1 it suffices to analyse the error term where the parameter γ ≥ 0 has to be chosen as γ = α − σ . To this end, for ν > 0 we consider the auxiliary function Lemma 4 Let ε > 0. The maximum of the function φ γ,L,ν in (15) on [0, ε] is bounded above by γ−ν and the strictly monotonically decreasing constant Proof The even and continuous function φ γ,L,ν in (15)  The first derivative of φ γ,L,ν (S) is given by so that the necessary condition for a maximum of φ γ,L,ν in (0, ε) reads Case 1: For 0 ≤ γ ≤ ν, equation (16) has no solution in (0, ε] and, actually, This shows that φ γ,L,ν is strictly monotonically increasing in (0, ε] and, thus, Case 2: For γ > ν, the unique positive solution of equation (16) is given by Moreover, φ γ,L,ν (S) < φ γ,L,ν (S * ) = 0 < φ γ,L,ν (s) ∀ 0 < s < S * < S so that φ γ,L,ν is strictly increasing on (0, S * ) and strictly decreasing on (S * , ∞). Because of for L ≥ L * ε . Let us finally regard the constant c γ,ν as a function c ν ≡ c ν (γ) of the parameter γ > ν, i.e., For the first derivative of c ν follows that and, consequently, the constant c γ,ν is strictly monotonically decreasing in γ > ν.
We are now prepared to derive H σ -error estimates for A C -windows with L p -derivatives and start with the case p = ∞, i.e., we assume that W (k) ∈ L ∞ ([−ε, ε]).
Theorem 4 (H σ -error estimate for A C -windows with L ∞ -derivatives) For α > 0, let f ∈ L 1 (R 2 ) ∩ H α (R 2 ) and, for k ∈ N and ε > 0, let W satisfy W Then, for 0 ≤ σ ≤ α, the H σ -norm of the inherent FBP reconstruction error e L = f − f L is bounded above by for α − σ > k and sufficiently large L > 0 with the strictly decreasing constant In particular, Proof Based on our assumptions and on the H σ -error estimate (13) from Theorem 1, i.e., it is sufficient to analyse the error term For k ≥ 2, we can apply integration by parts and get, for all S ∈ [0, ε], By iteratively applying integration by parts we obtain Hence, for any k ∈ N we have and, thus, the error term Φ α−σ ,W,ε (L) can be bounded above by and the strictly decreasing constant Combining the estimates completes the proof.
Note that the convergence rate of e L σ in Theorem 4 is determined by the difference between the smoothness α of the target function f and the order σ of the Sobolev norm in which the reconstruction error e L is measured, as long as α − σ ≤ k. But for α − σ > k the rate of convergence saturates at integer order O(L −k ). However, in this case the error bound still decreases at increasing α, since the constant c α−σ ,k is strictly decreasing in α − σ > k.
Finally, we remark that Theorem 4 shows that our estimates from [1, Theorem 3] and [2, Corollary 7.2] continue to hold under relaxed assumptions. In [1,2] we have proven these statements for band-limited filters under the more restrictive assumption W ∈ C k ([−1, 1]).
Theorem 5 (H σ -error estimate for A C -windows with L p -derivatives) For α > 0, let f ∈ L 1 (R 2 ) ∩ H α (R 2 ) and, for k ∈ N and ε > 0, let W satisfy Then, for 0 ≤ σ ≤ α, the H σ -norm of the inherent FBP reconstruction error e L = f − f L is bounded above by Proof Based on our assumptions and on the H σ -error estimate (13) from Theorem 1, i.e., it is sufficient to analyse the error term As in the proof of Theorem 4, we (iteratively) apply integration by parts and obtain and, thus, it follows that Consequently, the error term Φ α−σ ,W,ε (L) can be bounded above by Applying Lemma 4 shows that and strictly decreasing constant Combining the estimates completes the proof.
We observe that the convergence rate of the error bound in Theorem 5 is determined by the difference between the smoothness α of the target function and the order σ of the Sobolev norm in which the reconstruction error e L is measured, as long as α − σ ≤ k − 1 /p. But for α − σ > k − 1 /p the rate of convergence saturates at fractional order O(L −(k−1/p) ). However, in this case the error bound still decreases at increasing α, since the involved constant c α−σ ,k,p is strictly decreasing in α − σ > k − 1 /p. Thus, although the decay rate saturates, a smoother target function still allows for a better approximation, as expected.
We finally remark that the error estimate from Theorem 5 is consistent with the error estimate from Theorem 4 for windows W with W (k) ∈ L ∞ ([−ε, ε]). Indeed, we have as well as 1 (k − 1)! To close this section, we give an example of a window function W that satisfies the requirements of our error theory with ε = 1.
Case 2: If µ = k +η with k ∈ N and 0 < η ≤ 1, we have W ( j) ∈ A C ([−1, 1]) for all 0 ≤ j ≤ k and and, for 0 < η < 1, we have W (k+1) ∈ L p ([−1, 1]) for all 1 ≤ p < 1 1−η with When using the window W from (17) for the low-pass filter A L in the FBP method, our error theory from Theorems 4 and 5 predicts that, for 0 ≤ σ ≤ α, the H σ -norm of the inherent FBP reconstruction error e L = f − f L behaves like In particular, for L → ∞ and the rate of convergence saturates at order O(L −µ ). For the L 2 -case, i.e., σ = 0, this behaviour will be observed numerically in the following section.

Numerical simulations
In this section, we finally present selected numerical examples to evaluate the inherent FBP reconstruction error numerically and to validate our error theory. In the following, however, we restrict ourselves to the L 2 -case and set σ to 0. First note that the approximate FBP reconstruction formula assumes the data R f (t, θ ) to be available for all (t, θ ) ∈ R × [0, π). In practice, however, only finitely many Radon samples are given. We assume that the target function f is compactly supported with and that the Radon data are given in parallel beam geometry (cf. [13,14]), i.e., the input data are of the form where d > 0 is the spacing of 2M + 1 parallel lines per angle and N is the number of angles.
The reconstruction of f from discrete Radon data (18) requires a suitable discretization of the FBP method (9), i.e., We follow a standard approach and apply the composite trapezoidal rule to discretize the convolution * and back projection B. This leads to the discrete reconstruction formula in short, Note that the evaluation of the discrete reconstruction f D requires the computation of the values (F −1 A L * D R f )(x cos(θ k ) + y sin(θ k ), θ k ) ∀ 0 ≤ k ≤ N − 1 for each reconstruction point (x, y) ∈ R 2 . To reduce the computational costs, we evaluate the function h(t, θ k ) = (F −1 A L * D R f )(t, θ k ) for t ∈ R only at the points t = t l , l ∈ Z, and interpolate the value h(t, θ k ) for t = x cos(θ k ) + y sin(θ k ) using an interpolation method I . This leads us to the discrete FBP reconstruction formula For target functions f of low regularity it is sufficient to use linear spline interpolation. To exploit a higher regularity of f we use cubic spline interpolation in our simulations. Further, according to [13, Section V.1] we couple the discretization parameters d > 0 and M, N ∈ N with the bandwidth L via and choose L to be a multiple of π, i.e., L = π · K for some K ∈ N. We remark that the different discretization steps introduce additional discretization errors that are not covered by our error theory in §3 and §4. Analysing the discretization errors is beyond the aims and scopes of this paper and for work in this direction we refer to [5,13]. In our numerical experiments, we used the popular Shepp-Logan phantom, cf. [21]. For this test case, the Radon transform can be calculated analytically, so that no errors occur during the data acquisition and the observed errors are due to the discretized approximate reconstruction method. The phantom consists of ten ellipses of constant densities, but different sizes, eccentricities and locations, see Fig. 1(a). Its attenuation function f SL is given by where each function f j , 1 ≤ j ≤ 10, is of the form of the characteristic function of an ellipse given by The parameters of the ellipses used in the Shepp-Logan phantom can be found in [21]. Its sinogram, i.e., the plot of the phantom's Radon data in the (t, θ )-plane, is shown in Fig. 1(b). As explained at the end of §2, the function f SL belongs to the Sobolev space H α 0 (R 2 ) with α < 1 2 , which is an upper bound for the decay rate of the inherent FBP reconstruction error e L = f − f L in the L 2 -norm. To observe higher rates of convergence and saturation at the rates given in §4 we need a test case of higher regularity and with analytically computable Radon transform. To this end, we consider the radially symmetric function p ν : R 2 −→ R, given by with parameter ν ∈ R >0 , which is in H α 0 (R 2 ) for any α < ν + 1 2 . Adapting the approach in [18], we now define the smooth phantom of order ν via where each function f ν j , 1 ≤ j ≤ 3, is of the form The parameters used in the definition of the smooth phantom can be found in [18]. For illustration, Fig. 2 shows f ν smooth for ν = 3 (see Fig. 2(a)) and its sinogram (see Fig. 2(b)). The FBP reconstructions of both phantoms are displayed in Fig. 3, where we used the Shepp-Logan filter (see Table 1) with L = 40π. This corresponds to M = 40 and N = 160 so that (2M + 1)N = 12960 Radon samples are taken. To measure the reconstruction error, we used the standard root mean square error, which is defined for images with I × J pixels as In our numerical experiments, we evaluated the phantoms and their FBP reconstructions with different window functions and bandwidths on a square grid with 1024 × 1024 pixels.
In our first numerical simulations we have employed four commonly used low-pass filters: Consequently, these windows satisfy the assumptions of Theorem 4 with k = 2 and ε = 1 and our error theory states that, for any function f ∈ L 1 (R 2 ) ∩ H α (R 2 ) with α > 0, the L 2 -norm of inherent FBP reconstruction error is bounded above by Fig. 4 shows the RMSE of the FBP reconstruction of the Shepp-Logan phantom f SL as a function of the bandwidth L in logarithmic scales for different window functions. In addition to the popular Shepp-Logan filter (Fig. 4(a)), we applied the Hamming filter with parameter β = 0.92 (Fig. 4(b)), the Gaussian filter with parameter β = 4.9 (Fig. 4(c)) and the Parabola filter with parameter β = 0.59 (Fig. 4(d)). These parameters were chosen such that these filters have the same value for W L ∞ ([0,1]) as the Shepp-Logan filter. Hence, the corresponding reconstruction errors should behave similarly due to our error estimate (20).
As expected, we see that the RMSE for the Shepp-Logan filter, Hamming filter with β = 0.92, Gaussian filter with β = 4.9 and Parabola filter with β = 0.59 are nearly the same. Moreover, in all four cases we observe a decrease of the RMSE with rate L −0.5 . This is exactly the behaviour we expected due to our L 2 -error estimate (20), since we have f SL ∈ H α (R 2 ) for all α < 1 2 . Fig. 5 now shows the RMSE of the FBP reconstruction for the smooth phantom f ν smooth of order ν = 3, which belongs to H α (R 2 ) for all α < 7 2 . Hence, according to the estimate (20) the convergence rate of the RMSE should saturate at order L −2 . Indeed, this behaviour can be observed in our numericals results, see Fig. 5(a)-5(d). Further, the RMSE for the Shepp-Logan filter again coincides with the RMSE of the Hamming filter with β = 0.92, Gaussian filter with β = 4.9 and Parabola filter with β = 0.59. Note that this behaviour is more pronounced for the smooth phantom f ν smooth than for the Shepp-Logan phantom f SL . In conclusion, our numerical results for C 2 -windows totally comply with our L 2 -error theory with k = 2 and p = ∞.
In our second set of numerical simulations we used the generalized polynomial filter with window of order µ ∈ R >0 and with jump height β ∈ [0, 1), see Example 1. Recall that for this filter our error theory from Theorems 4 and 5 states that, for any function f ∈ H α (R 2 ) with α > 0, the L 2 -norm of the inherent FBP reconstruction error is bounded above by where, for fixed α and µ, the constant C α,µ,β > 0 decreases with increasing β ∈ [0, 1). In particular, the rate of convergence is predicted to saturate at fractional order L −µ . The numerical results for the reconstruction of the Shepp-Logan phantom f SL are displayed in Fig. 4(e)-4(h) and can be summarized as follows. For µ = 0.2, the convergence rate of the RMSE saturates at fractional order L −0.2 . Moreover, when increasing the parameter β from β = 0 to β = 0.2, the RMSE decreases (Fig. 4(e)-4(f)). But for µ ∈ {0.9, 2.7}, the RMSE behaves like L −0.5 , see Fig. 4(g)-4(h), where we always chose the jump height β = 0.8. Since f SL ∈ H α (R 2 ) for all α < 1 2 , this is exactly the behaviour we expected due to our L 2 -error estimate (21). In particular, for µ ∈ {0.9, 2.7} the rate of convergence is given by the smoothness of the target function f SL .  In contrast to that, our numerical results for the reconstruction of the smooth phantom f ν smooth of order ν = 3 show that the rate of convergence saturates for all our choices of µ, see Fig. 5(e)-5(h). Indeed, for µ = 0.2, the convergence rate of the RMSE again saturates at fractional order L −0.2 and, further, increasing the parameter β from β = 0 to β = 0.2 decreases the RMSE (Fig. 5(e)-5(f)). Also for µ ∈ {0.9, 2.7}, the RMSE behaves like L −µ , see Fig. 5(g)-5(h). But this was expected, since we have f ν smooth ∈ H α (R 2 ) for all α < ν + 1 2 . Consequently, our numerical results again totally comply with our L 2 -error theory and especially the saturation of the convergence rate at fractional order is observable.
In our final set of numerical simulations we considered the generalized ramp filter with window for |S| > 1 of width β ∈ (0, 1) and jump height γ ∈ [0, 1]. Note that choosing γ = 1 gives the classical Ram-Lak filter (see Table 1). Furthermore, these filters satisfy Assumption (A) and our error theory from Theorem 3 (with ε = 1) states that, for any function f ∈ H α (R 2 ) with α > 0, the L 2 -norm of the inherent FBP reconstruction error e L = f − f L is bounded above by In particular, the rate of convergence is always determined by the smoothness α of the target function f . Furthermore, for fixed L and f , we see that the L 2 -error bound decreases when increasing the window's width β ∈ (0, 1) or jump height γ ∈ [0, 1]. In all our simulations for the reconstruction of the Shepp-Logan phantom f SL we observe that the RMSE behaves like L −0.5 , as predicted by the theory, see Fig. 6. Moreover, increasing the width β of the window W leads to an decrease of the RMSE (Fig. 6(a)-6(b)). This can also be seen when fixing β and increasing the jump height γ (Fig. 6(c)-6(d)). Consequently, we exactly observe the behaviour predicted by the L 2 -error estimate (22).
When considering the FBP reconstruction of the smooth phantom f ν smooth of order ν = 3, we see that for all choices of the parameters β and γ the RMSE behaves like L −3.5 , see Fig. 7. Thus, the rate of convergence is determined by the smoothness of the target function and the numerical observations comply with our L 2 -error estimate (22). Further, as for the Shepp-Logan phantom, increasing the width β results in a decrease of the RMSE (Fig. 7(a)-7(b)). The same holds true when fixing β and increasing the jump height γ (Fig. 7(c)-7(d)).
In conclusion, our numerical results totally comply with our L 2 -error theory. In particular, we observe that the decay rate of the L 2 -error is indeed determined by the smoothness α of the target function f if Assumption (A) is satisfied. Finally, we remark that for both phantoms the RMSE is minimal for the Ram-Lak filter. This also complies with our error theory for band-limited low-pass filters with supp(W ) ⊆ [−1, 1], since in the L 2 -error estimate f − f L L 2 (R 2 ) ≤ Φ

Conclusion
We have proven Sobolev error estimates and convergence rates for the method of filtered back projection to approximate bivariate functions from fractional Sobolev spaces, where we required W ∈ L ∞ (R) and | · |W (·) ∈ L 2 (R) for the window W instead of supp(W ) ⊆ [−1, 1].
For target functions f ∈ L 1 (R 2 ) ∩ H α (R 2 ) of smoothness α > 0 we proved H σ -error estimates for all 0 ≤ σ ≤ α, where we, in particular, considered the following special case. If, for ε > 0, W satisfies W ∈ C k−1 ([−ε, ε]) with W (0) = 1 and W ( j) (0) = 0 for all 1 ≤ j ≤ k −1 as well as W (k) ∈ L p ([−ε, ε]) for 1 < p ≤ ∞, the H σ -norm of the FBP reconstruction error satisfies f − f L σ = O L − min{k−1/p,α−σ } for L → ∞ and, thus, the decay rate saturates at fractional order L −(k−1/p) . This generalizes our results in [3], where we considered windows W with supp(W ) ⊆ [−1, 1] and assumed regularity of W on [−1, 1] instead of [−ε, ε] for ε > 0. Hence, our results substantiate the conclusion that the flatness of W at zero determines the decay rate of the inherent FBP reconstruction error. Finally, we remark that our theoretical error estimates provide upper bounds on the inherent FBP approximation error which is incurred by filtering the Radon data. However, the numerical simulations demonstrate that the predicted saturation of the decay rate depending on the flatness of the filter's window at 0 is indeed intrinsic and not an artefact of the proof.