Inversion of the attenuated geodesic X-ray transform over functions and vector fields on simple surfaces
IInversion of the attenuated geodesic X-ray transform overfunctions and vector fields on simple surfaces
Fran¸cois Monard ∗ November 13, 2018
Abstract
We derive explicit reconstruction formulas for the attenuated geodesic X-ray transformover functions and, in the case of non-vanishing attenuation, vector fields, on a class ofsimple Riemannian surfaces with boundary. These formulas partly rely on new explicitapproaches to construct continuous right-inverses for backprojection operators (and, in turn,holomorphic integrating factors), which were previously unavailable in a systematic form.The reconstruction of functions is presented in two ways, the latter one being motivatedby numerical considerations and successfully implemented at the end. Constructing theright-inverses mentioned require that certain Fredholm equations, first appearing in [29], beinvertible. Whether this last condition reduces the applicability of the overall approach to astrict subset of simple surfaces remains open at present.
We consider explicit inversion formulas for the two-dimensional attenuated X-ray (or, equiv-alently in two dimensions, Radon) transform of a function or vector field. Let (
M, g ) a non-trapping Riemannian surface-with-boundary with unit circle bundle SM = { (x , v ) ∈ T M, g ( v, v ) = 1 } , and let us denote the geodesic flow ϕ t : SM → SM by ϕ t (x , v ) = ( γ x ,v ( t ) , ˙ γ x ,v ( t )), where γ x ,v ( t )denotes the basepoint and ˙ γ x ,v ( t ) is the (unit) velocity vector. For x ∈ ∂M , let ν x the unit innernormal and define the influx/outflux boundaries ∂ ± SM = { (x , v ) ∈ ∂SM, ± g ( ν x , v ) > } . Fix a a smooth enough function on M . Then for a function f ∈ L ( SM ), we define the attenuated(geodesic) X-ray transform of f by I a f (x , v ) := (cid:90) τ (x ,v )0 f ( ϕ t (x , v )) exp (cid:18)(cid:90) t a ( γ x ,v ( s )) ds (cid:19) dt, (x , v ) ∈ ∂ + SM, ∗ Department of Mathematics, University of Michigan, Ann Arbor MI, 48109; [email protected] a r X i v : . [ m a t h . A P ] D ec here τ (x , v ) denotes the first exit time of the geodesic γ x ,v ( t ) (the non-trapping conditionimplies that τ is uniformly bounded above on SM ). Such a transform can be regarded as theinflux trace I a f = u | ∂ + SM of the solution u to the transport equation on SMXu + au = − f ( SM ) , u | ∂ − SM = 0 , (1)and where X = ddt ϕ t | t =0 is the geodesic vector field. Cases of interest here are ( i ) the case where f is a function on M i.e. f (x , v ) = f (x) with applications to X-ray Tomography in media withvariable refractive index, a topic receiving much interest at the moment [22, 15], and ( ii ) thecase where f represents a vector field F on M via the relation f (x , v ) = g ( F (x) , v ), and whoseapplications to Doppler tomography justify its nickname of Doppler transform .Such a transform over functions was studied earlier in the Euclidean setting. Inversion meth-ods with known attenuation were obtained independently by Arbuzov, Bukhgeim and Kazantzevusing A -analytic function theory `a la Bukhgeim in 1998 [1], and by Novikov via complexifica-tion methods in 2002 [25, 23], see [4] for a joint study of both approaches. While both methodspresent two interesting ways of looking at this problem, the latter is of lesser complexity (i.e.,comparable to that of an inverse Radon transform) and is therefore more amenable to imple-mentation, as illustrated in [21, 7]. This latter method was extended to partial data and moregeneral integrands in [2], to the attenuated Radon transform in hyperbolic geometry and to thehorocyclic transform in [3], to even more general curves in [10], and more recently, to weightedRadon transforms in the Euclidean plane, where the weight has finite harmonic content in [24],with an implementation in [8]. Both methods were also adapted to the study of the attenuatedtransform over functions and vector fields in fan-beam geometry in [12], describing both a so-lution via A -analytic functions and an approach using fiberwise holomorphic solutions to thetransport equation (1) in the Euclidean setting. This last approach was then generalized by Saloand Uhlmann to simple Riemannian surfaces (i.e. surfaces with strictly convex boundary andno conjugate points) in [34], where a method reconstructing functions from knowledge of theirgeodesic X-Ray transform was developed. The Doppler transform was studied microlocally inthe Riemannian case in [11], and a range characterization in the Euclidean case was recentlygiven in [32]. Additional range characterizations of the Euclidean transform in convex domainsof R were provided in [33], and a study of the attenuated Euclidean transform over second-ordertensors was recently provided in [31].The X-ray transform can also be considered over vector-valued unknowns defined on SM (i.e.sections of bundles), with matrix-valued attenuations, connections and Higgs fields, for whichrecent results can be found in [27, 26, 6]. More general settings include the case of weighted X-ray transforms (of which the attenuated case is a particular example). In this setting, earlyEuclidean implementations have been provided in [14], and analytic microlocal analysis hasled in [5] to injectivity and stability of the attenuated transform over functions when both themetric and the attenuation coefficient are real-analytic. In dimension d ≥
3, local injectivityof weighted X-ray transforms near convex boundary points was recently established in [37],following methods in [36] where the unattenuated case was first treated.2hile laying the groundwork of inversions on surfaces with non-constant curvature in [34] byproving injectivity and providing the first reconstruction procedure, the authors there pointedout two open questions: ( i ) how to explicitely invert the unattenuated X-ray transform overfunctions and solenoidal vector fields (we call them I and I ⊥ here), and ( ii ) how to explicitelyconstruct (fiber-)holomorphic integrating factors when curvature is not constant. Here by holo-morphic integrating factor (HIF) for the attenuation a , we mean a function which is holomorphicon the fibers (see Sec. 2.1), and which turns the attenuated transport equation (1) into an unat-tenuated one. As will be seen below, after performing Fourier analysis on the fibers of SM , atransport equation such as (1) takes the form of a tridiagonal, doubly infinite system of partialdifferential equations on M , which is best understood when it can be made one-sided. As fiber-holomorphic functions have one-sided harmonic content, HIFs allow to fulfill this requirement toa certain extent. In all past aproaches analyzing attenuated transforms, a “holomorphization”step of integrating factors always occurs, be it on the fibers of SM when solving a transportequation [12, 34], in terms of a complexified parameter when formulating a Riemann-Hilbertproblem [25, 3, 2], or when considering sequence-valued systems [1, 31, 32, 33].In an attempt to address the open questions mentioned above, an answer to ( i ) was presentedby the author in [18], by implementing the reconstruction formulas proposed in [29, 13, 28]. Afirst feature of the present article is to address point ( ii ) by providing explicit ways of construct-ing holomorphic integrating factors, obtained by first constructing preimages of the adjointoperators I (cid:63) and I (cid:63) ⊥ explicitely. Second, we derive other reconstruction algorithms, some ofwhich are similar in spirit to those in [12], another one similar to the approach in [29, 19]. Inthe first approach, the main novelty here, partly motivated by the approach in [12], consists indecomposing v = ue − w (with u the solution of (1) and w a holomorphic solution of Xw = − a )into the sum of a holomorphic function and a function that is constant along geodesics. This canbe done as soon as one can construct explicit preimages of the operators I (cid:63) and I (cid:63) ⊥ . A secondapproach reconstructing functions, more similar to approaches in [29, 19], is then presented, anda numerical implementation is provided. Additionally, these reconstruction formulas are fastin that no full three-dimensional transport equation needs to be solved unlike [34]. The classof surfaces where the current approach is valid is that of simple surfaces where some Fredholmequations (see Equations (10) and (11) below), which first appeared in [29, Theorem 5.4], areinvertible. It remains open at present whether this latter additional requirement applies to allor a strict subset of simple surfaces, though past numerical experiments done in [18] by theauthor have showed that such Fredholm equations were invertible on some family of surfaceswhich could become arbitrarily close to non-simple.Recent work by the author with P. Stefanov and G. Uhlmann [20] shows that stable inversionof the attenuated ray transform should still be possible in some surfaces where conjugate pointsoccur at most in pairs. Adapting the current approach to this latter setting will be the objectof future work. 3 utline. The structure of the paper is as follows. We first introduce in Sec. 2 the basicsetting (Sec. 2.1), some new notation and operators (e.g., I ⊥ ) which play an important rolein the inversion process (Sec. 2.2), the construction of explicit, continuous right-inverses for I (cid:63) and I (cid:63) ⊥ (Sec. 2.3) and how they allow the explicit construction of holomorphic integratingfactors (Sec. 2.4). Section 3 presents the reconstruction formulas for functions and vector fields,following initial ideas in [12], first adapted in [34]. Section 4 proposes an alternative approachto reconstruction of functions, stating in passing L → L continuity and bounds for a certainfamily of operators, the proof of which is relegated to appendix A. Section 5 presents a numericalimplementation of inversions as set up in Theorem 4.2. We briefly recall the geometry and notation associated with the unit circle bundle. With (
M, g )as above, the geodesic flow ϕ t : SM → SM is defined on the domain D = { (x , v, t ) : (x , v ) ∈ SM, t ∈ ( − τ (x , − v ) , τ (x , v )) } . (2)with X the geodesic vector field as in the introduction, one may construct a global frame { X, X ⊥ , V } of T ( SM ), which encodes the geometry of M via the structure equations[ X, V ] = X ⊥ , [ X ⊥ , V ] = X, [ X, X ⊥ ] = − κV ( κ : Gaussian curvature) .V is sometimes referred to as the vertical derivative and X ⊥ the horizontal derivative . Forx ∈ M , we also denote S x := { v ∈ T x M, g ( v, v ) = 1 } the unit tangent circle at x.Though the proofs will be coordinate-free, a convenient set of coordinates is that of isothermalcoordinates (they can be made global as the simplicity assumption on M implies that it is simplyconnected), for which the metric g is scalar of the form g = e λ ( x,y ) ( dx + dy ), and where(x , v ) ∈ SM is parameterized by ( x, y, θ ) where x = ( x, y ) and v = e − λ ( x,y ) (cos θ∂ x + sin θ∂ y ) for θ ∈ S . In these coordinates, the frame { X, X ⊥ , V } takes the expression X = e − λ (cos θ∂ x + sin θ∂ y + ( − sin θ∂ x λ + cos θ∂ y λ ) ∂ θ ) , V = ∂ θ ,X ⊥ = − e − λ ( − sin θ∂ x + cos θ∂ y − (cos θ∂ x λ + sin θ∂ y λ ) ∂ θ ) . We use the Sasaki metric on SM for which the basis ( X, X ⊥ , V ) is orthonormal, with volumeform d Σ preserved by the frame (in isothermal coordinates, d Σ = e λ dx dy dθ ). Introducingthe inner product ( u, v ) = (cid:90) SM u ¯ v d Σ , u, v : SM → C , L ( SM, C ) decomposes orthogonally as a direct sum L ( SM, C ) = (cid:77) k ∈ Z H k , where H k := ker( V − ikId ) . (3)We also denote Ω k := C ∞ ( SM ) ∩ H k , and denote the corresponding decomposition u (x , v ) = (cid:80) k ∈ Z u k (x , v ) with u k ∈ H k . In isothermal coordinate, this corresponds to Fourier series expan-sions u (x , θ ) = ∞ (cid:88) k = −∞ u k (x , θ ) , where u k (x , θ ) = e ikθ ˜ u k (x) , ˜ u k (x) = 12 π (cid:90) S u (x , θ ) e − ikθ dθ. In particular, we denote either by u or π u the fiberwise average u (x) = π (cid:82) S x u (x , v ) dS ( v ),so π : L ( SM ) → L ( M ) is the projection onto H . Such functions admit an even/odddecomposition w.r.t. to the involution v (cid:55)→ − v (or θ (cid:55)→ θ + π ), denoted u = u + + u − , where u + := (cid:88) k even u k and u − := (cid:88) k odd u k . (4)An important decomposition of X and X ⊥ due to Guillemin and Kazhdan (see [9]) is given bydefining η ± := X ± iX ⊥ , so that one has the following decomposition X = η + + η − and X ⊥ = 1 i ( η + − η − ) , with the important property that η ± : Ω k → Ω k ± for any k ∈ Z , so that both X and X ⊥ mapodd functions on SM into even ones and vice-versa.In the harmonic decomposition above, a diagonal operator of particular interest is the so-called fiberwise Hilbert transform H : C ∞ ( SM ) → C ∞ ( SM ), whose action on each componentis described by Hu k := − i sign( k ) u k , k ∈ Z , with the convention sign(0) = 0 , and we denote H + / − the composition of H with projection onto even/odd Fourier modes. Wesay that a function u ∈ L ( SM ) is (fiber-)holomorphic if ( Id + iH ) u = u , i.e. if u has onlynonnegative Fourier components. An important identity first proved in [30] is the commutatorbetween the Hilbert transform and the geodesic vector field:[ H, X ] = π X ⊥ + X ⊥ π , [ H, X ⊥ ] = − π X − Xπ . Note also that H = − Id + π and that Hπ = π H = 0. Using these observations and thecommutators above, we write[ H , X ] = H [ H, X ] + [
H, X ] H = (cid:24)(cid:24)(cid:24)(cid:24) Hπ X ⊥ + HX ⊥ π + (cid:24)(cid:24)(cid:24)(cid:24) X ⊥ π H + π X ⊥ H On the other hand, [ H , X ] = [ − Id + π , X ] = π X − Xπ . Upon splitting into odd and evenparts, we arrive at the following equalities, to be used subsequently: π X = π X ⊥ H and Xπ = − HX ⊥ π . (5)5 .2 An alternative notation for the unattenuated case We now introduce some notation that emphasizes the L ( M )-duality arising in the Pestov-Uhlmann reconstruction formulas [29]. The main novelty below is the introduction of I ⊥ . Thegeneral unattenuated transform can be defined over functions f ∈ L ( SM ) as If (x , v ) = (cid:90) τ (x ,v )0 f ( ϕ t (x , v )) dt, (x , v ) ∈ ∂ + SM. (6)Considering integrands of the form f ∈ L ( M ) (i.e. f (x , v ) = f (x)), and X ⊥ h for some h ∈ H ( M ), let us define the unattenuated transforms I f := If, I ⊥ h := I ( X ⊥ h ) , (that is, in the definition of I , we identify f with its pullback by the canonical projection SM → M ). These transforms are continuous in the following settings when ( M, g ) is non-trapping I : L ( M ) → L µ ( ∂ + SM ) , I ⊥ : H ( M ) → L µ ( ∂ + SM ) , where L µ is a weighted L space with weight µ (x , v ) = g ( ν x , v ). Define the scattering relation α : ∂SM → ∂SM as follows: if (x , v ) ∈ ∂ + SM , α (x , v ) = ϕ τ (x ,v ) (x , v ) ∈ ∂ − SM , and if(x , v ) ∈ ∂ − SM , α (x , v ) = ϕ − τ (x , − v ) (x , v ) ∈ ∂ + SM . Recall the definitions of A ± : L µ ( ∂ + SM ) → L | µ | ( ∂SM ) and their adjoints A (cid:63) ± , introduced in [29]: A ± w ( x, θ ) := (cid:26) w ( x, v ) on ∂ + SM, ± w ◦ α ( x, v ) on ∂ − SM, A (cid:63) ± u := ( u ± u ◦ α ) | ∂ + SM . The fundamental theorem of calculus along a geodesic reads
IXu = ( u ◦ α − u ) | ∂ + SM = − A (cid:63) − u. (7)With ψ : SM → ∂ − SM the endpoint map ψ (x , v ) = ϕ τ (x ,v ) (x , v ) and h ∈ L µ ( ∂ + SM ), we denoteby h ψ := h ◦ α ◦ ψ the function extended by free geodesic transport to SM , i.e. solution of theequation Xu = 0 ( SM ) , u | ∂ + SM = h. Straightforward computations using Santal´o’s formula yield that for h ∈ L µ ( ∂ + SM ), I (cid:63) h = 2 π ( h ψ ) , I (cid:63) ⊥ h = − π ( X ⊥ h ψ ) . (8)6 he V ± decomposition. Let us define the antipodal scattering relation α : ∂ + SM → ∂ + SM as ( v (cid:55)→ − v ) ◦ α , and write L µ ( ∂ + SM ) as the direct sum V + ⊕ V − , where h ∈ V + (resp. V − )iff h is even (resp. odd) with respect to the involution α . Since a function of x only can beregarded as an even function of v on SM , and a vector field can be regarded as an odd functionof v on SM , it is straighforward to establish thatRange I ⊂ V + and Range I ⊥ ⊂ V − . Moreover, we have the following lemma.
Lemma 2.1.
The direct sum L µ ( ∂ + SM ) = V + ⊕ V − is orthogonal.Proof. For h ∈ L µ ( ∂ + SM ), define u = (cid:16) h √ τ (cid:17) ψ where τ = I ( ) ( denotes the constant functionequal to 1 on M ). Santal´o’s formula allows to show that the map h (cid:55)→ u is L µ ( ∂ + SM ) → L ( SM )continuous and (cid:107) u (cid:107) L ( SM ) = (cid:107) h (cid:107) L µ ( ∂ + SM ) . Moreover, u is even/odd in v whenever h ∈ V + / V − ,so that, if h ∈ V + and g ∈ V − , and upon calling u = (cid:16) h √ τ (cid:17) ψ and v = (cid:16) g √ τ (cid:17) ψ , we have( h, g ) L µ ( ∂ + SM ) = ( u, v ) L ( SM ) = 0 , hence the proof. Inversion of I and I ⊥ . We now revisit the inversion of the operators I and I ⊥ , previouslyestablished in [29], adapted here to the present notation. Recall the notation u f ( f ∈ L ( SM ))for the solution of a transport problem of the form Xu = − f ( SM ) , u | ∂ − SM = 0 , (9)and for f ∈ C ∞ ( M ), define W f = ( X ⊥ u f ) . It is established in [29] that W extends as asmoothing (hence compact) operator W : L ( M ) → C ∞ ( M ) and that the L -adjoint of W isgiven by W (cid:63) h = ( u X ⊥ h ) . Proposition 2.2.
Let ( M, g ) a simple surface with boundary. Then we have for every f ∈ L ( M ) and every h ∈ H ( M ) , f + W f = 12 π I (cid:63) ⊥ w, w = 14 A (cid:63) + H − A − I f, (10) h + ( W (cid:63) ) h = − π I (cid:63) w, w = 14 A (cid:63) + H + A − I ⊥ h. (11) Remark 2.3.
Formulas (10) and (11) differ slightly from [29, Theorem 5.4] because it is statedthere that the solution to the transport problem Xu = − f, u | ∂ + SM = w, is u f + w ψ , which is what the Fredholm equation there is based upon. The correct answer wouldbe u f + ( w − I f ) ψ , which in turn yields the modified formula (10) . roof. Proof of (10) . Let f ∈ C ∞ ( M ) and define u f as in (9) so that u | ∂ + SM = I f . Applying π X = π X ⊥ H (derived in (5)) to (9), we obtain f = π f = − π Xu f = − π X ⊥ Hu f = − π X ⊥ Hu f − = − ( X ⊥ Hu f − ) , so f can be obtained from the last equation if we can relate Hu f − to the known data I f . Inorder to do so, we use the commutator [ H, X ] to write a transport equation for Hu f − Xu f − = − f ⇒ X ( Hu f − ) = − (cid:24)(cid:24)(cid:24)(cid:24)(cid:24) X ⊥ ( u f − ) − ( X ⊥ u f − ) = − W f, so that Hu f − satisfies the transport problem X ( Hu f − ) = − W f, Hu f − | ∂ − SM = 12 H − A − I f | ∂ − SM := η, which means that Hu f − = u W f + w ψ , w = η ◦ α. Upon applying ( X ⊥ · ) , we obtain( X ⊥ Hu f − ) = W f − π I (cid:63) ⊥ w. Since we have established that f = − ( X ⊥ Hu f − ) , we conclude that f + W f = 12 π I (cid:63) ⊥ w, w = (cid:18)
12 ( H − A − I f ) | ∂ − SM (cid:19) ◦ α. In terms of A (cid:63) ± operators, we can also write w as w = A (cid:63) + − A (cid:63) − H − A − I f . Moreover, it canbe seen that the function A (cid:63) − H − A − I f has V + symmetry, so it is annihilated by I (cid:63) ⊥ . Thus,Equation (10) follows, extended to every f ∈ L ( M ) by density of C ∞ ( M ) in L ( M ). Proof of (11) . Let h ∈ C ∞ ( M ) with h | ∂M = 0, and let u X ⊥ h solve the transport equation Xu = − X ⊥ h ( SM ) , u | ∂ − SM = 0 . Upon projecting onto odd functions of v , we have Xu X ⊥ h + = − X ⊥ h . Direct manipulations andthe use of the commutator formula imply Xu X ⊥ h + = − X ⊥ h = ⇒ X ( Hu X ⊥ h + − h ) = − X ⊥ W (cid:63) h − (cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24) ( X ⊥ u X ⊥ h + ) = − X ⊥ W (cid:63) h. Moreover, the trace of u X ⊥ h + is given by u X ⊥ h + | ∂SM = A − I ⊥ h . Since the function Hu X ⊥ h + satisfies the transport problem X ( Hu X ⊥ h + − h ) = − X ⊥ W (cid:63) h, Hu X ⊥ h + | ∂ − SM = 12 H + A − I ⊥ h | ∂ − SM ,
8e deduce that Hu X ⊥ h + − h = u X ⊥ W (cid:63) h + w ψ , where w := (cid:18) H + A − I ⊥ h | ∂ − SM (cid:19) ◦ α. Upon averaging over θ and rearranging terms, we obtain h + ( W (cid:63) ) h = − π I (cid:63) w. In terms of A (cid:63) ± operators, w can also be written as w = A (cid:63) + − A (cid:63) − H + A − I ⊥ h . We can see that A (cid:63) − H + A − I ⊥ f has V − symmetry and as such is annihilated by I (cid:63) . Thus, Equation (11) follows,extended to every h ∈ H ( M ) by density. Remark 2.4.
Inspection of symmetries shows that A − I f in (10) is odd in v and A − I ⊥ f in (11) is even in v . This tells us that the expressions H − A − I f and H + A − I ⊥ f are redundant, asone could just replace H ± by H . This further emphasizes the similarity between formulas (10) and (11) , for both of which the operator A (cid:63) + HA − acts as a first step in the postprocessing ofdata, though on a different subspace depending on the formula. I (cid:63) and I (cid:63) ⊥ The question of surjectivity of I (cid:63) and I (cid:63) ⊥ have proved to be crucial for answering boundary rigid-ity questions (see [30] where a proof of surjectivity of I (cid:63) appears) and constructing holomorphicintegrating factors in a prior study of the attenuated transform (see [34] where Lemma 4.5 statesthat I (cid:63) ⊥ is surjective), which are so important to the present approach. Such proofs relied onpseudodifferential arguments on an extended simple compact manifold, and did not constructexplicit preimages of either operator. In order to derive and implement explicit inversions, con-structing explicit preimages becomes a necessity, and we notice here that, while writing formulas(10) and (11) in a way that emphasizes duality, we also notice that the right-hand sides involve I (cid:63) and I (cid:63) ⊥ directly. Under the assumption that the operators Id + W and Id + ( W (cid:63) ) areinvertible, this allows for an explicit construction of continuous right-inverses of I (cid:63) and I (cid:63) ⊥ . Remark 2.5.
Although it would be enough to show that Id + W is injective, which is open atpresent for general simple surfaces, it is shown in [13] that the operator W admits a bound ofthe form (cid:107) W (cid:107) L → L ≤ C (cid:107)∇ κ (cid:107) ∞ . This implies that if curvature is close enough to constant, theoperators Id + W and Id + ( W (cid:63) ) are invertible via Neumann series. Whether this qualitativeassumption covers the case of all simple surfaces remains open at present. Proposition 2.6.
Suppose the operators Id + W and Id + ( W (cid:63) ) are L ( M ) → L ( M ) in-vertible. Then for every k ∈ N , the operators R ⊥ : C k ( M ) → C k ( ∂ + SM ) and R : C k +1 ( M ) →C k ( ∂ + SM ) defined by R ⊥ := 18 π A (cid:63) + H − A − I ( Id + W ) − , R := − π A (cid:63) + H + A − I ⊥ ( Id + ( W (cid:63) ) ) − , (12)9 re continuous and satisfy I (cid:63) ⊥ R ⊥ f = f and I (cid:63) R f = f for f smooth enough.Proof. Suppose that Id + W and its adjoint are invertible. Then their inverses map any C k ( M )to itself. This is because the kernels of W, W (cid:63) are proved in [29] to be smooth so that, e.g., if f ∈ L ( M ) solves f + W f = f , where f ∈ C k ( M ) and W f is smooth, then f = f − W f is C k . Since the operators W, W (cid:63) are also L → C k continuous, then similar arguments allow toshow that ( Id + W ) − and its adjoint are C k ( M ) → C k ( M ) continuous.The relations I (cid:63) ⊥ R ⊥ = Id and I (cid:63) R = Id are straightforward to check, as a directly applica-tion of equations (10) and (11) and the invertibility of Id + W and Id + ( W (cid:63) ) .It remains to prove that (cid:107) R ⊥ f (cid:107) C k ( ∂ + SM ) ≤ C (cid:107) f (cid:107) C k ( M ) , (cid:107) R g (cid:107) C k ( ∂ + SM ) ≤ C (cid:107) g (cid:107) C k +1 ( M ) . (13)Looking at the compound expression of these operators, we see that A − and A (cid:63) + preserve C k norms since the scattering map is smooth and the function τ (x , v ) is smooth on ∂ + SM whenever ∂M is strictly convex (see [35, Lemma 4.1.1 p.115]), H ± preserve C k norms as convolutionoperators, and I : C k ( M ) → C k ( ∂ + SM ) and I ⊥ : C k +1 ( M ) → C k ( ∂ + SM ) are continuous sincethe geodesic flow is smooth. Remark 2.7.
A study of symmetries with respect to the involution α shows that p = R ⊥ f and q = R g thus constructed satisfy p ∈ V − and q ∈ V + . This is compliant with the continuitystatements (13) , as any component of p in V + would be annihilated by I (cid:63) ⊥ and any componentof q in V − would be annihilated by I (cid:63) . A crucial tool in the inversion of attenuated ray transforms is the construction of holomorphicintegrating factors , whose existence relies on the surjectivity of I (cid:63) and I (cid:63) ⊥ . In the simple Rie-mannian setting, it is proved in [26, Theorem 4.1] that the transport equation Xu = − F (for F ∈ C ∞ ( SM )) admits holomorphic solutions if and only if F is of the form F = f + X ⊥ f for some functions f , f ∈ C ∞ ( M ). Although uniqueness of such solutions may not hold (e.g.adding a constant to such a solution makes another solution), a constructive approach, inspiredin part by [26, Theorem 4.1], is to look for an ansatz, holomorphic by construction, of the form u = ( Id + iH ) p ψ + ( Id + iH ) q ψ , p is a smooth element in V − and q is a smooth element in V + , so that p ψ is odd and q ψ is even. Plugging this ansatz into Xu = − f − X ⊥ f , we obtain Xu = X ( Id + iH ) p ψ + X ( Id + iH ) q ψ = − i [ H, X ] p ψ − i [ H, X ] q ψ ( Xp ψ = Xq ψ = 0)= − i ( X ⊥ p ψ ) − i (cid:24)(cid:24)(cid:24)(cid:24)(cid:24) X ⊥ ( p ψ ) − i (cid:24)(cid:24)(cid:24)(cid:24)(cid:24) ( X ⊥ q ψ ) − iX ⊥ ( q ψ ) = i π I (cid:63) ⊥ p − i π X ⊥ I (cid:63) q (using (8)) . Therefore, a sufficient condition for u to solve Xu = − F is if p and q satisfy − i π I (cid:63) ⊥ p = f , and i π I (cid:63) q = f , which we may solve explicitly using the previous section. Using Proposition 2.6, we summarizethis construction in the following result, whose proof is straightforward and omitted. Proposition 2.8.
Under the hypotheses of Proposition 2.6, and given k ∈ N , f ∈ C k ( M ) and f ∈ C k +1 ( M ) , the function u = 2 πi [( Id + iH )( R ⊥ f ) ψ − ( Id + iH )( R f ) ψ ] , is a holomorphic solution of Xu = − f − X ⊥ f on SM , satisfying the estimate (cid:107) u (cid:107) C k ( SM ) ≤ C ( (cid:107) f (cid:107) C k ( M ) + (cid:107) f (cid:107) C k +1 ( M ) ) , where R and R ⊥ are defined in (12) . `a la Kazantzev-Bukhgeim
In [12], the authors provide reconstruction formulas for functions and vector fields from knowl-edge of their ray transforms in the case where the metric is Euclidean and the domain is theunit disk. The present section generalizes these ideas to the case of simple Riemannian surfaces.
In this first approach, we follow the idea in [12] that, if we can find a solution u (cid:63) , holomorphicwith u (cid:63) = 0 of Xu (cid:63) + au (cid:63) = − f , then projecting this equation onto Ω gives f = − η + u (cid:63) − − η − u (cid:63) − au (cid:63) = − η − u (cid:63) , u (cid:63) in terms of known data. Here and below,assuming that the operators Id + W and Id + ( W (cid:63) ) are invertible, we denote I , ⊥ the raytransform restricted to integrands of the form f + X ⊥ f , where f , f ∈ C ∞ ( M ). If D := I f + I ⊥ f = I ( f + X ⊥ f ), we can reconstruct f , f (and in turn, f + X ⊥ f ) from D bycarrying out the following steps:1. Compute the V + / V − decomposition of D . This decomposition will be given by ( D ± α (cid:63) D ) ∈ V ± , with α the antipodal scattering relation.2. Reconstruct f from the projection of D onto V + by inverting (10).3. Reconstruct f from the projection of D onto V − by inverting (11).In what follows, we summarize the above procedure by writing I − , ⊥ D = f + X ⊥ f . Theorem 3.1.
Let ( M, g ) a simple Riemannian surface with boundary and f, a ∈ C ∞ ( M, R ) .Then f is uniquely determined by its attenuated geodesic transform via the reconstruction formula f = 2 iη − ( I m( e w h (cid:48) ψ )) , h (cid:48) := 12 ( D − w (cid:48) ) | ∂ + SM , (14) where w is an odd, holomorphic solution of Xw = − a , D = (( Id − iH )( e − w ( u | ∂SM )) − with u | ∂ − SM = 0 and u | ∂ + SM = I a f , and w (cid:48) is a holomorphic solution of Xw (cid:48) = − I − , ⊥ ( A (cid:63) − D ) , givenby Proposition 2.8.Proof. Step 1: find u (cid:63) a holomorphic solution of Xu (cid:63) + au (cid:63) = − f . Call u the solutionto Xu + au = − f , u | ∂ − SM = 0 so that u | ∂ + SM = I a f . Let w a holomorphic, odd solution of Xw = − a . Then v = e − w u solves the transport problem Xv = e − w ( Xu − ( Xw ) u ) = e − w ( Xu + au ) = − e − w f, v | ∂ + SM = e − w I a f, v | ∂ − SM = 0 , (15)The right-hand side b := e − w f is holomorphic. This is because, since the product of holomorphicfunctions is holomorphic, so are (convergent) powers series of holomorphic functions. Moreover,plugging the expansion w = w + w + . . . ( w j ∈ Ω j ) into the exponential, we see that thefiberwise decomposition of b (x , v ) looks like b = f (1 + w + w . . . ) = f (cid:124)(cid:123)(cid:122)(cid:125) b + f w (cid:124)(cid:123)(cid:122)(cid:125) b + f w (cid:124) (cid:123)(cid:122) (cid:125) b + . . . , so b = f . We now want to decompose v into v = v (cid:63) + h (cid:48) ψ , where v (cid:63) is holomorphic and h (cid:48) ψ isconstant along geodesics, then we will have that v (cid:63) solves Xv (cid:63) = − b . We proceed as follows.12irst decompose v = ( v (+) + v ( − ) ), where v ( ± ) = ( Id ± iH ) v . The function v ( − ) solves thetransport equation Xv ( − ) = X ( v − iHv )= ( Id − iH ) Xv + i [ H, X ] v = − ( Id − iH ) b + i ( X ⊥ v ) + iX ⊥ v = − ( f − i ( X ⊥ v ) ) + iX ⊥ v . Integrating along geodesics, we deduce that I ( f − i ( X ⊥ v ) ) − iI ⊥ v = A (cid:63) − v ( − ) | ∂SM , where the right-hand-side A (cid:63) − ( v ( − ) | ∂SM ) = A (cid:63) − D , where D has the expression in the statementof the theorem and is known from data I a f . Therefore the right-hand side ( f − i ( X ⊥ v ) ) − iX ⊥ v can be reconstructed upon inverting I and I ⊥ , a relation which we denote( f − i ( X ⊥ v ) ) − iX ⊥ v = I − , ⊥ A (cid:63) − D. Let w (cid:48) a second holomorphic function such that Xw (cid:48) = − ( f − i ( X ⊥ v ) ) + iX ⊥ v , constructedfollowing Proposition 2.8, that is, w (cid:48) = 2 πi ( Id + iH )( p ψ + q ψ ) with p, q solving I (cid:63) ⊥ p = ( f − i ( X ⊥ v ) ) and I (cid:63) q = iv . With w (cid:48) thus constructed, we have X ( v ( − ) − w (cid:48) ) = 0, i.e. ( v ( − ) − w (cid:48) ) is constant along geodesics.Upon rewriting v as v = ( v (+) + w (cid:48) ) + ( v ( − ) − w (cid:48) ), we see that the first term is holomorphicand the second is constant along geodesics. In other words, v is of the form v = v (cid:63) + h (cid:48) ψ , where v (cid:63) = ( v (+) + w (cid:48) ) and h (cid:48) ψ = ( v ( − ) − w (cid:48) ). Additionally, with this choice of w (cid:48) , we have v (cid:63) = 12 ( v (+)0 + w (cid:48) ) = 12 ( v + 2 πi ( q ψ ) ) = 12 ( v + iI (cid:63) q ) = 0 . Finally, defining u (cid:63) = v (cid:63) e w , we see that u (cid:63) is holomorphic as the product of holomorphicfunctions, and, using the last equation, we see that u (cid:63) = ( v (cid:63) e w ) = v (cid:63) = 0, and that it solves Xu (cid:63) + au (cid:63) = − be w = − f. Projecting the equation above into Ω , we obtain − f = η + u (cid:63) − + η − u (cid:63) + au (cid:63) = ⇒ f = − η − u (cid:63) , so it remains to show how to compute u (cid:63) in terms of known data. Step 2: express u (cid:63) in terms of known data. We now write u (cid:63) = u (cid:63) − u (cid:63) − = u (cid:63) − ( u (cid:63) ) = 2 i ( I m( u (cid:63) )) = 2 i ( I m( u (cid:63) − u )) , u (cid:63) − = 0 and the last comes from the fact that u is real valued. We now write, by definition, u (cid:63) − u = ( v (cid:63) − v ) e w = − e w h (cid:48) ψ , so that I m( u (cid:63) − u ) = − I m( e w h (cid:48) ψ ), and we arrive at the expression u (cid:63) = − i ( I m( e w h (cid:48) ψ )) . Now looking to compute h (cid:48) , since we proved that h (cid:48) ψ = ( v ( − ) − w (cid:48) ), then h (cid:48) = 12 ( v ( − ) − w (cid:48) ) | ∂ + SM = 12 v ( − ) | ∂ + SM −
12 (2 πi )[( Id + iH )( p ψ + q ψ )] | ∂ + SM , where I (cid:63) ⊥ p = g, I (cid:63) q = iv , I g − iI ⊥ v = A (cid:63) − ( v ( − ) | ∂SM ) . (16)The proof is complete. Remark 3.2.
This proof generalizes the one completed in [12] in the Euclidean case when thedomain is a disk. In order to complete the argument there, it is required to relate the fiberwiseHilbert transform with the Hilbert transform of the domain for the so-called divergent beamtransform (see [12, Lemma 4.1]), which in turn uses the singular value decomposition (SVD) ofthat operator, established in earlier references (see [12] for detail and references there).While this SVD, specific to the choice of metric and domain, does not seem straightforwardto generalize systematically, the present construction of holomorphic solutions allows here tosimplify the proof by bypassing these steps altogether.
In isothermal coordinates (x , θ ), Equation (14) takes the following form: f = − η − u (cid:63) = − e − λ ∂ (cid:16)(cid:102) u (cid:63) e λ (cid:17) , (cid:102) u (cid:63) (x) = − i π (cid:90) π e − iθ I m( e w (x ,θ ) h (cid:48) ψ (x , θ )) dθ, (17)with ∂ := ( ∂ x + i∂ y ). Constructing both holomorphic solutions w, w (cid:48) using Proposition 2.8, wearrive at the following reconstruction procedure:1. Construct w = 2 πi ( I + iH )( R ⊥ a ) ψ , odd holomorphic solution of Xw = − a .2. Compute v | ∂ + SM = e w I a f and v | ∂ − SM = 0, then v ( − ) | ∂SM = ( Id − iH )( v | ∂SM ).3. Reconstruct g and v from I g − iI ⊥ v = A (cid:63) − ( v ( − ) | ∂SM ). (inversion of I and I ⊥ )4. Construct p and q from iI (cid:63) ⊥ p = g and I (cid:63) q = iv . (inversion of I (cid:63) and I (cid:63) ⊥ , see Sec. 2.4)5. Construct h (cid:48) from (16).6. Reconstruct (cid:102) u (cid:63) then f according to (17).14 .2 Reconstruction of vector fields As in [12], the method of proof above can easily be generalized to the case of reconstruction ofvector fields, which we now present. Such a problem has applications to Doppler tomographyin media with variable index of refraction, previously studied in [11, 12]. Integrands of the form f (x , θ ) = f (x)+ f (x) cos( θ − α ) are also considered in [2] in the Euclidean setting. In our context,a smooth real-valued vector field takes the form V (x , v ) = ( F + F − ) /
2, with F ± an element ofΩ ± and F = F − . In isothermal coordinates, we may write F (x , θ ) = ( f (x) − if (x)) e − λ (x) e iθ for some real-valued functions f , f , and this would correspond to integrating the vector field f (x) ∂ x + f (x) ∂ y . Theorem 3.3.
Let ( M, g ) a simple Riemannian surface with boundary and V (x , v ) = ( F + F − ) / a smooth vector field as above with F − = F . Then at every x ∈ M where a (x) (cid:54) = 0 , F − (and hence V ) can be uniquely reconstructed from data I a V via the formula F − = − iη − (cid:18) a η − ( I m( e w h (cid:48) ψ )) (cid:19) , h (cid:48) := 12 ( D − w (cid:48) ) | ∂ + SM , (18) where w is an odd, holomorphic solution of Xw = − a , D = (( Id − iH )( e − w ( u | ∂SM )) − ) | ∂SM with u | ∂ + SM = I a V and u | ∂ − SM = 0 , and w (cid:48) is a holomorphic solution of Xw (cid:48) = − I − , ⊥ ( A (cid:63) − D ) .Proof. Start from the equation Xu + au = − ( F + F − ) / , u | ∂ + SM = I a V , u | ∂ − SM = 0 . Step 1: find u (cid:63) a holomorphic solution of Xu (cid:63) + au (cid:63) = −V . Let w = w + w + . . . bean odd, holomorphic solution of Xw = − a and define v = e − w u . Then v satisfies the transportproblem Xv = − e − w ( F + F − ) / − G (x , v ) = − ( G − + G + G + G . . . ) , where G − (x , v ) = F − (x , v ) / G (x) = − w F − /
2. Split v into an holomorphic andantiholomorphic part, i.e. look at v = ( v (+) + v ( − ) ), where v ( ± ) = ( Id ± H ) v . v ( − ) satisfies thetransport equation Xv ( − ) = X ( v − iHv ) = ( Id − iH ) Xv + i [ H, X ] v = − ( Id − iH ) G + i ( X ⊥ v ) + iX ⊥ v = − ( G − i ( X ⊥ v ) ) + iX ⊥ v − G − . The data D defined in the statement is D = v ( − ) | ∂SM . Using the Hodge decomposition, we nowwrite G − = Xg + X ⊥ h for g, h two functions on M , where g fulfills the additional prescription g | ∂M = 0. Then the previous transport equation becomes X ( v ( − ) + g ) = − ( G − i ( X ⊥ v ) ) + X ⊥ ( iv − h ) . v ( − ) + g ) | ∂ + SM = v ( − ) | ∂ + SM = ( Id − iH )( e − w I a f ) | ∂ + SM is known, so that uponintegrating the transport equation along each geodesic we can see that the data gives us I ( G − i ( X ⊥ v ) ) − I ⊥ ( iv − h ) = A (cid:63) − D, with D defined as in the statement of the theorem. Now construct w (cid:48) a holomorphic solution of Xw (cid:48) = − ( G − i ( X ⊥ v ) ) + X ⊥ ( iv − h ) = − I − , ⊥ [ A (cid:63) − D ] , so that X ( v ( − ) + g − w (cid:48) ) = 0, i.e. there exists h defined on ∂ + SM such that v ( − ) + g − w (cid:48) = h ψ .We now decompose v = ( v (+) + w (cid:48) − g ) + ( v ( − ) − w (cid:48) + g ) = v (cid:63) + h (cid:48) ψ , where the firstterm is holomorphic and the second is constant along geodesics. In particular, we get that Xv (cid:63) = − G ( x, v ) where v (cid:63) is now holomorphic. Then defining u (cid:63) = e w v (cid:63) , we find that u (cid:63) isholomorphic and satisfies Xu (cid:63) + au (cid:63) = − ( F − + F ) / . Looking at the projections onto Ω − and Ω yields the relations η − u (cid:63) = − F − / , η − u (cid:63) + au (cid:63) = 0 , which implies the reconstruction formula, at each point where a does not vanish: F − = 2 η − (cid:18) a η − u (cid:63) (cid:19) . Step 2: obtain u (cid:63) from the measurements. This part is, again, similar to the proof ofTheorem 3.1. We write u (cid:63) = u (cid:63) − u (cid:63) − = u (cid:63) − ( u (cid:63) ) = 2 i ( I m( u (cid:63) )) = 2 i ( I m( u (cid:63) − u )) , where the first equality comes from the fact that u (cid:63) − = 0 and the last comes from the fact that u is real valued. Next we have the relation u (cid:63) − u = ( v (cid:63) − v ) e w = − e w h (cid:48) ψ , so it remains to compute h (cid:48) ψ , which after unrolling definitions, h (cid:48) = h (cid:48) ψ | ∂ + SM = 12 ( v ( − ) − w (cid:48) + g ) | ∂ + SM = 12 ( v ( − ) − w (cid:48) ) | ∂ + SM , where both terms are, again, computible from data: v ( − ) | ∂ + SM = D | ∂ + SM and w (cid:48) is a holomor-phic solution of Xw (cid:48) = − I − , ⊥ [ A (cid:63) − D ], a solution of which can be explicitely constructed followingProposition 2.8. This ends the proof. 16 A second approximate formula for functions, conditionally in-vertible via Neumann series
While Theorem 3.1 reconstructs functions exactly and, in some sense, in a “one-shot” fashion,it presents a couple of weaknesses: ( i ) it involves all values of w (x , θ ) throughout SM , whichin turn involves storing three dimensions of data when everything should be dealt with usingtwo-dimensional structures, and ( ii ) the effects of curvature (which need iterative corrections asin [18]) are not being corrected in the right place.We now propose an algorithm based on a more direct interplay of the Hilbert transform withtransport equations as in [19, 29], which leads to a Neumann-series based inversion, faster inimplementation, and valid when curvature is close enough to constant and the attenuation issmall enough in C norm. Such an algorithm is then implemented in Section 5.We first state a result about a certain family of operators generalizing the operator W f =( X ⊥ u f ) first defined in [29]. These operators appear as error operators of the next reconstruc-tion formula. We relegate the proof and some remarks about these operators to the appendix. Proposition 4.1.
Let h ∈ L ∞ ( SM ) such that V h ∈ C ( SM ) . Then the operator K h : L ( M ) → L ( M ) defined by K h f := ( X ⊥ u fh ) is well-defined and continuous, and there exists a constant C such that (cid:107) K h (cid:107) L → L ≤ C ( (cid:107)∇ κ (cid:107) ∞ (cid:107) h (cid:107) ∞ + (cid:107) V h (cid:107) C ) . (19)We now state the main result of the section. Theorem 4.2 (Iterative inversion) . If f ∈ L ( M ) and I a f denotes its attenuated GXRT withgiven attenuation a ∈ C ( M ) , then the function f satisfies the following equation f + Kf = 12 π I (cid:63) ⊥ η, η = 14 A (cid:63) + H ( e − w I a f ) − , (20) where ( e − w I a f ) − denotes extension of e − w I a f from ∂ + SM to ∂SM by oddness w.r.t. v (cid:55)→ − v , w is an odd, holomorphic function solving Xw = − a and K : L ( M ) → L ( M ) is a boundedoperator satisfying the following estimate (cid:107) K (cid:107) L → L ≤ C (cid:16) (cid:107)∇ κ (cid:107) ∞ + e (cid:107) a (cid:107) ∞ ( (cid:107) a (cid:107) ∞ (cid:107)∇ κ (cid:107) ∞ + (cid:107) a (cid:107) C ) (cid:17) . (21) Proof.
Proof of (20) . Let w a holomorphic, odd solution to Xw = − a . If u is the solutionto Xu + au = − f with boundary condition u | ∂ − SM = 0, then the function v = e − w u solvesthe transport problem Xv = − b := − f e − w with boundary condition v | ∂ − SM = 0, so that v isno other than u fe − w . Applying π X = π X ⊥ H (derived in (5)) to the equation Xv = − b , weobtain − f = − π b = π Xv = ( X ⊥ Hv ) = ( X ⊥ Hv − ) . Hv − more explicit. Hitting the transport equation satisfied by v withthe Hilbert transform H , we write a transport equation for HvX ( Hv ) = − Hb − X ⊥ v − ( X ⊥ v ) . As b is holomorphic, we have Hb = f H ( e − w ) = − if ( e − w −
1) and upon taking the even part ofthe last equation w.r.t. θ (cid:55)→ θ + π , the function Hv − solves the transport equation X ( Hv − ) = if (cosh w − − ( X ⊥ v ) , where we have used that w (x , θ + π ) = − w (x , θ ). This tells us that Hv − = u − if (cosh w − + u ( X ⊥ v ) + η ψ , where η ψ is constant along geodesics. We can deduce η from the fact that the first two terms inthe last r.h.s. vanish on ∂ − SM , so that η = ( Hv − | ∂ − SM ) ◦ α , and since v − is known from dataon ∂SM , so is η . More precisely, we have, for any (x , v ) ∈ ∂SMv − (x , v ) = 12 ( e − w I a f ) − (x , v ) := (cid:26) e − w (x ,v ) I a f (x , v ) if (x , v ) ∈ ∂ + SM, − e − w (x , − v ) I a f (x , − v ) if (x , v ) ∈ ∂ − SM. (22)Upon applying the operator π X ⊥ , we obtain f = − ( X ⊥ Hv − ) = i ( X ⊥ u f (cosh w − ) − W ( X ⊥ v ) − ( X ⊥ η ψ ) = i ( X ⊥ u f (cosh w − ) − W ( X ⊥ u fe − w ) − ( X ⊥ η ψ ) i.e. this equation takes the form f + Kf = − ( X ⊥ η ψ ) , η = ( Hv − ) | ∂ − SM ◦ α, where Kf := W ( X ⊥ u fe − w ) − i ( X ⊥ u f (cosh w − ) . (23)Note that upon rewriting u fe − w = u f + u f ( e − w − , the operator K can be rewritten as Kf = W ( X ⊥ u f ) + W ( X ⊥ u f ( e − w − ) − i ( X ⊥ u f (cosh w − ) = W f + W ( X ⊥ u f ( e − w − ) − i ( X ⊥ u f (cosh w − ) . Equation (23) now corresponds to (20) upon noticing that I (cid:63) ⊥ η = − π ( X ⊥ η ψ ) , and that I (cid:63) ⊥ η = I (cid:63) ⊥ (cid:18)
12 ( H ( e − w I a f ) − ) | ∂ − SM ◦ α (cid:19) = I (cid:63) ⊥ (cid:18) A (cid:63) + − A (cid:63) − H ( e − w I a f ) − ) (cid:19) = 14 I (cid:63) ⊥ A (cid:63) + H ( e − w I a f ) − , A (cid:63) − applied to an odd function on ∂SM makesa function with V + symmetry, which in turn is annihilated by I (cid:63) ⊥ . Proof of (21) . The theorem will be proved once we show that the operator K satisfies thebound (21). From the last equation, we see that K = W + W K h − iK h , where h = e − w − h = cosh w − K h j refers to Proposition 4.1. It is proved in [13] that (cid:107) W (cid:107) L → L ≤ C (cid:107)∇ κ (cid:107) ∞ for some constant C . By virtue of Proposition 4.1, we deduce that for j = 1 , (cid:107) K h j (cid:107) L → L ≤ C ( (cid:107)∇ κ (cid:107) ∞ (cid:107) h j (cid:107) ∞ + (cid:107) V h j (cid:107) C ) . Now we bound (cid:107) h j (cid:107) ∞ ≤ (cid:107) w (cid:107) ∞ e (cid:107) w (cid:107) ∞ and (cid:107) V h j (cid:107) C ≤ C (cid:107) V w (cid:107) C e (cid:107) w (cid:107) ∞ ≤ C (cid:107) w (cid:107) C e (cid:107) w (cid:107) ∞ for j = 1 ,
2, and together with the fact that w is constructed following Prop. 2.8, it satisfiesestimates of the form (cid:107) w (cid:107) C k ≤ C (cid:107) a (cid:107) C k , k = 0 , , , . . . Combining all these estimates together, we obtain estimate (21).
Remark 4.3.
In the absence of attenuation a ≡ , equation (20) is exactly the Fredholm equation (10) , in which case w (x , v ) ≡ so that Kf = W f and the right-hand side of (23) is the post-processing of unattenuated data. Corollary 4.4.
Under the hypotheses of Theorem 4.2, if (cid:107)∇ κ (cid:107) ∞ and (cid:107) a (cid:107) C are small enoughso that estimate (21) implies (cid:107) K (cid:107) L → L < , then f can be reconstructed from equation (20) viathe Neumann series f = 12 π ∞ (cid:88) n =0 ( − K ) n I (cid:63) ⊥ η, η = 14 A (cid:63) + H ( e − w I a f ) − . Remark 4.5.
As in the unattenuated case, this restriction on (cid:107)∇ κ (cid:107) ∞ and (cid:107) a (cid:107) C is of ratherqualitative nature and does not tell us whether all simple cases will work, and whether thisapproach would work for attenuations less than C . This is to be contrasted with successfulnumerical reconstructions below, which work for both discontinuous attenuations and cases ofmetrics arbitrarily close to non-simple. We now present a brief implementation of an inversion of (20) via a Neumann series approach.We use the code developed by the author in [18] for the unattenuated case, augmented withattenuation. The domain M is the unit disk { x = ( x, y ) , x + y ≤ } endowed with the metric g (x) = e λ (x) Id , where λ (x) = 0 . (cid:18) exp (cid:18) (x − x ) σ (cid:19) − exp (cid:18) (x + x ) σ (cid:19)(cid:19) , x = (0 . , . , σ = 0 . , and “high sound speed” near − x . The effecton geodesic curves can be seen Fig. 1. For such a domain, it can be computed that the boundaryis strictly convex and there are no conjugate points.The influx boundary ∂ + SM = S × ( − π , π ) is parameterized by x( β ) = (cid:0) cos β sin β (cid:1) and ingoingspeed direction θ ( β, α ) = β + π + α (in this case, β + π is the direction of the unit inner normal). M is represented into the unit square [ − , , discretized in an equispaced cartesian fashionwith N = 300 points along each dimension.Figure 1: Geodesics cast from the boundary point x( β ) = (cid:0) cos β sin β (cid:1) with, from left to right, β = 0 , π , π , π . The colorbar on the right indicates the values of the background sound speed c (x) = e − λ (x) .The function f and the attenuation a are displayed on Fig. 2. Note that both quantitiescontain jump singularities.Figure 2: From left to right: Unknown f , attenuation a , function n defined on ∂ + SM such that I (cid:63) ⊥ n = a (so that w = 2 πi ( Id + iH ) n ψ is a holomorphic solution of Xw = − a ). Axes for n are( β, α ) ∈ [0 , π ] × (cid:0) − π , π (cid:1) . Strategy for inversion.
Equation (20) is of the form f + Kf = L a I a f, (24)20here I a f represents forward data, L a D = π I (cid:63) ⊥ A (cid:63) + H ( e − w D ) − is an approximate inversionand we assume that K is a contraction. Once I a and L a are discretized (call their discretizedversions with the same name for simplicity), discretizing K separately and computing a finitesum of (cid:80) ∞ k =0 ( − K ) k L a I a f to reconstruct f may introduce additional numerical errors due to thefact that (24) may not be satisfied at the numerical level. A better approach is to set directly − K = Id − L a I a and to implement a finite sum of the series f = ∞ (cid:88) k =0 ( Id − L a I a ) k L a I a f. We now briefly explain how the operators I a and L a are implemented. Forward operator I a . The computation of the forward data is done by discretizing the in-flux boundary ∂ + SM into an equispaced family ( β i , α j ) ≤ i ≤ N, ≤ j ≤ N and for each data point,we compute I a f ( β i , α j ) by discretizing the system of ODEs over t ∈ [0 , τ ij ] where τ ij = τ (x( β i ) , θ ( β i , α j ))˙x = e − λ (x) (cid:18) cos θ sin θ (cid:19) , ˙ θ = e − λ (x) ( ∂ y λ (x) cos θ − ∂ x λ (x) sin θ ) , with initial conditions x(0) = (cid:0) cos β i sin β i (cid:1) , θ (0) = β i + π + α j , along which we compute I a f ( β i , α j ) = (cid:90) τ ij f (x β i ,α j ( t )) exp (cid:18)(cid:90) t a (x β i ,α j ( s )) ds (cid:19) dt, via a discrete sum, with x β i ,α j the solution of the ODE above. Computation of a holomorphic solution w of Xw = − a . Following Proposition 2.8, welook for w in the form w = 2 πi ( Id + iH ) n ψ , where n solves I (cid:63) ⊥ n = a . This requires implementing n = R ⊥ a , with R ⊥ = π A (cid:63) + H − A − I ( Id + W ) − . Following [18], ( Id + W ) − is computed viaa few iteration of a rapidly convergent Neumann series, I is the unattenuated X-ray transform.In the present case, n is represented on the right-hand plot of Figure 2. Once n is computed, asthe expression of L a only involves values of w on ∂SM , then we can compute w | ∂SM = (2 πi ( Id + iH ) n ψ ) | ∂SM = 2 πi ( Id + iH )( n ψ ) | ∂SM = 2 πi ( Id + iH ) A + n, where the Hilbert transform is processed via Fast Fourier Transform on the columns of A + n .As n has V − symmetry, A + n amounts to extending n to ∂ − SM by oddness w.r.t. θ (cid:55)→ θ + π (the latter is much more straightforward numerically).21 pproximate inverse L a . For D a data function defined on a discretization of ∂ + SM , wewish to compute L a D = π I (cid:63) ⊥ A (cid:63) + H ( e − w D ) − . The computation of w and H is explained in theprevious paragraph. A (cid:63) + h ( β, α ) is computed by combining values of h ( β, α ) and h at the endpointof the geodesic starting from coordinate ( β, α ). The main technical step is the computation of I (cid:63) ⊥ , which in isothermal coordinates can be simplified as follows (see [18, Sec. 3.1.2]): I (cid:63) ⊥ h (x) = e − λ (x) ∇ x · (cid:18) e λ ( · ) (cid:90) S (cid:18) − sin θ cos θ (cid:19) h ψ ( · , θ ) dθ (cid:19) , where ∇ x · is just divergence on a cartesian grid. Integrals in S can then be discretized usingfinite sums and the divergence is implemented using finite differences.Both computations of A (cid:63) + and I (cid:63) ⊥ require computing several geodesic endpoints, which is themain bottleneck of the code. An alternative option, trading memory for much shorter CPU time,is to compute and store all endpoints required at first, and reusing them in further Neumanniterations.Figure 3: Left: attenuated data I a f (data for Experiment 1) with f, a from Fig. 2. Right:unattenuated data I f for comparison. The differences may be observed in magnitude but notin support. Experiments.
We present two experiments, in which a and f refer to the functions displayedon Figure 2. 5 a refers to the function M (cid:51) x (cid:55)→ a (x). • Experiment 1. (low attenuation) Neumann series based reconstruction of f from I a f . • Experiment 2. (high attenuation) Neumann series based reconstruction of f from I a f .Experiment 1 successfully and stably reconstructs f within 3 Neumann iterations, as shown inFig. 4 (up to numerical inaccuracies, and given the fact that jumps can never be fully capturedexactly). As this example works even if a is discontinuous, this is to be contrasted with theregularity requirements on a from Theorem 4.2.22igure 4: From left to right: reconstruction f rc from I a f and pointwise error | f rc − f | after oneiteration, f rc from I a f and pointwise error | f rc − f | after three iterations.Experiment 2 displays a divergent Neumann series, due to the fact that attenuation 5 a is toohigh. This is in agreement with the smallness requirements of Corollary 4.4 on the attenuationcoefficient, as the operator norm of the error operator in (21) potentially grows like e (cid:107) a (cid:107) ∞ .Figure 5: From left to right: I a f (data for Experiment 2), pointwise error | f rc − f | after oneiteration, pointwise error | f rc − f | after two iterations. Divergence of the algorithm ensues dueto the appearing unstable modes. Acknowledgements
The author thanks Gunther Uhlmann for encouragement and support, Hart Smith, Colin Guil-larmou and Plamen Stefanov for helpful discussions, Larry Pierce for communicating references[22, 15], as well as the anonymous referees for valuable comments. Partial funding by NSF grantNo. 1265958 is acknowledged. 23
A certain family of operators - proof of Proposition 4.1
On simple surfaces, it is proved in [29] that the operator
W f = ( X ⊥ u f ) is smoothing on simplesurfaces, so that the equation reconstructing a function from its unattenuated ray transform isof Fredholm type.It has been observed that if f is now a function on SM , the corresponding operator nolonger has such properties. However, the operators K h introduced in Proposition 4.1 appearnaturally as error operators in Theorem 4.2, and they generalize W since W = K h with h ≡ K h may no longer be smoothing, but we can still obtain L ( M ) → L ( M ) continuitywith estimates on the norm in terms of h and the ambient curvature. We must recall some factsabout Jacobi fields (or variations of the exponential map), following [16]. For ξ ∈ T (x ,v ) SM , wemay decompose dϕ t ( ξ ) along the frame { X ( t ) , X ⊥ ( t ) , V ( t ) } at the basepoint ϕ t (x , v ) as dϕ t ( ξ ) = ζ (x , v, t ) X ( t ) + ζ (x , v, t ) X ⊥ ( t ) + ζ (x , v, t ) V ( t ) . The structure equations then provide us a differential system in t for the coefficients ζ i :˙ ζ = 0 , ˙ ζ + ζ = 0 , ˙ ζ − κ ( γ x ,v ( t )) ζ = 0 . In particular, we may express the variation fields dϕ t ( X ⊥ ) = aX ⊥ ( t ) − ˙ aV ( t ) and dϕ t ( V ) = − bX ⊥ ( t ) + ˙ bV ( t ) in term of two functions a , b defined on D , solving¨ a + κ ( γ x ,v ( t )) a = ¨ b + κ ( γ x ,v ( t )) b = 0 , (cid:104) a b ˙ a ˙ b (cid:105) (0) = [ ] . The assumption of simplicity implies that b does not vanish outside { t = 0 } and is thus positivefor all t >
0. Note the constancy of the Wronskian a ˙ b − b ˙ a ≡ Proof of Proposition 4.1.
For a smooth function φ (x , v ) on SM vanishing on ∂SM , we firstwrite, using the chain rule, X ⊥ (cid:90) τ (x ,v )0 φ ( ϕ t (x , v )) dt = (cid:90) τ (x ,v )0 ( a (x , v, t ) X ⊥ φ ( ϕ t (x , v )) − ˙ a (x , v, t ) V φ ( ϕ t (x , v ))) dt. We then rewrite (keeping (x , v ) implicit), for t (cid:54) = 0: aX ⊥ φ − ˙ aV φ = ab ( bX ⊥ φ − ˙ bV φ ) + (cid:32) − ˙ a + a ˙ bb (cid:33) V φ = − ab V ( φ ◦ ϕ t ) + 1 b V φ. Integrating this equality for t ∈ (0 , τ (x , v )) and integrating by parts over S x (using that φ vanishes at ∂SM ), we arrive at the conclusion that (cid:90) S x X ⊥ (cid:90) τ (x ,v )0 φ ( ϕ t (x , v )) dt dS ( v ) = (cid:90) S x (cid:90) τ (x ,v )0 (cid:18) V (cid:16) ab (cid:17) φ + 1 b V φ (cid:19) dt dS ( v ) . φ (x , v ) by f (x) h (x , v ) and using the fact that V ( f h ) = f V ( h ), we arrive at K h ( f )(x) = 12 π (cid:90) S x (cid:90) τ (x ,v )0 (cid:18) V (cid:16) ab (cid:17) h + 1 b V h (cid:19) f ( γ x ,v ( t )) dt dS ( v )= K h, f (x) + K h, f (x) , upon expanding the sum. The operator K h, is just as well-behaved as the operator W and forthe same reason: defining q (x , v, t ) := b (x ,v,t ) V (cid:16) a (x ,v,t ) b (x ,v,t ) (cid:17) , it is shown in [13] that | q (x , v, t ) | ≤ C (cid:107)∇ κ (cid:107) ∞ for every (x , v, t ) ∈ D . So we can rewrite K h, f (x) = 12 π (cid:90) S x (cid:90) τ (x ,v )0 q (x , v, t ) h ( ϕ t (x , v )) f ( γ x ,v ( t )) b (x , v, t ) dt dS ( v )= 12 π (cid:90) M q (x , v (y) , t (y)) h ( ϕ t (y) (x , v (y)) f (y) dM y , so that the kernel k h, (x , y) = q (x , v (y) , t (y)) h ( ϕ t (y) (x , v (y)) of K h, is bounded (hence in L ( M × M )), i.e. the operator K h, : L ( M ) → L ( M ) is continuous (in fact, compact) with an operatornorm controlled by C (cid:107)∇ κ (cid:107) ∞ (cid:107) h (cid:107) ∞ . On to the study of the second term K h, f (x) = 12 π (cid:90) S x (cid:90) τ (x ,v )0 q (x , v, t ) b (x , v, t ) f ( γ x ,v ( t )) b dt dS ( v ) , q (x , v, t ) := V h ( ϕ t (x , v )) . The function q satisfies (cid:82) S x q (x , v, dS ( v ) = (cid:82) S x V h (x , v ) dS ( v ) = 0, so that the integral isexpected to make sense as a principal value integral. Note that near t = 0, we have b (x , v, t ) = t + t c (x , v, t ) where c is smooth on D , and 1 + t c (x , v, t ) does not vanish on D since b does notvanish outside { t = 0 } by simplicity of the surface. More precisely, we write2 πK h, f (x) = (cid:90) S x (cid:90) τ (x ,v )0 q (x , v, t ) b (x , v, t ) f ( γ x ,v ( t )) b dt dS ( v ) = Af (x) + Bf (x) + Cf (x)where, upon writing b (x , v, t ) = t + t c (x , v, t ), we define Af (x) = (cid:90) S x (cid:90) τ (x ,v )0 q (x , v, t f ( γ x ,v ( t )) b dt dS ( v ) Bf (x) = (cid:90) S x (cid:90) τ (x ,v )0 q (x , v, t ) − q (x , v, t f ( γ x ,v ( t )) b dt dS ( v ) Cf (x) = − (cid:90) S x (cid:90) τ (x ,v )0 q (x , v, t ) c (x , v, t )(2 + t c (x , v, t ))(1 + t c (x , v, t )) f ( γ x ,v ( t )) b dt dS ( v ) . Upon changing variable ( v, t ) (cid:55)→ y( v, t ) = γ x ,v ( t ) (with Jacobian dM y = b dt dS ( v )), the C termbecomes an operator with bounded kernel, i.e. L → L bounded, with operator norm controlled25y (cid:107) q (cid:107) ∞ , i.e. (cid:107) V h (cid:107) ∞ . On to the B term, we may write | q (x , v, t ) − q (x , v, | ≤ Ct (cid:107) V h (cid:107) C uniformly on D , so that the kernel of B has an integrable singularity and the B term becomesan operator with integrable kernel, i.e. L → L bounded, with operator norm controlled by (cid:107) V h (cid:107) C . On to the A term, we assume without loss of generality to be working in isothermalcoordinates. We change variable ( v, t ) (cid:55)→ y( v, t ) = γ x ,v ( t ) to make appear Af (x) = (cid:90) M q (x , v (y) , d g (x , y)) f (y) dM y = (cid:90) M q (x , v (y) , d g (x , y)) e λ (y) f (y) d y , where d y now represents the Lebesgue measure on R . Expansions near x give that d g (x , y) = e λ (y) | x − y | + | x − y | d (x , y) , ˆ θ (y) = x − y | x − y | + | x − y | ˆ θ (x , y) , this allows to rewrite A as a Calder´on-Zygmund operator of the form Af (x) = (cid:90) M q (x , ˆ θ, | x − y | f (y) d y + (cid:90) M q (x , y) | x − y | f (y) d y , ˆ θ := x − y | x − y | , (25)where q is uniformly bounded by C (cid:107) V h (cid:107) C . By virtue of [17, Theorem XI.3.1], the first termtogether with the zero mean value condition (cid:82) S q (x , ˆ θ, dθ = 0 is an operator L → L continuous, with an operator norm bounded by C sup x ∈ M (cid:107) q (x , · ) (cid:107) L ( S ) , in turn bounded by C (cid:107) V h (cid:107) ∞ . The second term of (25) is another weakly singular operator whose operator normcan be bounded by C (cid:107) V h (cid:107) C as well. References [1]
E. V. Arbuzov, A. L. Bukhgeim, and S. G. Kazantsev , Two-dimensional tomographyproblems and the theory of A-analytic functions , Siberian Advances in Mathematics, 8(1998), pp. 1–20. 2, 3[2]
G. Bal , On the attenuated radon transform with full and partial measurements , InverseProblems, 20 (2004), pp. 399–418. 2, 3, 15[3] ,
Ray transforms in hyperbolic geometry , J. Math. Pures Appl., 84 (2005), pp. 1362–1392. 2, 3[4]
D. Finch , Inside Out, Inverse problems and applications , MSRI publications. CambridgeUniversity Press, 2004, ch. The attenuated X-ray transform: Recent developments. 2[5]
B. Frigyik, P. Stefanov, and G. Uhlmann , The X-ray transform for a generic familyof curves and weights , J. Geom. Anal., 18 (2008), pp. 81–97. 2266]
C. Guillarmou, G. Paternain, M. Salo, and G. Uhlmann , The X-ray transform forconnections in negative curvature , Communications in Mathematical Physics(to appear),(2015). 2[7]
J. Guillement, F. Jauberteau, L. Kunyansky, R. Novikov, and R. Trebossen , On single-photon emission computed tomography imaging based on an exact formula for thenonuniform attenuation correction , Inverse Problems, 18 (2002), pp. L11–L19. 2[8]
J.-P. Guillement and R. Novikov , Inversion of weighted radon transforms via fi-nite fourier series weight approximations , Inverse Problems in Science and Engineering,22 (2014), pp. 787–802. 2[9]
V. Guillemin and D. Kazhdan , Some inverse spectral results for negatively curved 2-manifolds , Topology, 19 (1980), pp. 301–312. 5[10]
N. Hoell and G. Bal , Ray transforms on a conformal class of curves , Inverse Problemsand Imaging, 8 (2014), pp. 103–125. 2[11]
S. Holman and P. Stefanov , The weighted doppler transform , Inverse Problems andImaging, 4 (2010). 2, 15[12]
S. G. Kazantsev and A. A. Bukhgeim , Inversion of the scalar and vector attenuatedx-ray transforms in a unit disc , J. Inv. Ill-Posed Problems, 15 (2007), pp. 735–765. 2, 3, 4,11, 14, 15[13]
V. Krishnan , On the inversion formulas of Pestov and Uhlmann for the geodesic raytransform , J. Inv. Ill-Posed Problems, 18 (2010), pp. 401–408. 3, 9, 19, 25[14]
L. A. Kunyansky , Generalized and attenuated radon transforms: restorative approach tothe numerical inversion , Inverse Problems, 8 (1992), p. 809. 2[15]
R. Manjappa, S. Makki, R. Kumar, and R. Kanhirodan , Effects of refractive indexmismatch in optical CT imaging of polymer gel dosimeters , Medical Physics, 42 (2015). 2,23[16]
W. J. Merry and G. P. Paternain , Lecture notes: Inverse problems in geometry anddynamics , March 2011. 24[17]
S. G. Mikhlin and S. Pr¨ossdorf , Singular Integral Operators , Springer-Verlag, 1980.26[18]
F. Monard , Numerical implementation of two-dimensional geodesic X-ray transforms andtheir inversion. , SIAM J. Imaging Sciences, 7 (2014), pp. 1335–1357. 3, 17, 19, 21, 222719] ,
On reconstruction formulas for the X-ray transform acting on symmetric differentialson surfaces , Inverse Problems, 30 (2014), p. 065001. 3, 17[20]
F. Monard, P. Stefanov, and G. Uhlmann , The geodesic X-ray transform on Rieman-nian surfaces with conjugate points , Communications in Mathematical Physics, 337 (2015),pp. 1491–1513. 3[21]
F. Natterer , Inversion of the attenuated radon transform , Inverse Problems, 17 (2001),pp. 113–9. 2[22]
N. Q. Nguyen and L. Huang , Ultrasound bent-ray tomography using both transmissionand reflection data , in Proc. SPIE 9040, Medical Imaging 2014: Ultrasonic Imaging andTomography, 90400R, 2014. 2, 23[23]
R. Novikov , On the range characterization for the two-dimensional attenuated x-ray trans-formation , Inverse Problems, 18 (2002), pp. 677–700. 2[24]
R. Novikov , Weighted radon transforms and first order differential systems on the plane ,Moscow mathematical journal, 14 (2014), pp. 807–823. 2[25]
R. G. Novikov , An inversion formula for the attenuated x-ray transformation , Ark. Math.,40 (2002), pp. 147–67. (Rapport de recherche 00/05–3 Universit´e de Nantes, Laboratoirede Math´ematiques). 2, 3[26]
G. Paternain, M. Salo, and G. Uhlmann , The attenuated ray transform for connec-tions and higgs fields , Geom. Funct. Anal. (GAFA), 22 (2012), pp. 1460–1480. 2, 10[27] ,
On the range of the attenuated ray transform for unitary connections , InternationalMath. Research Notices, (online) (2013). 2[28]
L. Pestov , Questions of well-posedness of the ray tomography problems (russian) , Sib.Nauch. Izd., Novosibirsk, (2003). 3[29]
L. Pestov and G. Uhlmann , On the characterization of the range and inversion formulasfor the geodesic X-ray transform , International Math. Research Notices, 80 (2004), pp. 4331–4347. 1, 3, 6, 7, 10, 17, 24[30] ,
Two-dimensional compact simple Riemannian manifolds are boundary distance rigid ,Annals of Mathematics, 161 (2005), pp. 1093–1110. 5, 9[31]
K. Sadiq, O. Scherzer, and A. Tamasan , On the X-ray transform of planar symmetric2-tensors , preprint, (2015). arXiv:1503.04322. 2, 3[32]
K. Sadiq and A. Tamasan , On the range characterization of the two-dimensional atten-uated Doppler transform , (2014). arXiv:1411.1923. 2, 32833] ,
On the range of the attenuated radon transform in strictly convex sets , Trans. Amer.Math. Soc., (2014). (to appear) 2, 3[34]
M. Salo and G. Uhlmann , The Attenuated Ray Transform on Simple Surfaces , J. Diff.Geom., 88 (2011), pp. 161–187. 2, 3, 4, 9[35]
V. Sharafutdinov , Integral geometry of tensor fields , VSP, Utrecht, The Netherlands,1994. 10[36]
G. Uhlmann and A. Vasy , The inverse problem for the local geodesic ray transform ,Inventiones Math.(online), (2015). 2[37]