ROI reconstruction from truncated cone-beam projections

Region-of-Interest (ROI) tomography aims at reconstructing a region of interest \begin{document} $C$ \end{document} inside a body using only x-ray projections intersecting \begin{document} $C$ \end{document} and it is useful to reduce overall radiation exposure when only a small specific region of a body needs to be examined. We consider x-ray acquisition from sources located on a smooth curve \begin{document} $Γ$ \end{document} in \begin{document} $\mathbb R^3$ \end{document} verifying the classical Tuy condition. In this generic situation, the non-trucated cone-beam transform of smooth density functions \begin{document} $f$ \end{document} admits an explicit inverse \begin{document} $Z$ \end{document} as originally shown by Grangeat. However \begin{document} $Z$ \end{document} cannot directly reconstruct \begin{document} $f$ \end{document} from ROI-truncated projections. To deal with the ROI tomography problem, we introduce a novel reconstruction approach. For densities \begin{document} $f$ \end{document} in \begin{document} $L^{∞}(B)$ \end{document} where \begin{document} $B$ \end{document} is a bounded ball in \begin{document} $\mathbb R^3$ \end{document} , our method iterates an operator \begin{document} $U$ \end{document} combining ROI-truncated projections, inversion by the operator \begin{document} $Z$ \end{document} and appropriate regularization operators. Assuming only knowledge of projections corresponding to a spherical ROI \begin{document} $C \subset B$ \end{document} , given \begin{document} $e >0$ \end{document} , we prove that if \begin{document} $C$ \end{document} is sufficiently large our iterative reconstruction algorithm converges at exponential speed to an \begin{document} $e$ \end{document} -accurate approximation of \begin{document} $f$ \end{document} in \begin{document} $L^{∞}$ \end{document} . The accuracy depends on the regularity of \begin{document} $f$ \end{document} quantified by its Sobolev norm in \begin{document} $W^5(B)$ \end{document} . Our result guarantees the existence of a critical ROI radius ensuring the convergence of our ROI reconstruction algorithm to an \begin{document} $e$ \end{document} -accurate approximation of \begin{document} $f$ \end{document} . We have numerically verified these theoretical results using simulated acquisition of ROI-truncated cone-beam projection data for multiple acquisition geometries. Numerical experiments indicate that the critical ROI radius is fairly small with respect to the support region \begin{document} $B$ \end{document} .


(Communicated by Gabriele Steidl)
Abstract. Region-of-Interest (ROI) tomography aims at reconstructing a region of interest C inside a body using only x-ray projections intersecting C and it is useful to reduce overall radiation exposure when only a small specific region of a body needs to be examined. We consider x-ray acquisition from sources located on a smooth curve Γ in R 3 verifying the classical Tuy condition. In this generic situation, the non-trucated cone-beam transform of smooth density functions f admits an explicit inverse Z as originally shown by Grangeat. However Z cannot directly reconstruct f from ROI-truncated projections. To deal with the ROI tomography problem, we introduce a novel reconstruction approach. For densities f in L ∞ (B) where B is a bounded ball in R 3 , our method iterates an operator U combining ROI-truncated projections, inversion by the operator Z and appropriate regularization operators. Assuming only knowledge of projections corresponding to a spherical ROI C ⊂ B, given > 0, we prove that if C is sufficiently large our iterative reconstruction algorithm converges at exponential speed to an -accurate approximation of f in L ∞ . The accuracy depends on the regularity of f quantified by its Sobolev norm in W 5 (B). Our result guarantees the existence of a critical ROI radius ensuring the convergence of our ROI reconstruction algorithm to an -accurate approximation of f . We have numerically verified these theoretical results using simulated acquisition of ROI-truncated cone-beam projection data for multiple acquisition geometries. Numerical experiments indicate that the critical ROI radius is fairly small with respect to the support region B.
in CT, several strategies have been explored such as sparsifying the numbers of xray projections or truncating the projections so that only x-rays intersecting a small Region-of-Interest (ROI) are acquired. Reconstructing a fully generic density f from its ROI-truncated projections is an ill-posed problem, and small perturbations of the projections may lead to significant reconstruction errors. For non-truncated projection data, many regularized reconstruction formulas have been introduced, such as the classical Filtered Back-Projection or the FDK algorithms [24,Ch.5]. When projections are truncated by eliminating all rays which do not meet a spherical ROI, the reconstruction problem may become severely ill posed and non-uniquely solvable [23]. For the interior problem where projection data are only known for rays intersecting a region strictly inside the support of the density f , there is no unique solution in general as shown in 2D by [23]. As a result, naive numerical reconstruction algorithms, such as direct applications of classical reconstruction formulas with all missing projection data set to zero, typically produce serious instability and unacceptable visual artifacts.
The problem of ROI recontruction in CT has been studied in multiple papers and using a variety of methods (see review [28]). Remarkable theoretical results have been obtained, for instance, by reducing the ROI interior problem to solving a family of one-dimensional inversions of truncated Hilbert transforms [1,5,32]. It was shown that unique and stable inversion from ROI truncated data is guaranteed when the density f is already known on some adequate subset inside the ROI [19,25,3,31]. Another group of more recent results show that unique and stable inversion from ROI truncated data can also be achieved without any precise subregion knowledge provided that the density f is an unknown polynomial of known degree [8,16,18,30]. To obtain practical numerical reconstructions, most such ROI CT methods rely on iterative algorithms combined with a regularization term which is frequently based on total variation minimization (e.g., [8,30]) or wavelet techniques (e.g., [18]). Even though iterative methods may be computationally intensive, especially for 3D data, advances in computational capabilities (e.g., [29]) and recent ideas from compressed sensing [34] offer powerful tools to overcome this limitation.
Here, we consider the ROI reconstruction problem where an unknown density f has to be reconstructed inside a spherical ROI C using only those projections acquired by rays intersecting C. For a fixed bounded ball B ⊂ R 3 and a smooth curve Γ ⊂ R 3 of ray sources exterior to B, we consider only density functions f ∈ L ∞ (B) with compact support included in B. For any spherical ROI C ⊂ B, we denote by R C the set of all half-lines (or 'half-rays') r intersecting C and emanating from arbitrary point sources located on Γ. The C-truncated cone beam projection operator D C maps any density function f ∈ L ∞ (B) into a function D C f defined, for each half-ray r ∈ R C , by integrating f over r.
If Γ and B verify the well known geometric Tuy condition (see Definition 2.1) then Grangeat's classical formula (as recalled below) provides an exact inverse Z of the non truncated cone beam projection D, verifying ZDf = f for every C 2density f with compact support included in B. For very specific types of curves Γ, the non truncated cone beam operator D can be inverted by curve-specific explicit operators Z defined on smooth densities. This is the case for sources on a circular helix Γ and axially truncated data where Z can be implemented by the Katsevich's formula [12,14].
For arbitrary source curves Γ and generic non smooth densities f , non uniqueness cannot be a priori avoided when inverting ROI-truncated cone beam data. Hence, for any given > 0, we have deliberately aimed to construct approximate inverses which are only -accurate in L ∞ norm. These -accurate inverses are naturally never unique but they can be rigorously constructed whenever the intersection of the ROI and the support f is not too small within the support of f .
In this paper we define and study iterative reconstruction algorithms which converge at exponential speed to -accurate inverses for cone beam data truncated by arbitrary spherical ROIs strictly included in a fixed sphere B. This has naturally led us to require that the source curve Γ in R 3 and the whole sphere B verify the Tuy condition, instead of requiring only that Γ and a specific ROI verify the Tuy condition. This slightly stronger hypothesis enables to derive uniform speeds of convergence for large classes of ROIs included in B. Our reconstruction algorithm from ROI-truncated cone-beam data is based on an explicit iterative scheme, which is easily implemented provided one knows how to implement a stable inverse Z for non truncated cone beam data. In R 3 , for any generic smooth source curve Γ and any fixed sphere B verifying the Tuy condition, Z can for instance given by the well-known Grangeat formula which inverts the non truncated cone beam operator when the unknown density f is in C 2 .
To prove convergence of our algorithm in those generic cases, we first had to prove strong stability in a Sobolev norm for the Grangeat operator Z. This condition led us to assume that f has compact support and belongs to the Sobolev space W 5 . More precisely, we consider densities f in W 5 (B) with support in a fixed ball B ⊂ R 3 and generic smooth sources curves Γ surrounding B and verifying Tuy condition. Starting from the C-truncated cone beam projection data D C f , we develop an explicit iterative ROI-reconstructions converging in L ∞ (B) to an -accurate inverse Z C f , provided the radius rad(C) is larger than a critical radius ρ < rad(B). Note that the critical radius depends on . Starting with an explicit operator Z inverting the non-truncated cone beam projection D for smooth densities, such as the generic Grangeat inversion operator, we construct a smoothing regularization operator τ such that U = Zτ (D − D C ) becomes a contraction on L ∞ (B). Then we set f 0 = Zτ D C f to define iteratively the approximating sequence of densities (f j ) by (1) f Our mathematical analysis proves that, as j → ∞, the sequence (f j ) converges to an -accurate reconstructionf = Z C f of f at exponential speed in L ∞ (C), provided rad(C) is larger than the critical radius ρ < rad(B). In particular we then have These results also extend to the academic case where sources are located over a whole sphere in R 3 containing B. Note that, for discretized ROI-truncated data, our ROI-reconstruction approach can be applied whenever the non-truncated cone beam projection operator D can be inverted by an implementable formula or a "blackbox algorithm" Z that is at least well defined on smooth densities. Convergence to an -accurate inverse then has to be checked numerically, since our proofs of convergence only apply to the ideal case of continuous data and smooth densities. Unlike other existing methods, our approach does not require any explicit restriction on the ROI location or any prior knowledge of the density inside the ROI, as long as a stable inverse Z of the nontruncated cone beam projection D exists. The existence of five derivatives in L 2 theoretically required for the density f here is more general than other restrictive assumptions found in the literature, such as that the density has to be polynomial. Note however that, in practical implementations of our algorithm, one does not need to evaluate the W 5 Sobolev norm of f . In fact, by construction, for non smooth densities, whenever our ROI reconstruction algorithm is observed to converge, the (necessarily smooth) limit will simply be close to a smoothed version of the unknown f . A key condition for convergence to an -accurate inverse Z C of the C truncated cone beam projector D C is that C must have radius rad(C) larger than a critical radius ρ( ) < rad(B). Our proof yields an explicit but certainly too large theoretical upper bound for ρ( ). However in all our numerical experiments on simulated ROI truncated data (see Sec. 6 below) our numerical evaluations of ρ( ) for = 0.10 showed that ρ( ) < rad(B)/2. This indicates that in practical ROI-truncated acquisitions, our reconstruction algorithm should efficiently handle truncation by fairly large spherical ROIs, which is a quite favorable situation for exposure reduction by ROI-truncation.
Our iterative scheme (1) is formally similar to other (mostly heuristic) iterative algorithms also proposed in the literature for ROI reconstruction such as the so-called Iteration Reconstruction-Reprojection (IRR) algorithm [22,17,37], the Ordered Subsets Convex algorithm [11] proposed to speed up CT reconstruction by reducing the number of projections and the iterative maximum likelihood (ML) algorithm proposed by Ziegler et al. [37]. We provide here both a rigorous mathematical proof of convergence and a very concrete iterative construction of -accurate inverses for ROI truncated cone beam transform with sources located on a generic 3D smooth curve satisfying the classical Tuy condition. To validate our approach in the discrete setting, we performed numerical experiments using four classical discrete acquisition geometries, with sources located on a sphere, a spiral curve, a circular curve and twin orthogonal circular curves. For each setting, we simulated ROI-truncated cone beam data acquisition using three different density functions in R 3 : a Shepp-Logan phantom, a mouse tissue density data sample, a human jaw density data sample. We performed extensive numerical tests using spherical ROIs with various centers and radii and found that the numerically computed critical ROI radius is relatively small as compared to the size of the support of f and essentially insensitive to the ROI location.
Paper outline. In Sec. 2, we recall the definitions of the ray and cone-beam transforms, and the Tuy condition for x-ray acquisition with sources on smooth 3D curves. In Sec. 3, we examine known inverse operators Z implementing the reconstruction of densities from non-truncated projection data and study the continuity properties of Z on Sobolev spaces defined on the space of rays R B . In Sec. 4, we define a class of smoothing approximations of the identity in the image and projection domains, and we indicate how to implement these regularization operators by 'small' mollification. In Sec. 5, we describe our iterative ROI reconstruction algorithm from ROI-truncated data and prove our main convergence results. In Sec. 6, we present numerical implementations of our iterative ROI reconstruction for discrete acquisition setups where sources are located on (1) a sphere, (2) a spiral curve, (3) a circular curve, (4) twin orthogonal circular curves, with simulated ROI-truncated x-ray data acquired from three density functions in R 3 : a Shepp-Logan phantom, a mouse tissue density and a human jaw density. We analyze the accuracy of reconstruction and the sensitivity to the ROI radius. Finally, we give concluding remarks in Sec. 7.
2. X-ray projections and the Tuy condition. We consider classical projection operators mapping density functions with domain in R 3 into linear projections defined on appropriates spaces of rays. The most prominent examples of such projection operators are the ray transform and the cone-beam transform [24]. Recall that a rayr(u, θ) in R 3 is a line passing through the point u ∈ R 3 and parallel to the vector θ ∈ S 2 , where S 2 is the unit sphere of R 3 . That isr(u, θ) = {u + tθ : t ∈ R}. A half-ray r(a, θ) in R 3 is a half-line originating at the point a ∈ R 3 and parallel to the vector θ ∈ S 2 . That is r(a, θ) = {a + tθ : t ≥ 0}.
2.1. The ray transform. The ray transform maps a function f ∈ L 1 (R 3 ) into its linear projections Xf obtained by integrating over rays at various locations and orientations, that is, for u ∈ R 3 and θ ∈ S 2 . Since Xf (u, θ) does not change if u is moved parallel to θ, it is sufficient to restrict u to the plane through the origin that is orthogonal to θ in R 3 , henceforth denoted by T (θ). Thus, Xf is a function on the tangent bundle of the sphere that we denote by T = {(u, θ) : θ ∈ S 2 , u ∈ T (θ)}. Note that the pairs (u, θ) and (u, −θ) give the same rayr(u, θ), so that the mapping (u, θ) →r(u, θ) is a double covering of T which can be viewed as a 4-dimensional Riemannian quotient manifold. The associated Riemannian volume element on T is du dQ(θ), where dQ(θ) is the surface area on S 2 and du is the Lebesgue measure on the plane T (θ).
We will consider the action of the mapping X on functions with compact support inside a fixed open ball B ⊂ R 3 of radius ρ centred at the origin. We denote by T B the subset of T associated with the rays passing through B, that is Thus, T B is an open submanifold of T with compact closure in T and the natural Riemannian volume element at (u, θ) ∈ T B is given by du dQ(θ). We denote as L p (T B ), 1 ≤ p ≤ ∞, the standard L p function spaces associated to this Riemannian volume.

2.2.
The cone-beam transform. The cone-beam transform maps a function f ∈ L 1 (R 3 ) into the function Df defined by for a ∈ R 3 and θ ∈ S 2 . Here, we view a as the source of the half ray r(a, θ) with direction θ. Hence, Df is a function on the space of the half-rays R = {(a, θ) : θ ∈ S 2 , a ∈ R 3 }. The space R has the structure of a smooth 5-dimensional Riemannian manifold, with natural local coordinates defined by a ∈ R 3 and standard spherical coordinates on S 2 . In particular R has a Riemannian volume element dµ = da dQ(θ), where dQ(θ) is the surface area on S 2 and da is the Lebesgue measure on R 3 .
We assume that all the density functions f have compact support inside an open ball B ⊂ R 3 of radius ρ centered at the origin. In the more realistic tomographic setups considered below, sources are located on a smooth bounded curve Γ ⊂ R 3 supported outside the ball B. We denote by R B the subset of all half-rays with sources on Γ and actually intersecting B that is We call R B the set of active rays. R B is a 3-dimensional manifold of class C ∞ with natural local coordinates defined by the arclength parametrization t of the curve Γ and the standard spherical coordinates on S 2 . Thus, R B is a submanifold of R with Riemannian volume element at (a, θ) ∈ R B given by dt dθ), where dt is the Lebesgue measure on R. The total volume m(B) = µ(R B ) is clearly finite. When the ray sources are located on a piecewise smooth curve Γ exterior to a bounded open ball B, the classical Tuy condition on Γ and B [27,24] ensures that every smooth function f with compact support included in B can be recovered from its non-truncated cone-beam projection Df .
Definition 2.1. Let B be an open ball of finite radius centered at the origin and Γ ∈ R 3 \B be a C ∞ curve of length L parametrized by t → γ(t) ∈ R 3 , for 0 ≤ t ≤ L, with non zero velocities γ (t). The curve Γ is said to verify the strong Tuy condition with respect to the ball B if, for every point x ∈ B and any unit vector θ ∈ S 2 , there is on Γ a source point Γ(x, θ) with C 1 -dependence on (x, θ), such that (3) θ, Γ(x, θ) = θ, x and θ, Γ (x, θ) = 0.
We will then write Γ( By the Implicit Function Theorem, the function λ is then of class C ∞ .
Geometrically, the last condition means that any plane H intersecting B must also intersect the curve Γ non tangentially at least once. Note that each H can have more than one intersection with Γ. For a given ball B in R 3 , and when Γ is a circular helix 'surrounding' B, the strong Tuy condition is satisfied provided the helix is long enough. For the helix cone beam acquisition setup, tangential intersections play a negligible role in Katsevich cone beam inversion formula [13]. When Γ is the union of two concentric circles of equal radius r, positioned on orthogonal planes in R 3 surrounding a given ball B, the Tuy condition is again satisfied for r large enough.
As mentioned above, we will consider in Sec. 6 discrete applications of the conebeam transform for different practical acquisition setups including the spherical case, where Γ is a sphere surrounding the target ball B, the spiral case, where Γ is a segment of circular helix, the C-arm case, where Γ is a circular arc. In all these cases there is a formula to reconstruct a compactly supported smooth density function f from its non-truncated projections. We also consider the twin orthogonal circles case, where Γ is composed of two concentric circles positioned on orthogonal planes in R 3 ; in this case we use a classical approximate reconstruction formula for non-truncated reconstruction.
3. Reconstruction from non-truncated projections. For the non-truncated projection operators considered above mapping density functions in R 3 into a full set of linear projections, it is possible in many cases to define a formal inverse operator.
For the non truncated ray transform X, when f is in the space S of functions on R 3 having fast decreasing derivatives of all orders and when Xf (u, θ) is known for all (u, θ) ∈ T , then there exists an explicit inverse operator Z such that f = Z.Xf (cf. [24,Sec. 2.2] or [9]). As noted by Helgason in [9, Ch. 1, Sec. 7.B], for any fixed convex compact set K, to reconstruct f (x) for all x in the complement R 3 \ K of K, the formula f = Z.Xf requires to know Xf only for rays which do not meet K. Of course this does not solve the reconstruction problem for X-ray data truncated by a convex compact ROI.
For the non truncated cone-bean transform D, if the source location Γ is a piecewise C ∞ curve exterior to a ball B and verifying strong Tuy condition, then for all f in C 2 (R 3 ) with compact support included in B, there is an inverse operator Z such that Z.Df (x) = f (x), and Z can be implemented by one of several variants of Grangeat's formula [7,24]. For example, in spiral tomography, where Γ is a segment of a circular helix, the inverse Z of the non truncated cone beam transform D can be computed either by a variant of Grangeat's formula or alternatively, provided the data are axially truncated, by the Katsevitch's formula [14,33]. In the setting of C-arm tomography, where Γ is an arc of circle, an approximate inverse operator Z of the non truncated cone beam operator D can be computed using again a variant of Grangeat's formula [27,36].
We point out that all these exact formulas inverting the non truncated cone beam operator require smoothness conditions on f to reconstruct f as a function.
As noted by Natterer [24], Tuy [27] and other authors, to define the most generic linear operator Z inverting the transform f → Df , one should consider f and g = Df as distributions instead of functions. However, numerical reconstructions of f from discretized projection data Df usually smooth the non truncated data Df before reconstruction. Therefore classical proofs of exact inversion formulas for non truncated cone beam data tend to focus on smooth density functions f . Indeed, in spiral tomography, where Γ is an helix, the original proof of Katsevich's inversion formula in [14] requires f ∈ C ∞ 0 (B); finite degree of smoothness can be achieved using more sophisticated arguments [15]. Similarly, when Γ is a smooth curve, the proof of Grangeat's inversion formula in [7] requires f ∈ C 2 (B). Also in the more academic setting of the ray transform, where the full set of projections (for all (u, θ) ∈ T ) is known, the inversion formulas in [24,9] require the Fourier transform of f to decrease rapidly at infinity.
In the following, we will define exact inverses Z of the non-truncated projection operators D as explicit linear operators acting on Sobolev spaces of densities. This definition will be useful to derive important continuity properties of Z.
We start by defining appropriate Banach spaces to handle the space of rays.
3.1. Banach spaces of smooth functions on manifolds. Let M be a Riemannian manifold of class C ∞ with volume element dµ and finite volume µ(M).
Here we consider only manifolds which are either compact or are the interior of a compact manifold with a C 1 -boundary. One can then find and fix a finite covering of M by open relatively compact sets U j , j ∈ J endowed with diffeomorphic local maps h j : V j → U j , where the V j are bounded open balls in R 3 , and each h j is the restriction to V j of a local map defined on an open Euclidean ball containing the closure of V j . Explicit such finite coverings U j are easily specified for the manifolds of rays R B given by (2), and for Γ × S 2 , where Γ is either a piecewise C ∞ bounded curve in R 3 or a whole sphere in R 3 . Fix as above a finite covering U j , j ∈ J, of M and the local maps h j : V j → U j . A function g on M is said to be uniformly bounded if and only is all the g • hj are bounded. For any r > 0, call C r (M) the space of all functions on M having continuous and uniformly bounded differentials of all orders up to r. For 1 ≤ p ≤ +∞, we denote L p (M) the usual Banach spaces of functions g on M such that |g| p is µ-integrable and µ is the Borel measure on M.

3.2.
Inversion of the non-truncated ray transform. An analytic inversion formula for the ray transform can be derived from the classical Fourier slice theorem.
In this section, we derive explicit continuity properties for this inversion formula. For any θ in S 2 , let T (θ) be the plane orthogonal to θ in R 3 and containing the origin. As seen in Sec. 2.1, any ray r(u, θ) is non-ambiguously indexed by θ ∈ S 2 and u ∈ T (θ) and the tangent bundle T of the unit sphere in R 3 is a 4-dimensional manifold with volume element du dθ.
Fix an open ball B of radius ρ in R 3 and let T B be the manifold of all rays intersecting B. For any function g(u, θ) on T B , let g θ be the function defined on T (θ) by g θ (u) = g(u, θ). For v ∈ T θ , the 2-dimensional Fourier transform of g θ on the tangent plane T (θ) is given by whenever the integral is well-defined. There is then a constant c = c(ρ) depending only on ρ such that for all g ∈ W 4 (T B ), for all θ ∈ S 2 , v ∈ T (θ), one has the inequality To prove this, note first that, when B is the unit ball U , then equation (5) is a completely standard inequality, holding for a universal numerical constant c = γ. Now for a generic ball B, and for any g for all x. This rescaling immediately entails and . Using these two last equations, and applying (5) to the function G, one easily concludes that inequality (5) For f ∈ L 2 (B), the usual 3-dimensional Fourier transform of f will be denoted by By the Fourier slice theorem (cf. [24,Sec. 2.2]), for any θ ∈ S 2 and z ∈ R 3 such that z, θ = 0, the ray transform g = Xf of f verifies provided the two Fourier transforms involved in the formula are well defined. As shown in [24,9] whenf (z) tends to zero at infinity faster than any polynomial in z, then equation (6) can be used to derive inversion formulas to reconstruct f from its non-truncated projections g.
Fix any ball B in R 3 and let T B be the submanifold of all rays intersecting B. Note that any function g = g(u, θ) in W 4 (T B ) can be extended to T by setting g = 0 on T \ T B . For all θ ∈ S 2 , recall that g θ is the function defined for u in the plane T (θ) orthogonal to θ by g θ (u) = g(u, θ). For the non-truncated ray transform X, the following proposition specifies how to construct a bounded linear inverse Z of X defined on W 4 (T B ).
Proposition 2. Let X be the non-truncated ray transform in R 3 . Select and fix a Borel measurable function z → θ(z) from R 3 to S 2 such that z, θ(z) = 0 for almost all z ∈ R 3 . Fix a ball B, define T B as above and recall that g θ is the function defined for u in the plane T (θ) orthogonal to θ by g θ (u) = g(u, θ). For almost all z ∈ R 3 , the plane T (θ(z)) must contain z, the function h = g θ(z) is defined on T (θ(z)), and for all v in T (θ(z)) the 2-dimensional Fourier transform F g θ(z) (v) of h is given by formula (4), which thus specifies the value F g θ(z) (z). One can then define the function Jg on R 3 by the following integral, which is necessarily finite for all vectors x ∈ R 3 , The restriction Zg = 1 B Jg of Jg to the ball B defines then a bounded linear operator Moreover, for any f ∈ W 4 (B), the non-truncated ray transform g = Xf verifies the identity f = ZXf .
Proof. Let c = c(ρ) be the constant involved in inequality (5). Due to (5), for all g ∈ W 4 (T B ) and all x ∈ R 3 , the integral Jg(x) defined by equation (7) is uniformly bounded by Hence g → Zg is a bounded linear operator from W 4 (T B ) into L ∞ (B) and hence also into L 2 (B). Moreover, for any f ∈ W 4 (B), the function g = Xf is in W 4 (T B ). Therefore the Fourier slice formula (6) combined with (7) shows that

3.3.
Inversion of the non-truncated cone-beam transform. We now construct an operator inverting the non truncated cone-beam transform D when the projections belong to a Sobolev space of rays. Fix a ball B and define T B as above. Let Γ ⊂ R 3 be a C ∞ curve with support exterior to the open ball B and parametrized by γ : [0, L] → R 3 . Assume that the strong Tuy condition is verified, i.e., there is a (3). Consider any function g ∈ C 2 (R B ) with compact support inside R B . We extend g to Γ × S 2 by setting g = 0 on Γ × S 2 \ R B and then extend g to a C 2 function G defined on Γ × R 3 by and, hence, for all x ∈ B, we define the function Zg by where, as above, T (θ) ⊂ S 2 is the set of all α ∈ S 2 such that α, θ = 0.
We have the following result.
. Moreover there is a constant c depending only on Γ and the radius of B such that, for all g ∈ C 2 (R B ) with finite Sobolev norm g W 4 (R B ) , we have In particular Z can be extended to a bounded linear operator from W 4 (R B ) into L ∞ (B) such that whenever g = Df is the non truncated cone-beam transform of f ∈ W 4 (R B ), one has the identity f = Zg = ZDf .
Proof. When g = Df with f ∈ C 2 (B), the assertion Zg = f is proved with different notations in [24, Sec. 5.5.2] using a variant of the Grangeat's inversion formula due to Zeng, Clack and Gullberg [35]. For a generic g in C 2 (R B ), the vector valued function K, given by (10), is continuous by construction and hence remains bounded in In the following, c, c 1 , c 2 , . . . denote positive constants which depend only on the radius of B and Γ but not on g.
In equation (11), the denominator den(x, θ) = θ, γ(λ(x, θ)) is continuous for x ∈B, θ ∈ S 2 and is never zero due to Tuy conditions, so that |den| ≥ c > 0 for some constant c. Then equation (11) readily provides a constant c 1 such that Set h(t, θ) = g(γ(t), θ). Equations (9) and (10) show that there is a constant c 2 such that for all g in C 2 (R B ), for all t ∈ [0, L], α ∈ S 2 . By definition of W 4 (R B ), there is a constant c 3 such that for any function g in W 4 (R B ) the function h verifies The Sobolev imbedding theorem in dimension 3 holds on the Riemannian manifold R B , relating the norm of h in C 2 (R B ) with its norm in W 4 (R B ), as explained in Appendix. It thus provides a constant c 4 such that, for any function h ∈ W 4 ([0, L]× S 2 ), all partial differentials d dt ∂ θ h of order ≤ 2 of h are bounded and continuous on [0, L] × S 2 and verify (16) sup Combining the inequalities (13) (14) (15) (16), we conclude that Remark 1. As mentioned above, in spiral tomography, the source curve Γ ⊂ R 3 is a circular helix, which can be parametrized as γ(t) = (a cos(t), a sin(t), kt) for some fixed positive a, k. Let f be a compactly supported density function f ∈ C ∞ 0 (B), with ρ < a (so that the helix is surrounding the support of f ). Katsevich proved [12] that, in this setting, f can be reconstructed from its non-truncated cone-beam projections Df by a formula which can be written as where U and V are explicit operators involving divergent integrals. This divergence is carefully analyzed by adequate approximations in [12] but this inversion formula is rather unwieldy and the continuity properties of (U + V ) are not easy to evaluate directly.
As mentioned above, even though an helix does not satisfy strong Tuy condition in general, this condition holds for a bounded circular helix as long as the support of the density is sufficiently small. Hence, we can define an explicit inverse operator Z of the non-truncated operator D according to Proposition 3. Our approach has the ability to reconstruct f from projection data Df ∈ W 4 (R B ) for all f ∈ W 4 (B) and to provide a precise Sobolev continuity property for the inverse operator Z.

4.
Regularization in the space of rays. We now define and construct a class of regularization operators in the space L 2 (R B ). We start by defining a notion of approximations of the identity on the manifolds considered in Sec. 3.
Definition 4.1. Let M be a C ∞ Riemannian manifold with volume element dµ and finite volume. As above, assume that M is either compact or is the interior of a compact manifold with a C 1 -boundary. For any integer r ≥ 1, we call C r approximation of the identity in L 2 (M, µ) any sequence of linear operators τ N : L 2 (M) → C r (M) verifying the following conditions.
(i) There is a constant c such that for all g ∈ L 2 (M) and all integers N (ii) For any g ∈ L 2 (M), lim N →∞ g − τ N g L 2 (M) = 0.
(iii) For each integer 2 ≤ p ≤ (r + 1) there is a constant c such that for all g ∈ W p (M), g − τ N g W p−1 (M) ≤ c/N g W p (M) .
(iv) Whenever g has compact support then τ N g also has compact support.
Note that N → τ N will remain a C r approximation of the identity in L 2 (M, ν) for any positive Borel measure ν on M such that both densities dν dµ and dµ dν are bounded.
We next show how to construct C r approximations of the identity in L 2 (M, µ).

4.1.
Approximation of the identity by small convolutions. When the manifold M is a bounded open Euclidean ball in R k , one can generate a C r approximation of the identity as follows. Select any fixed C r function w ≥ 0 on R k with compact support and Lebesgue integral equal to 1 and, for f ∈ L 2 (R k ), define the "small" convolutions by . Standard results on convolutions show that the sequence (σ N ) verifies the properties (i),(ii) and (iv) of Definition 4.1. The proof of property (iii) is the following. For z ∈ R k and 1 ≤ p < ∞, we define φ p (z) = (1 + |z| 2 ) p/2 . Denoting by W N the Fourier transform of w N , we have that W N (z) =ŵ(z/N ), whereŵ is the Fourier transform of w. Hence, for all z ∈ R k , (18) |W for all z ∈ R k . From the last inequality and by the definition of Sobolev norms, we then get: This proves property (iii). Small convolutions are easily and explicitly extended to the manifolds of active rays R B , as we now show by patching together local small convolutions through appropriate local maps.
Proposition 4. Let M be a C ∞ Riemannian manifold of finite volume. As above, assume that M is either compact or is the interior of a compact manifold with C 1boundary. Then, for any integer r, one can construct explicitly a C r approximation of the identity N → τ N on L 2 (M).
Proof. As observed above, on M we can select an open finite covering U j , j ∈ J, and local maps h j : V j → U j to construct a finite partition of unity by C ∞ functions u j with compact supports included in U j and verifying 0 ≤ u j ≤ 1 and j∈J u j = 1. On each Euclidean ball V j , select a C r approximation of the identity N → σ N (j) in L 2 (V j ), for instance by small convolutions as indicated above.
For any g ∈ L 2 (M), let g j = g u j and define the operators g → τ N g by (19) Figure 1. ROI-truncated cone-beam acquisition: projections are restricted to half-rays intersecting the ROI, which is a ball C included in the target ball B.
where G j = σ N (j)(g j • h j ) Since each map h j can be smoothly extended to a neighborhood ofV j , each h j has bounded derivatives of any order. Hence the mapping g j → g j •h j is a bounded linear operator from L 2 (U j , µ) to L 2 (V j ), and from W r (U j ) to W r (V j ). Similar boundedness properties hold for the linear operators mapping G j → G j • h −1 j and g → gu j . The explicit formula (19) and the fact that the σ N (j) are C r approximations of the identity in L 2 (V j ) then implies directly that the sequence (τ N ) satisfies the properties (i)-(iv) of Definition 4.1. Hence (τ N ) is a C r approximation of the identity in L 2 (M, µ).
Proposition 4 applies in particular to manifolds of active rays R B associated to an open ball B and a smooth set of sources Γ exterior to B.

Reconstruction from ROI-truncated projections.
Fix a bounded open ball B ⊂ R 3 and a smooth set of sources Γ verifying strong Tuy condition. Let C be a spherical ROI strictly included in B. As illustrated in Figure 1, the C-truncated cone-beam transform D C f of f is the restriction of Df to the manifold R C ⊂ R B of half-rays which intersect C. Denoting the indicator function of a set G by 1 G , one can then write by Proposition 1, D C is a bounded linear operator from L ∞ (C) into L ∞ (R C ). As observed above, the non-truncated operator D can be inverted by a linear operator Z such that f = ZDf for all densities f ∈ W 4 (B) having compact support. However the C-truncated operator D C cannot be directly inverted by applying Z to D C f , even within the region C. We formally define 'approximate' inverses Z C for D C . For any > 0, we say that the ROI-truncated cone-beam transform D C admits an for all f ∈ W 5 (B).

5.1.
Reconstruction from ROI-truncated projection data. Let B ⊂ R 3 be a spherical region. Let Γ be a smooth curve surrounding B and verifying the strong Tuy condition. As seen above, the non truncated cone-beam transform D with sources on Γ has then an inverse Z : . The theorem below shows how to obtain an accurate inverse of the ROI truncated operator D C .  with f 0 = Zτ g. Then as j → ∞ the sequence (f j ) converges at exponential speed in L ∞ (B) to a limit Z C g. This defines a bounded linear operator Z C from L ∞ (R B ) into L ∞ (B). For any unknown density f in W 5 (B), one can apply Z C to the ROItruncated data g = D C f to compute an approximation of f given byf = Z C D C f . Finally for η( ) small enough,f will verify That is, Z C is an -accurate inverse of the ROI truncated cone beam operator D C .
Remark 2. According to Theorem 5.2, the truncation region C is strictly included in B and must be large enough for the theorem to hold. The proof given below gives an explicit upper bound η( ) for rad(B) − rad(C) guaranteeing the existence of an -accurate inverse Z C for the C-truncated cone beam projector D C . Our theoretical upper bound for η( ) is proportional to 8 , which is an overly pessimistic estimate. The main reason for this is that inequality (23) controls the L ∞ norm of the reconstruction error f −f on the whole sphere B rather than controlling the error only within the ROI sphere C ⊂ B. We conjecture that the reconstruction error on the ROI C is already of the order of as soon as rad(B) − rad(C) is of the order of 2 . Indeed, in all the practical situations we have studied numerically (see Sec. 6), our simulations indicate that, for = 0.10 one could actually take η( ) = rad(B)/2, so that for all spheres C ⊂ B with same center as B and radius rad(C) > rad(B)/2, our ROI-reconstruction algorithm (22) does converge to an -accurate inverse Z C of the truncated cone beam projector D C . Such situations are thus clearly quite favorable for radiation exposure reduction by ROI truncated data acquisition combined with our ROI reconstruction algorithm. As seen above, -accurate inverses are generally not unique. For instance the following variant of theorem 5.2 outlines alternative constructions of Z C , where two regularization operators are involved, one in the space of rays and another one in the space of densities.
We have the following variant of our ROI-reconstruction theorem.
Theorem 5.3. Notations and hypotheses are the same as in theorem 5.2. Fix two C 4 -approximations of the identity: (τ N ) in L 2 (R B ) and (σ N ) in L 2 (B). Given any > ε 0, one can find two regularisation operator τ = τ N () ε and σ = σ N () ε close enough to the identity, and η() ε < rad(B) such that the operator becomes a contraction of L ∞ (B) for all ROIs C verifying (21) Then for each g ∈ L ∞ (B), the sequence (f j ), given by the recurrence (22) with U as just defined and f 0 = σZτ g, converges in L ∞ (B) to a limit Z C g, where Z C is an -accurate inverse of the ROI-truncated cone beam transform D C .

Proof of the Theorem 5.2.
Before presenting the proofs, we need the following two lemmata. For the remaining of this section, let Γ, B, C, D, D C be given as above.
Lemma 5.4. There is a constant c determined by Γ and B only such that, for any sphere C ⊂ B and any f ∈ L ∞ (B), the linear operator where rad(C) is the radius of C.
Proof. Let s be any source position on the curve Γ. Denote by z(C) the center of C. Call H(s, C) ⊂ S 2 the set of all θ ∈ S 2 such that the half-ray r(s, θ) intersects C. The set of all these half-rays is a cone of revolution with vertex s, axis [s, z(C)], and half-aperture angle 0 < α(s, C) < π/2. The area of the spherical cap H(s, C) is hence classically given by (25) area(H(s, C)) = 2π(1 − cos(α(s, C)).
From the last equation, using (27), we obtain where L is the length of Γ. Hence we obtain the estimate ). Due to Proposition 1, there is a constant c 0 such that, for any f ∈ L ∞ (B), one has where c 1 = c 0 (4πLc) 1/2 . This achieves the proof when Γ is a curve of length L and does not intersect the closure of B.
As seen above, due to the Tuy condition, the non-truncated cone beam transform D can be inverted through a bounded linear operator Z : W 4 (R B ) → L ∞ (B) given by the Grangeat formula (11). The lemma below shows how to construct a contraction on L ∞ (B) starting from Z and D − D C . Lemma 5.5. Let (τ N ) be a C 4 approximation of the identity in L 2 (R B ) as in Definition 4.1. Then there is a constant k such that for any C ⊂ B with radius verifying rad(B) − rad(C) < k/N 8 , the linear operator Proof. By inequality (24), there is a constant c such that, for all C ⊂ B and all Applying inequality (17) By applying inequality (12) to the function τ N Y C f , we then obtain a new constant c 2 such that, for all C ⊂ B, all f ∈ L ∞ (B) and all N , Set c 3 = c 2 c 1 c. Then, provided rad(B) − rad(C) ≤ c 3 /N 8 , the linear operator U N is a contraction of L ∞ (B) with operator norm U N L ∞ (B) ≤ 0.9.
We can now prove Theorem 5.2. Select and fix a C 4 approximation of the identity N → τ N ) in L 2 (R B ), so that the bounded operators τ N : L 2 (R B ) → W 4 (R B ) verify definition 4.1. We will use the following shorthand notations for the norms of various linear operators T : is the norm of T : Due to Propositions 1 and 3, Definition 4.1 and equation (11), there is a constant c > 0 such that for all integers N Given an > 0, fix an integer N = N ( ) by Since N = N ( ) is now fixed, we will write τ = τ N . By Lemma 5.5, there is a constant k such that for any spherical region C ⊂ B verifying the inequality We now fix where the constant c t is determined by k and c, and we assume that the ROI radius verifies rad(B)−rad(C) ≤ η( ). As just seen above, this forces U to be a contraction with norm inferior to 0.9. Given any g in L ∞ (R B ), define a sequence (f j ) ⊂ L ∞ (B) by the algorithm (22). By construction one has then Hence, as j → ∞, the sequence (f j ) converges at exponential speed in L ∞ (B) to a limit Z C g ∈ L ∞ (B). By (22), we must have This defines a linear operator Z C : L ∞ (R B ) → L ∞ (B). We now show that Z C has bounded operator norm. Since U L ∞ (B) ≤ 0.9, the operator (I − U ) : L ∞ (B) → L ∞ (B) has a bounded inverse given by the converging series ∞ j=0 U j , which yields Equation (31) yields that, for all g ∈ L ∞ (R B ), Hence, due to the bounds (29), where m(B) is the finite Riemannian volume of R B . So the operator norm of Z C is bounded. For We then can write Combining this observation with the bounds given by (29) for the norms of Z, (I − τ ) and D, we obtain, for all f ∈ W 5 (B), From (34) and (35) we then get that for all f ∈ W 5 (B) and, hence, since (I − U ) −1 L ∞ (B) ≤ 10, we conclude that For any f ∈ W 5 (B), the expression of Z C h = Z C D C f given by equation (32) Our choice of N in (30) forces 10 c 3 N ( ) ≡ . Thus, for all f ∈ W 5 (B) and all regions C ⊂ B verifying rad(B) − rad(C) ≤ η( ), we get By definition 5.1, Z C is thus an -accurate inverse of D C . This completes the proof of Theorem 5.2.
Remark 3. The preceding proof introduces the condition rad(B) − rad(C) ≤ η( ) with η( ) ≡ c t 8 . This explicit upper bound is definitely too pessimistic as shown by all our simulations below. The key reason for this is that we actually control the L ∞ reconstruction error on the whole region B instead of controlling it only on a specific ROI C. The main advantage of our strong constraint on rad(B) − rad(C) is to yield a uniform control of L ∞ reconstruction errors for any such C and for our very wide range of regularization operators.

5.2.
Proof of Theorem 5.3. This proof is similar to the argument used in the preceding proof and will just be sketched. Select two C 4 approximations of the identity (σ N ) in L 2 (B) and (τ N ) in L 2 (R B ) and set V N = σ N Zτ N Y C . By applying Lemma 5.5 to V N we see that, given > 0, there is N = N ( ) large enough and η( ) small enough to ensure that, for all regions C verifying rad(B) − rad(C) ≤ η( ), the operator U = V N () ε becomes a contraction of L ∞ (B) with norm less than 0.9. Set then σ = σN () ε and τ = τ N () ε . As j → ∞, the iterative sequence of approximate densities f j+1 = f 0 + U f j initialized by f 0 = σZτ g will then converge in norm to a limit Z C g ∈ L ∞ (B). As in the preceding proof, one checks thatZ C is an -accurate inverse of D C . 5.3. Extensions. Theorems 5.2 and 5.3 can be similarly extended to the case where the set of sources Γ is a sphere strictly including B rather than a curve. The constant c in (24) then depends on the area of the sphere Γ instead of the length of the curve Γ. Our ROI-reconstruction algorithm to approximately invert the ROI-truncated cone beam transforms D C can also be extended using nearly identical arguments to generate approximate inverses of the ROI-truncated ray transforms X C .
For discrete data D C f generated by a numerical ROI-truncated cone beam transform of the unknown density f , our ROI reconstruction algorithm (22) can formally be applied, starting with any discretized reconstruction algorithm Z specific to the acquisition geometry at hand and known to perform well on non-truncated data. From a formal point of view, Z can even be a 'black-box' inversion software dedicated to inversion of non truncated data. Of course the proofs Theorem 5.2 do not a priori cover such cases from a theoretical point of view. Yet, our extensive numerical tests in Sec. 6 indicate that our ROI-reconstruction approach does perform well in many practical CT setups. 6. Numerical experiments. In this section, we present numerical experiments to evaluate the performance of our ROI reconstruction algorithm in multiple discrete settings with sources located on a smooth curve or a sphere. Since we had no access to ROI-truncated data acquired with an actual CT device enabling cone-beam ROI truncation, our numerical experiments simulated ROI-truncated cone-beam acquisition. We used four classical acquisition geometries and multiple spherical ROIs with different locations and sizes applied to several 3D density data. The goal of these numerical experiments was to numerically quantify accuracy of our ROI reconstruction algorithm within the ROI and to investigate how the ROI radius impacts this measure of accuracy.
For a given cone-beam acquisition setup with target ball B ⊂ R 3 , a point z in B and a number ν > 0, our Theorems 5.2 and 5.3 imply the existence of a critical radius ρ such that, for any spherical region C ⊂ B with center z and radius rad(C) > ρ, and for any density f with f L ∞ (B) ≤ ν, our ROI reconstruction algorithm from truncated data will converge within C to a good approximation of the unknown f . Our numerical experiments provide practical evaluations of this critical radius ρ.
6.1. Simulations of ROI-truncated acquisition. We have simulated ROItruncated cone-beam acquisition for three discretized densities f : a 3D Shepp-Logan phantom; a 3D scan of mouse tissue; a 3D scan of a human jaw. For each discretized density f , given by a 3D image of size 256 3 voxels, we have first computed discrete non trucated cone-beam projections Df by simulating discrete acquisition and used these data to generate ROI-truncated 3D projections for four distinct acquisition geometries with the following parameters: (i) Spherical tomography with sources on a full spherical surface. Ray discretization: 3 degrees in the polar direction, 5 degrees in the azimuthal direction; scanning radius: 400 voxels; number of detector rows: 256; source-detector distance = 900 voxels. For each one of these four acquisition setups, we selected four concentric spherical ROI C with radius values (in voxels) equal to 45, 60, 75, 90. Each such ROI C was used to truncate the discretized projection data Y = Df to the rays intersecting C and thus to generate a discretized version of the ROI-truncated data Y C = D C f . Note that non truncation corresponded to a much larger spherical radius (221 voxels) covering the entire 3D density volume.

Numerical implementations of our ROI reconstruction algorithm.
For each 3D discrete density function f , each cone-beam acquisition setup and each spherical ROI C, we have implemented our iterative ROI reconstruction algorithm to compute a reconstruction Z C f of f using only the ROI-truncated data D C f . We first select a regularization operator σ acting on the space of L 2 densities, and close to the identity. We then choose an explicit discrete inverse Z of the nontruncated cone-beam transform, and Z will naturally depend on the acquisition setup. According to Theorems 5.2 and 5.3, we then set U = σZ(D − D C ) and we apply our iterative ROI reconstruction formula f j+1 = f 0 +U f j , with f 0 = σZD C f . If C is large enough within B, we expect to observe numerical convergence of the f j to an ROI-reconstruction Z C f of f .
For each one of our simulated cone-beam acquisition setups, the inverse Z of the non-truncated cone-beam transform was implemented as follows: (i) Spherical tomography: Inversion of non-truncated cone-beam transform by filtered back-projection (FBP) [24]. (ii) Spiral tomography: Inversion of non-truncated cone-beam transform by Katsevich inversion formula [14]. (iii) C-arm tomography: Inversion of non-truncated cone-beam transform by a standard FDK algorithm [6]. (iv) Two circles tomography: Inversion of non-truncated cone-beam transform by a discretized Grangeat formula outlined in Defrise and Clack [4].
6.2.1. Choice of a regularization operator σ on the euclidean ball B. The generic ROI-reconstruction scheme covered by our theoretical results involves the use of linear regularization operators σ : L 2 (B) → W 4 (B) such as convolutions by smooth kernels with compact support. These small convolution operators can be well approximated by Gaussian convolutions and, in a few of our simulated acquisition setups, we did validate the possibility of implementing σ by an adequately discretized Gaussian convolution kernel. However, for computational efficiency and for better overall reconstruction performance, in all of our simulated acquisition setups, we have preferred to implement discretized versions of the regularization operator σ : L 2 (B) → W 4 (B) based on wavelet thresholding, even though thresholding is clearly not a linear operation. In [26,Ch.8], we have examined several regularization techniques and found that wavelet hard thresholding produces the best accuracy among the methods we considered. So to compute σh for any h in L 2 (B), we first expand h using standard Daubechies wavelets Daub4 [21] in R 3 to generate the wavelet decomposition of h h = m,n,i a(m, n, i)φ m,n,i .
Each wavelet φ m,n,i in this family is a C 4 function indexed by the discretized position (m, n) of its compact support and by an integer scale parameter i ≥ 0. At the coarsest scale i = 0, no truncation or shrinkage was applied to the wavelets coefficients a(m, n, i). At finer scales i ≥ 1, the wavelets coefficients a(m, n, i) were set to zero whenever |a(m, n, i)| < T HR i . where the thresholds T HR i were selected to discard 90% of wavelet coefficients. We found that the performance of the algorithm is not very sensitive to the choice of the percentage of discarded wavelet coefficients, that is, performance would remain essentially the same by discarding 70-90% of wavelet coefficients. Additional details can be found in [26].
The new wavelet expansion generated by this coefficients truncation defined the function σh which obviously belonged to W 4 (B). This operator σ is non linear but can be well approximated by linearized versions which implement smooth shrinkage of the wavelets coefficients instead of abrupt truncation [26].

6.2.2.
Stopping rules for our ROI reconstruction. Let C be the spherical ROI. As a stopping criterion, we adopted a standard rule so that the algorithm (22) is to stop the iteration over the index j when f j and f j+1 become close enough within C; in particular, as long as f j+1 − f j L 1 (C) ≤ b for some small tolerance b, e.g. b = 0.02. We automatically stop the ROI iterative reconstruction at j = 40 to avoid unnecessary computation as we found that, for all our numerical experiments, as soon as the radius of C was slightly superior to a critical radius, 40 iterations were amply sufficient to achieve convergence. 6.3. Performances of numerical ROI reconstructions. For each one of our three densities f , each one of our four x-ray acquisition setups, and each one of our selected spherical ROIs C, our simulations generated a discrete version g = D C f of the ROI-truncated cone-beam projections of f . Then the numerical application of our iterative ROI reconstruction algorithm to these truncated data g provided a discretized approximation Z C f of the "unknown" f . To assess the accuracy of our discretized ROI reconstruction Z C f , we have evaluated an ROI Relative L 1 Error of Reconstruction within C defined by the following ratio RL 1 of two discretized L 1 (C) norms For each one of the 48 ROI reconstruction cases indicated above (each density data set has size 256 3 ), we have recorded the Relative L 1 reconstruction error computed within the ROI in Table 1. As indicated above, our iterative algorithm (22) uses different numerical routines to implement the non-truncated inverse operator Z depending on the acquisition geometry. The number of iterations needed to achieve convergence of our ROI reconstruction algorithm was bounded above by 40 but the algorithm was found to converge (according to the stopping rule given above) with a much smaller number of iterations, typically between 10-12 iterations for sources on a curve. For each one of the various combinations of density data and acquisition setup, these accuracy results yield an estimate of the critical ROI radius ρ enabling a relative ROI reconstruction accuracy inferior or equal to 10%. Such critical radius estimates are displayed in Table 2.
The best performances of our ROI reconstruction algorithm naturally occur for spherical acquisition. Indeed for the somewhat academic spherical setup, the number of projections available is much larger than for the three other setups where sources are located on a curve.
For the twelve situations evaluated here, we obtain a range from 52 to 82 voxels for the critical radius ρ yielding a 10% accuracy in ROI reconstruction. This compares very favourably to the maximal ROI radius corresponding to non truncation. Indeed when one goes from non truncation to a critical spherical ROI, the reduction in irradiated volume ranges from 70% to 98%, indicating a quite strong "formal" reduction in x-ray exposure, while the loss in relative reconstruction accuracy is only of the order of 7%.
Note also that the actual critical radius estimates obtained here by simulations are much smaller that the theoretical upper bounds used in the proof of Theorem 5.2.
For a fixed ROI radius, the ROI Relative Reconstruction error is lower for the Shepp-Logan phantom than for Mouse Tissue or Human Jaw density data. Indeed, when the ROI radius is larger than the critical radius, our iterative ROI reconstruction essentially converges within the ROI to a regularization σf of f . The ROI reconstruction error in L 1 (C) can roughly be viewed as the sum of two terms, a 'convergence' error Z C f − σf L 1 (C) and a 'regularization' error f − σf L 1 (C) . To highlight the regularization effect, we have computed the relative 'regularization error' within C given by For an ROI radius of 70 voxels, this regularization error is equal to 1.1% for the 3D Shepp-Logan phantom, and to 2.4% for the Mouse Tissue and Human Jaw 3D data, because the piecewise constant Shepp-Logan phantom density can be approximated by our wavelet-based regularization operator much more effectively than the more textured Mouse Tissue and Human Jaw densities. So for the Human Jaw data the regularization error contributes about half of the relative L 1 reconstruction error.
To illustrate visually the overall performance of our iterative ROI reconstruction from truncated cone-beam data, we have include several examples. For all these examples, the size of the images is 256 3 voxels and the ROI radius is 45 voxels. In Figures 2-3 we show horizontal, coronal, and sagittal planes from our 3D reconstruction of the Shepp-Logan 3D Phantom and Mouse Tissue data using simulated Twin Circle acquisition. Results from our algorithm are compared against the ground truth and a 'naive' reconstruction that is obtained by directly applying the reconstruction formula to the smoothed truncated data and all missing projection data set to zero. The smoothing of truncated data was performed using a Gaussian kernel with standard deviation 0.1. We observed heuristically that larger values of the standard deviation increase the reconstruction error. For better comparison, we also include line profiles of the reconstructions against the ground truth.
In Figures 4-5, we show horizontal sections from the reconstructed volumes of the Shepp-Logan 3D Phantom and Mouse Tissue data using simulated spiral and C-arm acquisitions.
In all these figures, the comparison of results from our iterative ROI reconstructions with the naive reconstruction approach described above and based on the reconstruction formula devised for non-truncated cone-beam data show that, as expected, the one-step inversion formulas for non-truncated data perform poorly when applied to ROI-truncated cone-beam data, and display multiple visual artifacts especially near the ROI boundary. By contrast, our iterative ROI reconstruction results are very satisfactory even for relatively small ROI radii. Compared to the ground truth, our ROI reconstruction shows some blurring which is due to the wavelet-based regularization step. Since our wavelet filters have finite support and length 4 pixels, they have a rather limited impact on spatial resolution. The blurring effect is consistent with our theoretical prediction since our algorithm generates an approximation of the exact solution which is a smoother version of the true image.
In practical situation, acquired projected images are affected by measurement errors generating noise. Even though the accurate modeling of the noise is beyond Figure 2. Visual comparison of ROI reconstruction for 3D Shepp-Logan phantom using simulated Twin Circles acquisition and truncation of projection data. ROI radius = 45 voxels. Middles sections are shown from the xy, yz and xz planes. From left to right: inversion by one-step Grangeat formula; our iterative ROI reconstruction; ground truth. The last column shows intensity profiles corresponding to the middle row of the images. Green: one-step Grangeat formula; blue: our algorithm; red: ground truth. the scope of this paper, it is still natural to ask what is the effect of measurement errors on the stability of the proposed algorithm. The analysis of the stability to noise was carried out in [26,Ch. 9] using a standard Gaussian additive noise model. Results show that our reconstruction algorithm is very stable in handling noisy projection data. 7. Conclusion. This paper studies ROI tomographic reconstruction in R 3 from ROI-truncated cone-beam data, a highly relevant problem in many applications. For both our theoretical and numerical analysis, we have considered fairly generic cone-beam acquisition setups, with sources located on arbitrary bounded smooth curves Γ in R 3 verifying the classical Tuy condition. In these cases, the non-trucated cone-beam transform Df of any C 2 density f with compact support admits an explicit inverse ZDf = f given by Grangeat formula, but the Grangeat operator Z (or its variants) cannot directly reconstruct f from ROI-truncated data.
For this generic cone beam acquisition setup, we have developed and rigorously proved convergence of a new iterative ROI reconstruction from ROI-truncated cone beam data, valid for densities f in the Sobolev space W 5 (B), where B is a bounded ball in R 3 . Figure 3. Visual comparison of ROI reconstruction for Mouse Tissue data using simulated Twin Circles acquisition and truncation of projection data. ROI radius = 45 voxels. Middles sections are shown from the xy, yz and xz planes. From left to right: inversion by one-step Grangeat's formula; our iterative ROI reconstruction; ground truth. The last column shows intensity profiles corresponding to the middle row of the images. Green: one-step Grangeat formula; blue: our algorithm; red: ground truth.
We focus on ROIs which are spherical regions C ⊂ B. The ROI truncated cone beam operator D C involves only the half-rays intersecting C and emitted by sources located on the curve Γ. Let Z be an exact stable inverse of the non truncated cone beam projector D, such as the operator given by Grangeat formula. After a proper choice of the regularization operator (or mollifier) τ acting on L 2 (R 3 ), we introduce the linear operator U = Zτ (D − D C ) and we recursively compute the densities f j+1 = f 0 + U f j where f 0 = Zτ D C f . For each > 0, we show the existence of a critical radius ρ < rad(B) such that, provided ρ < rad(C) and τ is close enough to the identity, the sequence f j converges to at exponential speed in L ∞ (C) to a limitf = Z C f verifying Our construction thus defines a linear operator Z C on W 5 (B) which we call anaccurate inverse of the ROI-truncated cone beam projector D C . Such approximate inverses are not unique, since we only look for an approximate solution consistent with the ROI-truncated data D C f . . Visual comparison of ROI reconstruction for 3D Shepp-Logan phantom and mouse tissue using simulated spiral acquisition and truncation of projection data. A representative horizontal section from the 3D reconstructed volume is shown. From left to right: inversion by one-step Katsevich formula; our iterative ROI reconstruction; ground truth.
In our introduction, we have discussed other existing iterative methods for ROI reconstruction in CT, but to the knowledge of the authors, no published theoretical result seems to have proved in R 3 and for X-ray sources on a generic smooth curve the existence of a critical ROI radius ensuring the convergence of an iterative ROI reconstruction scheme.
A practical implementation of our ROI-reconstruction is straightforward and does not require any a priori knowledge of the W 5 Sobolev norm of the unknown density f . We have implemented our method for discontinuous densities such as the Shepp-Logan phantom in R 3 and numerically validated our theoretical results using simulated 3D acquisition of ROI-truncated cone-beam data for four classical acquisition geometries (spherical, spiral, circular arm, twin orthogonal circles), combined with three different density functions in R 3 . For each one of these 12 configurations we have studied multiple spherical ROIs, with various radii and locations. This led us to examine about 120 numerical implementations of our ROI reconstruction in 3D. Since our theoretical upper bound estimates for the critical radius ρ < rad(B) were deliberately conservative in order to achieve simpler proofs, we had to determine ρ numerically in each one of our 12 benchmark configurations. This was done for = 0.10 and for all these 12 examples, we did have ρ(0.10) inferior to rad(B)/2. Our numerical results thus confirm that strong reduction of radiation exposure and good accuracy of ROI reconstruction are simultaneously achievable by ROI truncated acquisition combined with our iterative ROI reconstruction. How to obtain Figure 5. Visual comparison of ROI reconstruction for 3D Shepp-Logan phantom and mouse tissue using simulated C-arm acquisition and truncation of projection data. A representative horizontal section from the 3D reconstructed volume is shown. From left to right: inversion by one-step FDK algorithm; our iterative ROI reconstruction; ground truth. an accurate theoretical estimate for the critical radius ρ is still open and will be the object of future work.
Appendix: Sobolev imbeddings. In Section 3.1, to show the regularity of the linear operator Z, we make use of Sobolev imbedding theorems. We quote a special case of a result by Aubin (Theorem 2.34 in [2]) in this context. and this imbedding is compact if k − α > n/2.
By the compactness of the embedding, it is also continuous, meaning that if k − α > n/2 then there exists c > 0 such that for each f ∈ W 4 (M), In particular, if n = 3 and k = 4, then we can choose α = 2 and have the following result.