Error analysis for filtered back projection reconstructions in Besov spaces

Filtered back projection (FBP) methods are the most widely used reconstruction algorithms in computerized tomography (CT). The ill-posedness of this inverse problem allows only an approximate reconstruction for given noisy data. Studying the resulting reconstruction error has been a most active field of research in the 1990s and has recently been revived in terms of optimal filter design and estimating the FBP approximation errors in general Sobolev spaces. However, the choice of Sobolev spaces is suboptimal for characterizing typical CT reconstructions. A widely used model are sums of characteristic functions, which are better modelled in terms of Besov spaces $\mathrm{B}^{\alpha,p}_q(\mathbb{R}^2)$. In particular $\mathrm{B}^{\alpha,1}_1(\mathbb{R}^2)$ with $\alpha \approx 1$ is a preferred model in image analysis for describing natural images. In case of noisy Radon data the total FBP reconstruction error $$\|f-f_L^\delta\| \le \|f-f_L\|+ \|f_L - f_L^\delta\|$$ splits into an approximation error and a data error, where $L$ serves as regularization parameter. In this paper, we study the approximation error of FBP reconstructions for target functions $f \in \mathrm{L}^1(\mathbb{R}^2) \cap \mathrm{B}^{\alpha,p}_q(\mathbb{R}^2)$ with positive $\alpha \not\in \mathbb{N}$ and $1 \leq p,q \leq \infty$. We prove that the $\mathrm{L}^p$-norm of the inherent FBP approximation error $f-f_L$ can be bounded above by \begin{equation*} \|f - f_L\|_{\mathrm{L}^p(\mathbb{R}^2)} \leq c_{\alpha,q,W} \, L^{-\alpha} \, |f|_{\mathrm{B}^{\alpha,p}_q(\mathbb{R}^2)} \end{equation*} under suitable assumptions on the utilized low-pass filter's window function $W$. This then extends by classical methods to estimates for the total reconstruction error.


Introduction
We consider the classical inverse problem of reconstructing a function f : Ω → R, Ω ⊂ R 2 , from its line integrals, which is the mathematical model underlying X-ray computerized tomography (CT). The line integrals of f are defined by the Radon transform Rf ≡ Rf (t, θ), given by where the set t,θ = {(x, y) | x cos(θ) + y sin(θ) = t} ⊂ R 2 denotes the unique straight line that is orthogonal to the unit vector n θ = (cos(θ), sin(θ)) and has (signed) distance t to the origin, i.e., t,θ passes through (t cos(θ), t sin(θ)) ∈ R 2 .
The assumption of a continuous f is necessary in order to ensure that the inversion formula holds pointwise, cf. [15,3]. However, the continuity assumption is not needed for error estimates concerning regularized reconstruction algorithms as discussed in the sequel of this paper. The analytical inversion is unstable with respect to highly oscillating variations of g = Rf , which motivates the introduction of a regularized inversion formula, the socalled method of filtered back projection (FBP). FBP is based on introducing a window function with W ∞ < ∞, which has either bounded support or decays fast enough at infinity. The window function W is scaled by a parameter L and we introduce the low-pass filter A L (S) = |S| W ( S /L) for S ∈ R to replace the Fourier multiplier |S| in the Riesz potential leading to the approximate Here, L serves as regularization parameter and will be adapted depending on the noise level in the data. Window functions of typical low-pass filters are displayed in Table 1.
In case of noisy data g δ with g δ − Rf ≤ δ we compute the reconstruction f δ L (x, y) = 1 2 B(F −1 [A L (S)F(g δ )(S, θ)])(x, y) and the total reconstruction error splits into an approximation error whose analysis is the main target of the present paper, and a data error f L − f δ L . Such FBP methods are the most widely used reconstruction algorithms in computerized tomography (CT). Studying the resulting reconstruction error has been a most active field of research in the 1990s and has recently been revived in terms of optimal filter design and estimating the approximation errors in general Sobolev spaces [2,3,4]. We will review the state of research in more detail in the next section.
For motivation of the present paper, we note that the choice of Sobolev spaces is suboptimal for characterizing typical CT reconstructions. A widely used model for CT reconstructions are sums of characteristic functions, which are better modelled in terms of Besov spaces B α,p q (R 2 ). In particular, B α,1 1 (R 2 ) with α close to 1 is a preferred model in image analysis for describing natural images [6,17]. Therefore, in this paper we focus on extending the analysis of the approximation error of FBP reconstructions to target functions f ∈ L 1 (R 2 ) ∩ B α,p q (R 2 ) with positive α ∈ N and 1 ≤ p, q ≤ ∞. We prove that the L p -norm of the inherent FBP approximation error f − f L can be bounded above by under suitable assumptions on the utilized low-pass filter's window function W . This then extends by classical methods to estimates for the total reconstruction error.
The transition from error estimates in Sobolev spaces to Besov spaces requires substantially different techniques. The Sobolev space estimates in [2,3,4] are implicitly based on Plancherel's formula, which states that the Fourier transform is an isometry between L 2 -spaces, which is not the case for L p -spaces with p = 2. Hence, the definition of Besov spaces in terms of moduli of smoothness requires different analytical tools for estimating the reconstruction error with respect to the target function's Besov norm. The paper is organized as follows. In Section 2 we review the state of the art concerning approximation errors of FBP reconstruction for functions on unbounded domains. We then introduce the definition of Besov spaces used in this paper along with some technical Lemmata, which will be needed later for estimating f − f L L p (R 2 ) . Section 4 then contains the main results of the paper. The proofs are split into the cases 0 < α < 1, 1 ≤ p, q ≤ ∞ and 1 < α, α / ∈ N, 1 ≤ p, q ≤ ∞. We also include a straight forward result on how these results on e L = f − f L extend to an estimate of the total approximation error for noisy data.
Section 5 contains some numerical experiments confirming the theoretical findings of the previous section.

State of the art
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. We shortly summarize the available literature on estimating total reconstruction and approximation errors for FBP reconstructions. For a general reference we refer to the standard textbooks [15,16] and to the introductory chapters of [3], which contains an in-depth description and comparison of the results by Madych, which are most relevant for our approach. Indeed, the description of the state of the art in [3] serves as the main reference for the following summary.
FBP algorithms and their approximation properties were explicitely or at least implicitely addressed already in the very first papers and textbooks [7,15] on the mathematics of computerized tomography. Arguably the first paper addressing an analysis of e L = f − f L in a classical function space setting is [18]. There, Popov showed pointwise convergence however with a restriction to a small class of piecewise smooth functions. Pointwise convergence and L ∞ -error estimates for e L are also discussed by Munshi et al. in [13,12]. Their results are supported by numerical experiments in [14].
The approach of Rieder and Schuster [23] leads to L 2 -convergence with suboptimal rates for compactly supported Sobolev functions. In contrast to this, in [20,21] 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 [19] showed convergence without rates in the L 2 -norm for compactly supported L ∞ -functions and in points of continuity under additional assumptions. Note that [19] deals with the continuous problem, while [20,23,21] discuss discrete settings.
More relevant for our present paper is the approach described by Madych, who proves error bounds on the L p -norm of e L in terms of L p -moduli of continuity of the target function f , see [11]. Madych chooses a convolution kernel K : R 2 → R as an approximation of the identity and computes the convolution product f * K L to approximate the target function f , where, for L > 0, the scaled kernel K L is given by If K is chosen to be a uniform sum of ridge functions, the convolution f * K L can be expressed in terms of the Radon data Rf as in the approximate FBP formula (1), see [11,Proposition 1]. The assumptions on K for proving these results are rather restrictive, in particular they require continuous filter functions, which e.g. excludes the case of a ramp filter. Using an essentially different approach, [2,3,4] then proved Sobolev space estimates in a more general setting with substantially weaker assumptions on the filter, including all classical choices.
To some extend, the approach of Madych is a special case of the mollifier approaches used in [9,10,22]. However, neither Sobolev nor Besov space error estimates for the continuous case are derived in these papers.

Besov spaces and technical Lemmata
The focus of the paper is on analysing the L p -norm of the inherent FBP reconstruction for α = n + θ with n ∈ N and 0 < θ < 1. As a general reference for properties of Besov spaces we refer to [8].
We start by proving some technical Lemmata, which will be needed in the subsequent sections for estimating e L . The critical part in these estimates is to control the L p -, resp. L ∞ -norm of the modulus of smoothness in the definition of the Besov semi-norm.
We start with an L p -estimate.

6
The proof of this Lemma is mostly technical and has been moved to Appendix A. We now use the previous Lemma for proving an L ∞ -estimate, which is equivalent to an embedding into the L p -setting.
Lemma 2 Let 1 ≤ q < ∞ and α > 0. Further, let g : [0, ∞) → [0, ∞) be an increasing function. Then, for any c > 1, The proof of this lemma can be found in Appendix B. The classical estimates are only concerned with fixed α and qualitative estimates of the constant involved in the Besovnorm estimates. However, the asymptotic behaviour for α 0 is needed for a refined analysis in the next section.
Lemma 3 Let 1 ≤ q < ∞ and let g : [0, ∞) → [0, ∞) be increasing and bounded from above. Further, assume that there exists some σ ∈ (0, 1) such that Then, Again, the proof of this lemma has been moved to the appendix, see Appendix C.

Approximation error in Besov spaces
We now turn to estimating the approximation error e L = f − f L of FBP reconstructions in L p -norms under the assumption that f ∈ L 1 (R 2 ) ∩ B α,p q (R 2 ). We will discuss the cases 0 < α < 1 and α > 1 separately. These estimates are then used for deriving a bound on the total FBP reconstruction error f − f δ L for noisy data g δ . As already stated, we consider the approximate filtered back projection (FBP) formula with a given low-pass filter of finite bandwidth L > 0 and with even window function W ∈ L ∞ (R) satisfying Recall that for target functions f ∈ L 1 (R 2 ) the approximate FBP reconstruction f L ∈ L 1 loc (R 2 ) is defined almost everywhere on R 2 and can be rewritten as where the convolution kernel K L ∈ L 2 (R 2 ) is given by

Error Estimate for
Several papers, see e.g. [6,17], have argued, that natural images including cross sections of the human body can be modelled by Besov spaces with α < 1. This also includes the case of functions which are superpositions of characteristic functions of smooth domains, which serves as a standard model for simulation in tomography. Hence, we start with analyzing the case 1 ≤ p, q ≤ ∞ and 0 < α < 1.
Then, the L p -norm of the inherent FBP reconstruction error e L = f − f L is bounded above by Proof First note that due to f ∈ B α,p q (R 2 ) ⊂ L p (R 2 ) and K ∈ L 1 (R 2 ), we have K L ∈ L 1 (R 2 ) and Furthermore, K L and W are related via Thus, for (x, y) ∈ R 2 holds that Thus, with the L ∞ -modulus of continuity Thus, with the L p -modulus of continuity where we use the scaling property Recall further that the convolution kernel K is radially symmetric, i.e., there exists a univariate function k : R → R such that Thus, transforming to polar coordinates gives For q = ∞ we can conclude that Now, let 1 ≤ q < ∞. Since the L p -modulus of continuity is monotonically increasing in δ > 0, we can apply Lemma 2 to obtain It remains to optimize the constant For c > 1, we have as well as Consequently, the unique minimizer of C α,q on R >1 is given by c * = exp ((2αq) −1 ) and which completes the proof.
Note that for fixed 1 ≤ q < ∞ the constant c α,q in Theorem 4 goes to 0 for α 0 as α 1 /q . However, an application of Lemma 3 to the L p -modulus of continuity ω p (f, ·) shows that in this case the Besov semi-norm |f | B α,p q (R 2 ) goes to ∞ for α 0 as α − 1 /q and, in particular, we have Thus, for α 0, the L p -error estimate in Theorem 4 reduces to and, for q → ∞, we obtain which is consistent with simply applying Young's inequality, as for f ∈ L p (R 2 ) we have

Error Estimate for α > 1
Convergence results in the general regularization theory for inverse problems typically depend on additional smoothness assumptions on f . Hence, we now consider the case α > 1, i.e., functions f which are slightly smoother than sums of characteristic functions.
Then, the L p -norm of the inherent FBP reconstruction error e L = f − f L is bounded above by Proof To start with, we recall that due to f ∈ B α,p q (R 2 ) ⊂ L p (R 2 ) and K ∈ L 1 (R 2 ) we have Hence, for (x, y) ∈ R 2 follows that To continue, we first assume that f ∈ B α,p q (R 2 ) ∩ C ∞ (R 2 ). Then, Taylor's theorem gives where we set f for the sake of brevity. By using the assumed moment conditions on K and its scaling property Consequently, we obtain Thus, for the L p -norm of the inherent FBP reconstruction error e L = f − f L follows that where, for 1 ≤ p < ∞, we applied Minkowski's integral inequality as in the proof of Theorem 4.
Recall that the convolution kernel K is radially symmetric, i.e., there exists a univariate function k : R → R such that Hence, transforming to polar coordinates gives For q = ∞ we can conclude that and, thus, For 1 ≤ q < ∞ we can apply Lemma 2 and obtain, as in the proof of Theorem 4, To prove the result also for functions f ∈ B α,p q (R 2 ) that are not smooth, let ϕ ∈ C ∞ c (R 2 ) be a standard mollifier function, i.e., let ϕ ≥ 0 satisfy supp(ϕ) ⊆ B 1 (0) and Moreover, for ε > 0, we define Furthermore, we have already proven that For all L > 0 we now have Thus, taking the limit ε 0 gives for any f ∈ L 1 (R 2 ) ∩ B α,p q (R 2 ) and the proof is complete.

Error estimates for noisy data
We now consider the case of noisy Radon data. To this end, we assume that the Radon data Rf ∈ L p (R × [0, π)), 1 ≤ p ≤ ∞, is known only up to a noise level δ > 0 so that we have to reconstruct the target function f from given noisy measurements g δ ∈ L p (R × [0, π)) satisfying Rf − g δ L p (R×[0,π)) ≤ δ. By applying the approximate FBP formula to the noisy data g δ , we obtain the reconstruction and the overall FBP reconstruction error e δ L = f − f δ L can be split into an approximation error term and a data error term, In the following, we assume that f is supported in a compact set Ω ⊂ R 2 and analyse the L p -norm of the overall FBP reconstruction error e δ L on Ω with respect to the noise level δ. By the triangle inequality, we have and, consequently, we can treat the approximation error and the data error separately.
The analysis of the data error f L − f δ L is based on the fact that the back projection B defines a mapping The case 1 ≤ p < ∞ is discussed in the following lemma.
We are now prepared to analyse the data error f L − f δ L in the L p -norm for target functions f ∈ L p (Ω) satisfying Rf ∈ L p (R × [0, π)), where with noisy measurements g δ ∈ L p (R × [0, π)).

Theorem 7 (Data error)
For compact domain Ω ⊂ R 2 and 1 ≤ p ≤ ∞ let f ∈ L p (Ω) satisfy Rf ∈ L p (R × [0, π)). Furthermore, let W ∈ L ∞ (R) be even such that the corresponding filter A ≡ A 1 satisfies A ∈ L 1 (R) ∩ L 2 (R) as well as F −1 A ∈ L 1 (R). Finally, for δ > 0, let g δ ∈ L p (R × [0, π)) be given with Rf − g δ L p (R×[0,π)) ≤ δ. Then, the L p -norm of the data error f L − f δ L on Ω is bounded above by Proof We first consider the case 1 ≤ p < ∞. By the linearity of the back projection B and Lemma 6 we obtain where we set e δ = Rf − g δ . By definition of the low-pass filter A L , for t ∈ R we have Consequently, an application of Young's inequality gives By combining the estimates we can conclude that Now, let p = ∞. Following along the lines from before, we obtain and the proof is complete.
By combining the above result for the data error with our previous findings for the approximation error we can now estimate the L p -norm of the overall FBP reconstruction . Corollary 8 (Convergence rates for noisy data) Let f ∈ L 1 (R 2 ) ∩ B α,p q (R 2 ) for α = n + θ with n ∈ N, 0 < θ < 1 and 1 ≤ p, q ≤ ∞ be supported in a compact domain Ω ⊂ R 2 . Furthermore, let W ∈ L ∞ (R) be even with W (0) = 1 such that the corresponding filter A ≡ A 1 satisfies A ∈ L 1 (R) ∩ L 2 (R) as well as F −1 A ∈ L 1 (R) and the convolution kernel K ≡ K 1 satisfies K ∈ L 1 (R 2 ) and R 2 x j 1 y j 2 K(x, y) d(x, y) = 0 ∀ 1 ≤ |j| ≤ n as well as Finally, let g δ ∈ L p (R × [0, π)) be given with Rf − g δ L p (R×[0,π)) ≤ δ. Then, the L p -norm of the overall FBP reconstruction error e δ L = f − f δ L on Ω is bounded above by where the bandwidth L is chosen as and the involved constants are given by and c W,α,q = c θ,q Γ(θ + 1) Γ(α + 1) In particular, Proof According to Theorem 5 the L p -norm of the approximation error f − f L is bounded above by and by Theorem 7 the L p -norm of the data error f L − f δ L can be estimated in terms of the noise level δ via Thus, the L p -norm of the overall FBP reconstruction error f − f δ L can be bounded above by ,Ω,p L δ and choosing the bandwidth L as as stated.
Note that the decay rate of the L p -error bound in Corollary 8 is independent of 1 ≤ p ≤ ∞, as where the filter's bandwidth L > 0 has to go to ∞ as the noise level δ > 0 goes to 0 with rate To close this section we give an example of a filter function A L,ν depending on a parameter ν ∈ N 0 , which fulfils the assumptions of Theorems 4, 5 and 7 for suitable ν. We remark that these assumptions especially imply continuity of the window W on R so that they are not satisfied for the typical choices summarized in Table 1.

Example 9
We define the smooth filter of order ν ∈ N 0 as A L,ν = | · | W ν ( · /L) with It is easy to verify that W ν ∈ L ∞ (R) is an even function with W (0) = 1 and A L,ν ∈ L 1 (R) ∩ L 2 (R) for all ν ∈ N 0 and L > 0. In the following, we analyse the inverse Fourier transform of A L,ν and the corresponding convolution kernel K L,ν . The inverse Fourier transform of A L,ν involves the generalized hypergeometric function 1 F 2 and is given by It can be verified that F −1 A L,ν ∈ L 1 (R) for all ν ∈ N, but F −1 A L,ν ∈ L 1 (R) for ν = 0. The corresponding convolution kernel K L,ν can be computed as for (x, y) R 2 = 0, where J ν+1 denotes the Bessel function of the first kind of order ν + 1 defined by For fixed α ≥ 0 we now develop a condition on ν ∈ N 0 such that the integral is finite, where we set K ν ≡ K 1,ν for the sake of brevity. To this end, first observe that Since |J ν+1 (r)| r α−ν is bounded on [0, η] for all η > 0 due to [1, 9.1.7 & 9.1.60], we have According to [1, 9.2.1], the Bessel function J ν+1 may be written as Therefore, choosing η sufficiently large yields | cos(r − 1 /2 π ν − 3 /4 π)| r ν−α+ 1 /2 dr and the latter integral converges for ν > α + 1 /2. We now prove divergence of I α,ν for ν ≤ α + 1 /2. For sufficiently large N ∈ N we have In particular, we have proven that K ν ∈ L 1 (R 2 ) if and only if ν > 1 /2. It remains to determine the maximal n ∈ N such that By transforming to polar coordinates the above integral can be rewritten as r ν−|j| dr = 0 ∀ 1 ≤ |j| ≤ n and, moreover, 2π 0 cos j 1 (ϕ) sin j 2 (ϕ) dϕ = 0 ∀ 1 ≤ |j| ≤ n ⇐⇒ n = 1.
Consequently, the smooth filter of order ν ∈ N satisfies the moment conditions only for n = 1.
Summarizing the results of Example 9, the smooth filter of order ν ∈ N 0 satisfies the assumptions of our error theory in Theorems 4, 5 and 7 for all α < 2 iff ν > α + 1 /2.

Numerical experiments
We now present selected numerical examples to illustrate our theoretical results. To this end, we assume that the target function f is compactly supported with and that the Radon data are given in parallel beam geometry where d 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 requires a suitable discretization of the approximate FBP reconstruction formula We follow a standard approach [16] and apply the composite trapezoidal rule to discretize the convolution * and back projection B, leading to the discrete convolution * D and discrete back projection B D , respectively. Moreover, we apply an interpolation method I to reduce the computational costs. This yields 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 we apply cubic spline interpolation. Furthermore, we couple the 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.
In our numerical experiments, we use the Shepp-Logan phantom with attenuation function where each function f j 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 [24]. For illustration, the Shepp-Logan phantom and its sinogram are shown in Figure 1. According to [5], the function f SL belongs to the Besov space B α,p p (R 2 ) for α < 1 /p, which determines the decay rate of the FBP approximation error e L = f − f L in the  L p -norm according to Theorem 4. To observe higher rates of convergence we consider the function with parameter σ > 0, which is in B α,p p (R 2 ) for α < σ + 1 /p. Adapting the approach in [20], we then define the smooth phantom of order σ via where each function f σ j is of the form f σ (x, y) = p σ (x a,b,h,k,ϕ (x, y)).
The parameters used in the definition of the smooth phantom can be found in [20]. For illustration, Figure 2 shows the smooth phantom of order σ = 1 along with its sinogram.    The FBP reconstructions of both phantoms are displayed in Figure 3, where we use the smooth filter from Example 9, i.e., with ν = 5 and L = 100π. This corresponds to M = 100 and N = 315 so that (2M + 1)N = 63315 Radon samples are taken. In our numerical experiments, we evaluate the phantoms and reconstructions on a square grid with 1024 × 1024 pixels. We start with illustrating our theoretical results concerning the approximation error To this end, Figure 4 shows the discrete L p -norm of the FBP approximation error of the Shepp-Logan phantom for p ∈ {1, 4 /3, 2, 4} as a function of the bandwidth L in logarithmic scales for the smooth filter with ν ∈ {5, 7}. In this case, Theorem 4 gives In all cases, the plots in Figure 4 show that the discrete L p -norm of the FBP approximation error decreases with increasing bandwidth L with rate L − 1 /p for both ν = 5 and ν = 7. This is exactly the behaviour we expect as we have f SL ∈ B α,p p (R 2 ) for any α < 1 /p. Moreover, we see that the error is smaller for ν = 5 than for ν = 7. This observation also fits to our expectations because the constant is smaller for ν = 5 for all corresponding values of α ∈ {1, 3 /4, 1 /2, 1 /4}, see Table 2. Note that this behaviour of the error can also be observed for other choices of 1 ≤ p < ∞ and the smooth filter with alternative parameter ν ∈ N.   logarithmic scales for the smooth filter with ν ∈ {5, 7}. In this case, Theorem 5 gives where θ = α − α is the fractional part of α > 1. In all cases, the plots in Figure 5 show that the discrete L p -norm of the approximation error decreases as L −(1+ 1 /p) for both ν = 5 and ν = 7. This is exactly the behaviour we expect, as we have f 1 smooth ∈ B α,p p (R 2 ) for any α < 1 + 1 /p. Moreover, we see that the error is smaller for ν = 5 than for ν = 7. This again fits to our expectations as the constant c α,K is smaller for ν = 5 for α ∈ {2, 7 /4, 3 /2, 5 /4}, see Table 2.
Recall that the smooth filter of order ν ≥ 3 satisfies the moment conditions R 2 x j 1 y j 2 K(x, y) d(x, y) = 0 ∀ 1 ≤ |j| ≤ n only for n = 1 so that Theorem 5 can only predict a decrease of the error with rate L −2 at most. Thus, our theory predicts saturation of the error decay rate if the smoothness α of the target function is larger than 2. To illustrate this saturation numerically, we use the smooth phantom of order σ = 2, which satisfies f 2 smooth ∈ B α,p p (R 2 ) for any α < 2 + 1 /p. Figure 6 shows the discrete L p -norm of the corresponding FBP approximation error exemplarily for p ∈ {1, 4}, where we use the smooth filter with ν = 5. We indeed observe that the error decreases as L −2 . This is true for all 1 ≤ p < ∞ and ν ≥ 3. Hence, the predicted saturation of the decay rate is observable in numerical experiments.  In summary, our reported numerical results totally comply with our theoretical findings concerning the FBP approximation error.
In our second set of numerical experiments we investigate the FBP data error denotes the approximate FBP reconstruction from noisy Radon measurements g δ with noise level δ > 0. To this end, we assume that we are given noisy measurements More precise, in our numerical simulations we use additive white Gaussian noise with noise-level is the arithmetic mean of the absolute values of the Radon samples of f . Moreover, we again use the smooth filter with parameter ν ∈ {5, 7}. Recall that, according to Theorem 7, the L p -norm of the data error can be estimated as where we have F −1 A L 1 (R 2 ) = 0.2976 for ν = 5 and F −1 A L 1 (R 2 ) = 0.2541 for ν = 7. Consequently, we expect the error to be smaller for ν = 7 and that the error increases with increasing L with rate L 1 , which is independent of the integrability parameter p. Figure 7 shows the discrete L p -norm of the FBP data error for the Shepp-Logan phantom for p ∈ {1, 4 /3, 2, 4} as a function of the bandwidth L in logarithmic scales. The results for the smooth phantom of order σ = 1 are summarized in Figure 8. In all cases, we observe that the data error increases with L with rate L 1 /2 , which is indeed independent of p. However, the growth rate is overestimated in Theorem 7, where the error bound increases with rate L 1 . Moreover, our numerical experiments show that the data error is indeed smaller for ν = 7 than for ν = 5, as suggested by our error estimate.
We wish to remark that the reported error behaviour can also be observed in numerical experiments with other choices of p, ν and δ. Moreover, we note that the correct growth rate L 1 /2 of the FBP data error was derived in [3] for the case p = 2.

Conclusion
In this paper, we have analysed the approximation and the total reconstruction error of the FBP method for CT reconstructions from parallel beam data. Our results depend on the smoothness and the flatness of the filter's window function near 0. Moreover, the rate of convergence depends on the smoothness of the true solution f ∈ B α,p q (R 2 ). The convergence results cover the cases 1 ≤ p, q ≤ ∞ and 0 < α < 1, resp. 1 < α < ∞, α / ∈ N. Integer values for α require slightly different definitions of the Besov spaces. This case is not covered by our analysis.
The numerical examples with phantoms of different Besov smoothness confirm the theoretical convergence rates of the approximation error, which require to link the flatness of the window function at 0 with the smoothness of the target function for optimal convergence rates. The data error, however, is overestimated by our theory and requires further in-depth research.