Efficient tensor tomography in fan-beam coordinates

We propose a thorough analysis of the tensor tomography problem on the Euclidean unit disk parameterized in fan-beam coordinates. This includes, for the inversion of the Radon transform over functions, using another range characterization first appearing in [Pestov-Uhlmann, IMRN 2004] to enforce in a fast way classical moment conditions at all orders. When considering direction-dependent integrands (e.g., tensors), a problem where injectivity no longer holds, we propose a suitable representative (other than the traditionally sought-after solenoidal candidate) to be reconstructed, as well as an efficient procedure to do so. Numerical examples illustrating the method are provided at the end.


Introduction
The present work proposes a detailed account of the tensor tomography problem (TTP) on the Euclidean unit disk in fan-beam coordinates, including inversion procedures (modulo kernel) and full range description of the operator over tensors of arbitrary order.
The tensor tomography problem, posed on a non-trapping Riemannian surface (M, g), consists of (i) determining what part of a symmetric m-covariant tensor f (m ≥ 0) is reconstructible from knowledge of its geodesic X-ray transform (XRT) If (the collection of its integrals along all geodesics passing through that surface), and (ii) how to reconstruct a faithful representative of f from If modulo the kernel of I. Such a problem arises for its applications to imaging sciences as well as for its ties to integral geometric problems such as spectral and boundary rigidity, and Calderón's inverse conductivity problem. Particular applications are Computerized Tomography (m = 0), Ultrasound Doppler Tomography (m = 1, see [1]), deformation boundary rigidity (m = 0, 2, see [41]) and imaging of elastic media with slightly anisotropic properties (m = 4, see [39]).
Regarding question (i), for any m ≥ 1, this problem is famously non-injective, as the inner derivative operator d (see [39] for details) generates an element in the kernel of the XRT out of any symmetric tensor vanishing at the boundary of M . In this regard, a way of reformulating question (i) is to ask whether these elements are the only ones in the kernel, a question which received positive answers in the case of surfaces of revolution satisfying Herglotz' condition in [40] (in fact, all dimensions ≥ 2 there), and in the case of simple 1 surfaces in [31]. Neither case is a subset of the other though both cover the Euclidean disk, studied in the present article. Let us mention that further dynamical tools have recently allowed results on cases with trapped geometry [11,12,13], Fourier analysis and representation theory have led to progress on Radon transforms on compact Lie groups and finite groups [20,18,19], and microlocal analytic approaches have led to results on general manifolds with conjugate points [27,43,44,45,47,17].
On studying question (ii), the solenoidal-potential decomposition of an m tensor f = dh+f s (with h an m − 1 tensor vanishing on ∂M ), true for any L 2 tensor f and where each summand is continuous in terms of f , suggests that, since If = If s , the "solenoidal" tensor f s is a good candidate as a representative to be reconstructed from data If . In this direction, Sharafutdinov provided a reconstruction of f s in Euclidean free space in any dimension in [39,Theorem 2.12.2], Kazantsev and Bukhgeim proposed in [22] a reconstruction algorithm for a tensor of arbitrary order in the case of the Euclidean unit disk, see also [21,46,8,6] for further works on the topic and the recent work [7]. Another approach is to make use of the theory of A-analytic functionsà la Bukhgeim, which has successfully led to range characterizations and inversion formulas, most recently in [36,38,37] for the ray transform in Euclidean convex domains over functions, vector fields and two-tensors, including attenuation coefficients.
The solenoidal representative is arguably not an easy quantity to work with: Sharafutdinov's formulas involve iterated use of the non-local operator ∆ −1 d 2 , where ∆ denotes componentwise Laplacian, and the expression of d 2 depends on the tensor order. A first salient feature of this paper is to exploit another decomposition, arising naturally for instance when performing Fourier analysis on tangent fibers: doing this comes with a different inner product structure than the one traditionally used on tensor fields, and in particular, this motivates another tensor decomposition than the solenoidal-potential one. This decomposition has proved to be more natural in earlier contexts such as the search for conformal Killing tensor fields (see for instance [41,42] or [5,Theorem 1.5]). Here we exploit this other decomposition to provide, for any tensor field, an element to be reconstructed with some advantages: its Fourier components are made up of a function to be inverted for via "traditional" inverse X-ray transform, plus additional terms which can be easily expressed in terms of analytic functions, the reconstruction of which via ad hoc Cauchy integrals is straighforward and efficient; in addition, the X-ray transforms of each component live on L 2 -orthogonal subpaces of data space, and the reconstruction formulas given below encode projection onto each subspace before reconstruction, without requiring intermediate processing. In the case of tensors of odd order, special care is given to the reconstruction of solenoidal one-forms supported up to the boundary, for which inversions in [32] only covered the case of compactly supported ones. Some numerical reconstructions are presented at the end of the article. We also briefly discuss the fact that the decomposition presented requires choosing a "central frequency", the choice of which may also yield other decompositions whose reconstruction would be equally fast to implement.
As a second feature of this article, understanding the behavior of data space under this decomposition has also shed some light on the range characterization of the XRT over functions (call it I 0 for the case m = 0) in fan-beam coordinates. More specifically, we prove an equivalence between a range characterization of I 0 given by Pestov and Uhlmann in [32] on simple Riemannian surfaces, and the classical moment conditions due to Gelfand and Graev [10] and Helgason and Ludwig [15,23], translated into fan-beam coordinates. The latter conditions can also be regarded as countably many algebraic conditions characterizing the range of I 0 over functions in the so-called parallel geometry. Translating these conditions into fan-beam coordinates has been and continues to be an object of active study [2,4], in particular for their many applications to the imaging modalities of Computerized tomography and Positron Emission Tomography: motion artifact reduction as in, e.g., [49,51]; consistent extrapolation of truncated data as in, e.g., [48,50]; monitoring of CT systems for faulty detector channels [29]; see also the detailed introduction in [3] on the matter. In the present article, we show that another range characterization provided in [32] is equivalent 2 to the moment conditions for functions with compact support. Additionally, appropriate boundary operators (some of them, labelled as P ± first appearing in [32]; others to be introduced below, found by the author) allow to enforce in a fast way all moment conditions simultaneously. Such operators involve the interplay between the scattering relation and the Hilbert transform on the tangent circles at the boundary of M .
It is the author's hope that such an equivalence bridges a gap between the two systems of coordinates and offers a new viewpoint on the X-ray transform and the moment conditions in fan-beam coordinates, whose theoretical focus in the Euclidean setting is often shifted to parallel coordinates. This is because parallel geometry, unlike fan-beam, enjoys the Central Slice Theorem, which allows for a deeper understanding of discretization issues and regularization theory via accurate and efficient convolution methods in data space [28,9]. This progress in fan-beam coordinates is partly motivated by the fact that dealing with such coordinates becomes almost unavoidable when considering the TTP on Riemannian surfaces without symmetries.
Finally, the expliciteness of the present results owes much to the fact that the scattering relation on a Euclidean disk admits a very explicit expression, thereby simplifying the description of the action of the boundary operators mentioned above. A generalization of these results to simple Riemannian surfaces is currently under study and will appear in future work.
We now discuss the main results and give an outline of the remainder of the article.

Statement of the main results
We consider M = {x = (x, y) ∈ R 2 : x 2 + y 2 < 1} the Euclidean unit disk, with boundary ∂M = {e iβ , β ∈ S 1 }. Here and below, we may identify at will a point x = (x, y) ∈ M with its complex representative z = x + iy. Define the usual Sobolev spaces L 2 (M ), H 1 (M ) and H 1 0 (M ), as well asḢ The domain M has a unit circle bundle SM ∼ = {(x, θ) ∈ M × S 1 } where θ identifies the unit tangent vector cos θ sin θ at x. SM has influx/outflux boundaries 3 where ν x is the unit outer normal at x ∈ ∂M (here we may identify ν x = x). We parameterize the influx boundary in fan-beam coordinates, that is, in terms of a couple (β, α) ∈ S 1 × − π 2 , π 2 , where β describes the point at the boundary x(β) = cos β sin β , and α denotes the direction of the tangent vector with respect to ν x , i.e. v = cos(β+π+α) sin(β+π+α) . We consider the X-ray transform (XRT) I : While this transform is usually made continous in the L 2 (SM ) → L 2 (∂ + SM, cos α) setting, we prove in Lemma 8 that it is in fact L 2 (SM ) → L 2 (∂ + SM ) continuous. It later becomes more convenient to construct Hilbert bases in this unweighted codomain.
Transport equations on SM . The ray transform can be represented as the ingoing trace If = u| ∂ + SM of a solution u(x, θ) to the following transport equation where we have defined the geodesic vector field on T (SM ) by and where ∂ := 1 2 (∂ x − i∂ y ). Using notation from the canonical framing of the unit circle bundle SM , we define The S 1 action on SM induces a Fourier decomposition on L 2 (SM, dx dθ): a function u ∈ L 2 (SM ) can be written uniquely as In this decomposition, a diagonal linear operator of interest is the fiberwise Hilbert transform H : L 2 (SM ) → L 2 (SM ), whose action on each harmonic component is described as We write the decomposition H = H + + H − , where H +/− stands for the composition of H with projection onto even/odd harmonics.
We now explain how definition 1 contains the case where one integrates tensor fields. The correspondence between tensors and functions on SM with finite harmonic content. In two dimensions, any symmetric m-covariant tensor may be represented as with σ denoting the symmetrization operator. For a curve (x(t), θ(t)) in SM , one traditionally defines the ray transform of f as the integral along this curve of f paired m times with the speed vector cos θ(t) sin θ(t) . This amounts to integrating the following function on SM along said curve, in one-to-one correspondence with the tensor f This function can then be uniquely represented in the Fourier decomposition in θ solely in terms of the following harmonics f = f 0 + f ±2 e ±2iθ + · · · + f ±m e ±imθ (m even), f ±1 e ±iθ + f ±3 e ±3iθ + · · · + f ±m e ±imθ (m odd).
For convenience, we still call such elements m-tensors and define S m the subspace of m-tensors in L 2 (SM ). In this identification, we have the following correspondences df ↔ Xf , and for g ∈ H 1 (M ), dg ≡ dg ↔ Xg and dg ↔ X ⊥ g ( : Hodge star operator). In particular, the Helmholtz-Hodge decomposition of one-forms reads as follows: a function w ∈ L 2 (SM ) of the form w = w 1 + w −1 decomposes uniquely in the form Main results. The decomposition problem mentioned in the introduction is decoupled between even and odd tensors. For an integer k, we denote: Theorem 1 (Decomposition). For any m-tensor f ∈ L 2 (SM ), there exists a unique m-tensor g ∈ L 2 (SM ) such that If = Ig, and g is of the following form: (i) if m = 2n is even, then g = g 0 + g ±2 e ±i2θ + · · · + g ±2n e ±i2nθ with g 0 ∈ L 2 (M ) and (g ±k e ±ikθ ) ∈ ker ±k η ∓ for all 2 ≤ k ≤ n, k even.
In both statements above, there exists a constant C m independent of f, g such that As mentioned earlier, this decomposition appears naturally in the context of finding conformal Killing tensor fields on manifolds, see e.g. [41,42] or [5,Theorem 1.5]. An advantage of this representation is that the ray transforms of each component of g live on orthogonal subspaces of L 2 (∂ + SM ), as explained in the following theorem. Theorem 2 (Range decomposition). For any natural integer m, we have the following range decompositions, orthogonal in L 2 (∂ + SM ): Another characterization of the range of the X-ray transform over tensors was provided in [30] on simple surfaces, where the sum there need not be direct as in the decomposition above. In the range decomposition above, the largest subspaces are the ranges of I 0 and I ⊥ , whose description benefits from the Pestov-Uhlmann range characterizations appearing in [32]. More precisely, upon defining operators P ± : C ∞ (∂ + SM ) → C ∞ (∂ + SM ) in terms of the Hilbert transform and the scattering relation, it is proved in [32] that the ranges of I 0 and I ⊥ match those of P − and P + in smooth topologies, respectively (see Section 4.4 for detail). An earlier characterization of the range of I 0 in parallel coordinates, in the form of moment conditions, was found by Gelfand, Graev [10], and separately by Helgason and Ludwig [15,23]. We prove here that, translated into fan-beam coordinates for compactly supported functions, they become equivalent to the Pestov-Uhlmann range characterization. Both theorems 2 and 3 make use of explicit bases of L 2 (∂ + SM ) introduced in Section 4.1. The orthogonality of the ranges in Theorem 2 makes it especially simple to separate the components in data space, and to reconstruct them separately and efficiently. More specifically, the corresponding reconstruction procedure is given below (see Sec. 4.3 for definition of A ± and A ± , and Eq. 34 for the definition of I 0 and I ⊥ ). In the statement below, we denote bẏ Theorem 4 (Reconstruction). Let f ∈ L 2 (SM ) an m-tensor and g as in Theorem 1, and denote the data D = If = Ig.
(i) If m = 2n is even, the functions g 0 , g ±2 , · · · , g ±2n can be uniquely reconstructed by means of the following formulas: (ii) If m = 2n + 1 is odd, the functions g 0 , g ±3 , . . . , g ±(2n+1) can be uniquely reconstructed by means of the formulas: decomposing the function g 0 into g 0 = g (0) + g ( Remark 5. (i) The main ideas in deriving 5 and 8 first appeared in [32], later modified using the operator A + HA − in [26, Proposition 2.2]. Such formulas only inverted I ⊥ for functions vanishing at the boundary. Inversion of I ⊥ over integrands which may be supported up the the boundary in the present paper leads us to introduce the factor Id + (A − HA − ) 2 in 8.
(ii) Each formula above contains a projection onto the appropriate subspace before reconstruction of the corresponding independent component, so that it is not necessary to pre-process the data D before applying these formulas.
(iii) The main bulk of the reconstruction, containing all possible visible singularities, is contained in the g 0 mode, involving the inversion of I 0 or I ⊥ . The additional analytic terms account for amplitude discrepancies in the resulting data, though they do not contain any singularities. These additional terms, unless identically zero, are never compactly supported inside M , even when the initial tensor f is, which was a feature already present in the case of the solenoidal representative discussed in the introduction.
(iv) The decompositions presented above rely on the choice of a particular "central harmonic" g 0 for the reconstructible representative. We discuss in Section 6 the possibility of other decompositions and representatives based on changing the central harmonic, and we briefly discuss how to reconstruct them.
Outline. Section 3 covers the proof of Theorem 1. Section 4 studies the range decomposition of the X-ray transform: we first introduce appropriate bases for L 2 (∂ + SM ) in Sec 3 Decomposition of a tensor -Proof of Theorem 1.
The proof of Theorem 1 relies on viewing the X-ray transform as the influx trace of a solution to certain transport equation posed on SM as in the introduction, and changing both solution and source term without altering the boundary values. We start by stating the following lemma, whose proof is standard and omitted.
(i) There exists a unique couple (v, g) ∈ H 1 0 (M ) × L 2 (M ) such that f = ∂v + g and ∂g = 0, satisfying the stability estimate (ii) By complex conjugation, there also exists a unique couple (v , g ) ∈ H 1 0 (M ) × L 2 (M ) such that f = ∂v + g with ∂g = 0, satisfying the same estimate as above.
Remark 7. In the context of Riemannian surfaces, the actual decomposition to look at is, for every ±k > 0 the splitting of a smooth element f ∈ ker(∂ θ − ikId) into f = η ± v + g with v| ∂SM = 0 and η ∓ g = 0, see [14]. A reader familiar with these decompositions may notice that the proof of Theorem 1 holds for geodesic ray transforms on any Riemannian surface where such decompositions exist.
The proof of Theorem 1 is then based on an iterative use of these decompositions.
Proof of Theorem 1. Proof of (i). We treat a tensor with only non-negative harmonic content and explain how to generalize, i.e., we prove the case of an even tensor of the form f = f 0 + f 2 e i2θ + · · · + f 2n e i2nθ ∈ L 2 (SM ). Recall the transport equation Using lemma 6(i), since f 2n ∈ L 2 (M ), we decompose f 2n = ∂v 2n−1 + g 2n with v 2n−1 | ∂M = 0 and ∂g 2n = 0. We then rewrite this as an equality of functions on SM : so that the transport equation may be rewritten as where ∂g 2n = 0, f 2n−2 − ∂v 2n−1 ∈ L 2 (M ) and all components are controlled by f L 2 (SM ) .
Since v 2n−1 | ∂M = 0, integrating this equation along all straight lines gives the data If , where the right-hand-side now has a top term satisfying ∂g 2n = 0. One may repeat this process for the term of harmonic order 2n − 2 and n − 2 more times, to arrive at a tensor g satisfying Ig = If , of the form where ∂g 2k = 0 for 0 < k ≤ n with an estimate g L 2 (M ) ≤ C f L 2 (M ) , for a constant C independent of f . Now if f has negative harmonic terms, we may use lemma 6(ii) to decompose the negative harmonics in a similar fashion and arrive at the desired result.
Proof of (ii). Similarly to the first case, if f is of the form f = f ±1 e ±iθ +· · ·+f ±(2n+1) e ±i(2n+1)θ , we may use Lemma 6 iteratively to construct a unique with C independent of f and ∂g −(2k+1) = ∂g 2k+1 = 0 for 1 ≤ k ≤ n. We then apply the Helmholtz-Hodge decomposition to the one-form g −1 e −iθ + g 1 e iθ to write it uniquely in the form so that, the transport equation Xu = −f may be recast as where (u + v)| ∂SM = u| ∂SM . The right-hand side is the desired decomposition. The vanishing of v at ∂M implies that the right-hand side of 13 has the same ray transform as f . Theorem 1 is thus proved.

Range decomposition -Proofs of Theorems 2 and 3
For (β, α) ∈ ∂ + SM , denoting ϕ β,α (t) = (e iβ + te i(β+π+α) , β + π + α) the geodesic flow, Santaló's formula reads in our context: This formula usually suggests that the ray transform is naturally defined in the setting I : L 2 (SM ) → L 2 (∂ + SM, cos α), see e.g. [39]. However, when considering the fiberwise Hilbert transform later, the weight cos α may be bothersome in the sense that on L 2 (S 1 , | cos α|), the Hilbert transform cannot be made continuous [34]. We will in fact establish that this weight can be removed so that we can work in a setting where the operator H is well-behaved. Proof. Using Cauchy-Schwarz inequality, we write that Integrating the inequality over ∂ + SM and using Santaló's formula, we obtain that where inequality becomes equality in the case f ≡ 1. Hence the result.
In fan-beam coordinates, L 2 (∂ + SM ) can be regarded as L 2 S 1 × − π 2 , π 2 , for which we define a Hilbert orthonormal basis For any couple (k, l) ∈ Z 2 , we also define φ k,l (β, α) = e iα φ k,l (β, α), and define a second orthonormal basis of L 2 (∂ + SM ) The main reason for introducing two bases is that we will later need to extend functions on ∂ + SM into functions on ∂SM by evenness or oddness w.r.t. the involution α → α + π, and depending on the representation, such processes simply consist of extending the expression at hand to α ∈ S 1 . The change of basis B ↔ B is non-local: indeed, the function α → e iα decomposes into the basis B as follows so that, using the property that φ p,q φ p ,q = φ p+p ,q+q , the formulas for changing basis are given by
Note that S : ∂ ± SM → ∂ ∓ SM and that S A is the composition of S with the antipodal map α → α + π. Both relations are involutions of their respective domains and we consider their associated pull-backs S A f := f • S A and S g := g • S for f ∈ L 2 (∂ + SM ) and g ∈ L 2 (∂SM ). The operators S A : L 2 (∂ + SM ) → L 2 (∂ + SM ) and S : L 2 (∂SM ) → L 2 (∂SM ) are self-adjoint, so when using these operators, we reserve the notation for pull-backs. For basis elements of B, B , it is straightforward to establish the identities, When the expressions φ p,q and φ p,q are regarded as elements on L 2 (∂SM ), we also have the identities Identity 15 makes S A an isometry of L 2 (∂ + SM ). Moreover, we have the following orthogonal splitting Out of the bases B and B , we can then extract two distinct bases for each space V + and V − by successively applying the operators Id ± S A and removing redundancies in (p, q). Given the formulas 15, let us define We obtain two bases for each subspace: V + = span (u p,q , p ≤ 2q) = span (u p,q , p < 2q + 1), All definitions considered are valid for all (p, q), though with the following redundancies: We now summarize how these bases transform under complex conjugation: using the property φ p,q = φ −p,−q true for every (p, q), we have where the pairs of indices in the last column are in the reduced (p, q) ranges as in 18. The next identities summarize how the basis elements transform under Id − S : 4. 3 The operators A ± , A ± , P ± and C ± Define the following operators A ± of even/odd extension with respect to S, i.e. for f smooth on ∂ + SM , Such operators extend continuously to the L 2 (∂ + SM, cos α) → L 2 (∂SM, | cos α|) setting (see e.g. [32]), though here, since the boundary is strictly convex, A ± are also continuous in the L 2 (∂ + SM ) → L 2 (∂SM ) setting. Moreover, because |Jac S| ≡ 1 here, the adjoints in either functional setting are given by In terms of these operators, we now define the operator P := A − HA + , which upon splitting the Hilbert transform into H = H + + H − , yields the decomposition P = P + + P − with P ± = A − H ± A + . A straightforward study of symmetries shows that In matrix notation in the V + ⊕ V − decomposition, one may think of P as P = 0 P − P + 0 . In terms of P ± operators, the following range characterizations were proved in [32] (see also [33]): for w ∈ C ∞ (∂ + SM ), we have • w ∈ Range I 0 if and only if w = P − h for some h ∈ C ∞ (∂ + SM ) with A + h ∈ C ∞ (∂SM ).
• w ∈ Range I ⊥ if and only if w = P + h for some h ∈ C ∞ (∂ + SM ) with A + h ∈ C ∞ (∂SM ).
Such results hold for simple surfaces with boundary, although in that context, it remains open at present how to use this characterization effectively. In the present Euclidean context, we now explain how this characterization can be used.
The bases introduced in Section 4.1 make it easy to compute the singular value decompositions of P ± , which in turn tells us about the harmonic content of the ranges of I 0 and I ⊥ .
(In all cases of nonzero spectral values below, all vectors have norm either 1 or √ 2 though we normalize them with the notation to avoid ambiguity in the spectral values.) Proposition 9. The singular value decomposition for P + : V + → V − is given by: for every (p, q) satisfying p ≤ 2q: if q > 0 and p = q, i v p,q if q = 0 and p < 0, 0 otherwise.
The singular value decomposition for P − : V − → V + is given by: for every (p, q) satisfying p < 2q + 1: Proof. Let E ± : L 2 (∂ + SM ) → L 2 (∂SM ) the operators of even/odd extension with respect to the antipodal map α → α + π. It is important to note for the sequel that (iii) Applying E + to an expression in the basis B consists of extending this expression to α ∈ S 1 .
(iv) Applying E − to an expression in the basis B consists of extending this expression to α ∈ S 1 .
Using the observations (i)-(iv) above, we are then able to derive the following alternative definitions for P ± : These expressions, together with remarks (iii) − (iv) above, suggest that the action of P + is best described in the basis B while that of P − is best described in the basis B . Using identities 20, we then compute Similarly, we have for P − Studying the values of these signs by splitting cases yields the result. Figure 1 locates the ranges of P ± in the (p, q) planes describing the spaces V ± . Of practical interest is the ability to project noisy data onto the range of I 0 . While the operator P − characterizes the range of I 0 , it is not a projection operator and in fact annihilates the range of I 0 (indeed, P 2 − = 0 since P − (V + ) = 0 and P − (V − ) ⊂ V + ). For projection purposes, let us now introduce the operator C : This time, a direct inspection of symmetries shows that In matrix notation in the V + ⊕ V − decomposition, one may think of C as C = C − 0 0 C + . By similar considerations as in the proof of Proposition 9, we can then prove that for any couple (p, q), Splitting cases according to (p, q), we arrive at the following For any (p, q) with p < 2q + 1, i v p,q if q < 0 and p < q, −i v p,q if q > 0 and p > q,

Equivalence of range characterizations of I 0 -Proof of Theorem 3
The classical moment conditions state, in parallel coordinates (see, e.g., [9,28]), that some data R is the Radon transform of some function if and only if R(s, θ) = R(−s, θ + π) and for every integer n ≥ 0, the moment function p n (θ) = R s n R(s, θ) dθ is a homogeneous polynomial of degree n in (cos θ, sin θ). This can be recast as a set of orthogonality conditions in the following form S 1 R s n e ±ikθ R(s, θ) ds dθ = 0, ∀ k > n, k − n even.
We now prove Theorem 3, i.e., that this set of orthogonality conditions is strictly equivalent to proving that the frequency content of D is contained in the range of P − . The proof of Theorem 3 makes use of the following lemma, whose proof (e.g., by induction) is left to the reader: Lemma 11. For every natural integer n ≥ 0, we have Proof of Theorem 3. According to 21, D satisfies the moment conditions if and only if, for every k > 0, D ⊥ span {sin n α cos αe ±ik(α+β) , n = 0, · · · , k − 1, k − n even}.
The previous two calculations describe this set exactly, p-slice by p-slice. Theorem 3 is proved.
On the projection of noisy data onto the range of I 0 . In light of the past two sections, we now propose two approaches in order to project noisy data onto the range of I 0 . Both approaches are extremely fast as they are based solely on the scattering relation and the fiberwise Hilbert transform (which can be computed slice-by-slice using Fast Fourier Transform), and they allow, by virtue of Theorem 3 to enforce the moment conditions at all orders at once.
1. For inversion purposes, the reconstruction formula 5 for functions already contains a projection step. Indeed, applying A + H − A − to the data before backprojection, as prescribed by the formula amounts to applying (up to sign) the adjoint of P − = A − H − A + . By Proposition 9, applying this operator to the data will remove the content which is supported on the orthocomplement of range P − , while moving data from the space V + to V − , a necessary step before applying backprojection I ⊥ .
2. Based on Proposition 10, we have that Id + C 2 − is exactly the projection onto the range of P − , i.e. the range of I 0 .

The operator I ⊥
We have defined I ⊥ :Ḣ 1 (M ) → L 2 (∂ + SM ) while the operator "I 1 " defined in [32] is the ray transform restricted to one-forms. By virtue of the Helmholtz-Hodge decomposition for a square-integrable one-form ω: we have that "I 1 "(ω) = I(Xf ) + I(X ⊥ g) = I ⊥ g, so that the range of these transforms is the same, characterized in frequency content by the operator P + . While the reconstruction formula inverting I ⊥ in [32] only applies to functions in H 1 0 (M ), we now state by how much the spaces H 1 0 (M ) andḢ 1 (M ) differ. Lemma 12. The following direct decomposition holds: Proof of Lemma 12. Let g ∈Ḣ 1 (M ). Then by trace theorems, g| ∂M ∈ H 1 2 (∂M ), and g| ∂M , with zero average at the boundary, can be written uniquely as We now define, inside the domain It is easy to see that g (±) ∈Ḣ 1 (ker 0 η ± , M ), and that one may then decompose g uniquely into The lemma is proved.
We now show that upon applying I ⊥ to g, the three resulting summands live on orthogonal subspaces of V − , each of which corresponding to a different nonzero spectral value of P + , and to the 0, i 2 and −i 2 -eigenspaces of C + . Proposition 13. The mapping I ⊥ :Ḣ 1 (M ) → V − maps the direct decomposition 24 into three orthogonal subspaces of V − , corresponding to the three distinct nonzero spectral values of P + (as described in Proposition 9). Moreover, Proof of Proposition 13. Write g ∈Ḣ 1 (M ) as g = g (0) + g (+) + g (−) according to 24. We first analyze the summands g (+) and g (−) . On to the summand g (−) , we have 0 = η − g (−) = (X − iX ⊥ )g (−) , so that, in particular, Using the fundamental theorem of calculus along each line of integration, this equals so I ⊥ g (−) is in the subspace of V − corresponding to the spectral value −i of P + , and according to Proposition 10, C + I ⊥ g (−) = −i 2 I ⊥ g (−) . Similary for g (+) , we have 0 = η + g (+) = (X + iX ⊥ )g (+) , so that, in particular, X ⊥ g (+) = iXg (+) , thus I ⊥ g (+) (β, α) = iI[Xg (+) ](β, α), which yields so I ⊥ g (+) is in the subspace of V − corresponding to the spectral value i of P + , and according to Proposition 10, C + I ⊥ g (+) = i 2 I ⊥ g (+) . We conclude by showing that the range of I ⊥ restricted to H 1 0 (M ) is orthogonal to both these subspaces. This fact mainly derives from the following lemma, whose proof we relegate at the end of this proof.
Lemma 14. For any f ∈ H 1 0 (M ), then I ⊥ f satisfies: Assuming Lemma 14 is proved, it is easy to see that, if f ∈ H 1 0 (M ), and using the fact that for any k = 1, 2, . . . , v k,k = −(Id − S A )φ k,0 and S A I ⊥ f = −I ⊥ f , we arrive at, Similarly using that v −k,0 = (Id − S A )φ −k,0 , we arrive at Since the harmonic content of the ranges of P + and I ⊥ must match, the range of I ⊥ restricted to H 1 0 (M ) cannot but be spanned by the space corresponding to the spectral value −2i of P + , on which C + vanishes identically. Proposition 13 is proved.
Proof of Lemma 14. It is enough to prove this for f ∈ C ∞ c (M ) and the result follows by density. We compute upon using the chain rule and the fact that f | ∂M = 0. Now since f has compact support inside M , for α close enough to ± π 2 , f is identically zero along the segment {e iβ + te i(β+α+π) , t ∈ (0, 2 cos α)}, so the last right-hand-side above is zero. Lemma 14 is proved.

4.6
The ranges over ker m η ± Let us first show how the characterization of the range of I 0 helps us understand the lack of injectivity over tensors. Fix m ∈ Z and f ∈ L 2 (M ) and suppose we want to compute I[f (x)e imθ ] = I m f . Then we obtain: This implies that, for instance if m = 2, the range of I 2 amounts to translating the range of I 0 in the (p, q)-plane by (2, 1), so that the overlap between the ranges of I 0 and I 2 is rather large. The next observation is that, in fact, the only part of the range of I 2 which does not lie in that of I 0 is precisely I(ker 2 η − ) !
We now formulate the range characterizations of I(ker m η ± ). To this end, if {u k } ∞ k=0 is a sequence of orthogonal vectors in a Hilbert space such that there exists C ≥ 1 with C −1 ≤ u k ≤ C for every k, we define: Proposition 15. For any m ∈ Z, we have Proof of Proposition 15. Study of ker 0 η ± . Suppose f ∈ L 2 (M ) is a holomorphic function, i.e. f ∈ ker 0 η − = ker ∂. Then we may write f as f (z) = ∞ k=0 a k z k for some numbers {a k } k such that Then for any (β, α) ∈ ∂ + SM , This implies straightforwardly that If f ∈ L 2 (M ) is antiholomorphic, of the form f (z) = ∞ k=0 a k z k , then f is holomorphic and we can compute, using 33, where we have used 19 in the last equality. We therefore deduce that Study of ker m η ± for m = 0. From 27, recall that and f (x)e imθ ∈ ker m η ± if and only if f ∈ ker 0 η ± so that multiplication by e im(β+α) is an isomorphism between I(ker 0 η ± ) and I(ker m η ± ). In terms of our bases B and B , this is framed as follows: Even case: we have that I[f (x)e 2imθ ] = φ 2m,m I 0 f , and using the fact that φ 2m,m u p,q = u p+2m,q+m for any p, q, m, we deduce 29 and 30.
The proof of Proposition 15 is complete.

Reconstruction formulas -Proof of Theorem 4
In light of Theorems 1 and 2, every m-tensor f has an equivalent m-tensor g such that If = Ig and such that the ray transforms of each component of g live on orthogonal subspaces of L 2 (∂ + SM ). We now explain how to reconstruct each component. Section 5.1 covers the inversion of I 0 and I ⊥ , while Section 5.2 covers the reconstruction of elements in any space ker m η ± for m ∈ Z.

Inversion of I 0 and I ⊥
Inversion formulas for I 0 were long known in the parallel geometry (see [35,28]), and reconstruction formulas in fan-beam geometry were obtained by changing variable in the former inversion [16]. Seen in a context of transport equations on Riemannian surfaces, the first inversion formulas for the operators I 0 and I ⊥ (restricted to smooth functions with compact support) appeared in [32], though not being represented using the A + HA − operator as we presently do. As explained earlier, adding this operator allows to project data onto the ranges of I 0 and I ⊥ . For a proof of the inversion formulas justifying the use of the operator A + HA − , see [26,Proposition 2.2] (note that the error operators W, W appearing there account for curvature and vanish identically in the present Euclidean setting). Above, we have defined the backprojection operators as follows: for D(β, α) defined on ∂ + SM and x ∈ M , Remark 16. Although we favored the unweighted space L 2 (∂ + SM ) for reasons of range description here, it is to be noted that when deriving reconstruction formulas, the operators I 0 and I ⊥ appear naturally and they are the adjoints of I 0 and I ⊥ when considering the space L 2 (∂ + SM, cos α) as their codomain, instead of the unweighted one.
Inversion of I ⊥ overḢ 1 (M ). As the previous section only covered the case of function vanishing at the boundary, we must now consider the inversion over the decomposition g (0) + g (+) + g (−) as in Section 4.5.
It turns out that applying the inversion formula above to I ⊥ (g (0) + g (+) + g (−) ) does not only pick up g (0) but also a fraction of g (+) and g (−) . In order to remove this coupling, we derive independent reconstruction formulas for g (±) exploiting their analyticity, and we introduce an additional boundary operator which will remove I ⊥ (g (+) + g − ) from I ⊥ g to reconstruct g (0) separately.
Proof of Equations 9 and 10. Call D = I ⊥ (g (0) +g (+) +g (−) ) and let us reconstruct the functions g (±) from D. Using the calculation 25, we have If the data is consistent in the sense that Id−S A 2 D = D, then we can simplify this formula into 9, hence the formula.
On to the function g (+) , using the fact that I ⊥ g (+) = I ⊥ g (+) and the fact that g (+) is holomorphic, we can use the previous formula to arrive at 10.
Calling D = I ⊥ (g (0) + g (+) + g (−) ), all we need to prove is that (Id + (A − HA − ) 2 )D = I ⊥ g (0) , or in other words, that where C + is defined in Section 4.3. This is an immediate consequence of Proposition 13.
Combining this with 35, we arrive at which is what we had to prove.
If h ∈ ker m η + , the proof of the reconstruction formula is a similar combination of 27 and 36.

On the possibility of other decompositions
Let us restrict the present discussion to the case of even tensors, though the case of odd tensors is similar. The decomposition described in Theorem 1 is based on a particular choice of "main harmonic" g 0 . However, for a 2n-tensor f with data If , one could very well pick any integer −n ≤ k ≤ n and construct a 2n-tensor g with If = Ig such that, • For −n ≤ < k, g 2 e i2 θ ∈ ker 2 η + , reconstructible from If using the formulas derived in the present article.
• For k < ≤ n, g 2 e i2 θ ∈ ker 2 η − , reconstructible from If using the formulas derived in the present article.
• The "main harmonic" is of the form e 2ikθ g 2k (x) with g 2k ∈ L 2 (M ), and is to be reconstructed from the transform I[e 2ikθ g 2k ], for which the author has given inversion formulas in [25,Theorem 5.2]. These formulas are similar in spirit to the case k = 0 except that they involve conjugations of the Hilbert transform by factors e ±i2kθ .
Such decompositions are also continous in the sense that g L 2 (SM ) ≤ C f L 2 (SM ) and enjoy the same efficiency in implementation. The main difference with the case k = 0 is that the harmonic g 2k would now be the one which contains all visible singularities. Figure 2 illustrates, on a 2-tensor reconstruction problem, two different ways in which to view the frequency content of the transform of a 2-tensor, to reconstruct two different candidates.

Numerical experiments
We now show some numerical illustrations of the reconstruction formulas, one for a tensor of even order, one for a tensor of odd order, using the Matlab code previously documented by the author in [24] in the context of Riemannian surfaces. The underlying grid is 300 × 300 cartesian, the discretization of the data space S 1 × (− π 2 , π 2 ) is 600 × 300 equispaced. All computations terminate within seconds on a regular personal computer.

Experiment 1: Reconstruction of a second order tensor
We first illustrate Theorem 4 with an example of a second-order tensor f = f 0 +f 2 e 2iθ +f −2 e −2iθ . For simplicity of display, we make f real-valued by imposing the constraints f 0 = f 0 and f 2 = f −2 . In this case, if we write f 2 = f r 2 + if i 2 , the tensor f takes the form to two different candidates g = g −2 e −i2θ + g 0 + g 2 e i2θ .
In tensor notation, this corresponds to the symmetric 2-tensor The functions f 0 , f r 2 , f i 2 are given in Fig. 3. As prescribed in Theorem 4, we reconstruct the equivalent tensor g = g 0 +g 2 +g −2 according to the formulas there. One should note in particular that since the data is real-valued, the reconstructed g satisfies g 0 = g 0 and g 2 = g −2 , so that we may represent g just like f , that is, in terms of functions g 0 , g r 2 and g i 2 as in 37. These three functions are represented in Fig. 4.
The forward data If , as well as the pointwise difference |If − Ig| are given on Fig. 5. Some comments are in order: 1. Even if the tensor has compact support (which is the case here), the reconstructed tensor may not have compact support.
2. The inability to separate the singularities of f 0 , f 2 , f −2 in the reconstruction is not a weakness of the method: the lack of injectivity implies the loss of such an information.  As seen in section 6, one could choose to reconstruct another h = h −2 + h 0 + h 2 with h −2 ∈ ker −2 η + , h 0 ∈ ker 0 η + and h 2 ∈ containing all reconstructed singularities, and this would make another perfectly acceptable candidate.

Experiment 2: Reconstruction of a solenoidal vector field
In the case of tensors of odd order, there is an additional numerical technicality: as the Cauchy integral type of inversion formulas may become unstable too close to the boundary, we need to cut off the reconstruction near the boundary. This in turn would create artifacts when recomputing I ⊥ of the reconstructed function and comparing it to the initial data, which makes this comparison method unreliable.
As an alternative, we first show that the inversion procedure of I ⊥ overḢ 1 (M ) is successful almost up to the boundary by comparing directly the reconstruction with the initial function (we can do that because this problem is injective). Once this is done, the next section will give an example of a reconstruction on an example of a third-order tensor.
In this example we take f = f (0) +f ∂ as the sum of a compactly supported function f (0) (sum of peaked gaussians) and a non-compactly supported harmonic term f ∂ = (z 3 ), as in Figure  6 (left). The forward data I ⊥ f is on the right of Figure 6. We first apply Id + (A − HA − ) 2 to the data to extract I ⊥ f (0) (Fig. 7, left) and reconstruct f (0) (see Fig. 8, left) from it. For convenience, I ⊥ f ∂ , extracted from the data, is visualized Fig. 7 (right). We then apply formulas 9-10 to the data to reconstruct f ∂ (see Fig. 8, middle). The pointwise error inside the centered disk of radius 0.99 is given on Fig. 8 (right).

Experiment 3: Reconstruction of a third-order tensor
We now give an example of a third-order tensor f = f ±1 e ±iθ + f ±3 e ±i3θ , where again, to simplify the exposition, we assume f real-valued so that f −3 = f 3 and f −1 = f 1 . If we write f k = f r k + if i k for k = 1, 3, we obtain that the tensor f takes the form   which in tensor notation corresponds to the tensor An example of such functions is given Figure 9, with the ray transform of the corresponding 3-tensor given in Figure 10. By Theorem 1, we reconstruct an equivalent real-valued 3-tensor of the form g = X ⊥ g 0 + 2g r 3 cos(3θ) − 2g i 3 sin(3θ), g 0 ∈Ḣ 1 (M ), (g r 3 ± ig i 3 )e ±i3θ ∈ ker ±3 η ∓ , according to Theorem 4. The reconstructed functions are given in Figure 11. Note that in this case, the function g 0 contains all visible singularities, though all of them smoother by 1 degree than the ones of the initial tensor f . This is because the contribution of g 0 in the equivalent tensor is in fact X ⊥ g 0 and not g 0 itself.   3: Left to right: the functions g 0 , g r 3 , g i 3 of the reconstructed tensor g = X ⊥ g 0 + 2g r 3 cos(3θ) − 2g i 3 sin(3θ), whose ray transform equals If .