Inverting Spherical Radon Transform by a Closed-form Formula: A Microlocal Analytic Point of View
aa r X i v : . [ m a t h . C A ] J u l INVERTING SPHERICAL RADON TRANSFORM BY ACLOSED-FORM FORMULA: A MICROLOCAL ANALYTICPOINT OF VIEW
LINH V. NGUYEN
Abstract.
Let R be the restriction of the spherical Radon transform to theset of spheres centered on a hypersurface S . We study the inversion of R bya closed-form formula. We approach the problem by studying an oscillatoryintegral, which depends on the observation surface S as a parameter. Wethen derive various microlocal analytic properties of the associated closed-forminversion formula. Introduction
Let S be a smooth hypersuface in R n and f : R n → R . We define the (restricted)spherical Radon transform R ( f ) of f by the formula:(1) R ( f )( z, r ) = Z S r ( z ) f ( y ) dσ ( y ) , ( z, r ) ∈ S × R + . Here, S r ( z ) ⊂ R n is the sphere of radius r centered at z and dσ is the surfacemeasure on S r ( z ). The transform R plays an important role in thermo/photo-acoustic tomography (TAT/PAT) (see, e.g., [FPR04, FHR07, KK08, KK10]). InTAT/PAT, f is the image of the biological tissue of interest, which needs to bereconstructed. The function R ( f ) is the available data, which is (roughly) thepressure wave recorded by the transducers located on the observation surface S .The main goal of TAT/PAT is to invert R , i.e., to find f from R ( f ). This problemalso appears in several other imaging modalities, such as ultrasound tomography(see, e.g., [Nor80b, Nor80a, NL81, NL79, AGZL10]), SONAR (see, e.g., [QRS11,LQ00]) and SAR (see, e.g., [Che01, NC04, SU12]). As a result, it has attracted asubstantial amount of work.In this article, we are interested in inverting R by closed-form formulas. Severalsuch formulas have been found when S is a sphere, cylinder, hyperplane, ellipse,and polygon with certain symmetries (e.g., [FPR04, XW05, FHR07, Kun07, Ngu09,KK08, BK78, NR10, Pal11, Nat12, Hal12b, Sal12, Kun11]). The obtained formulaslook very different and, in many cases, they only coincide in the range of R (see[Ngu09] the their relations in the case of spherical surface S ).Whether closed form inversion formulas exist for a general surface S is still anopen question. The approach by [Pal11] gives an inversion formula up a compactoperator. However, the nature of that compact operator is still not understood.The approach by [Nat12, Hal12b] gives an inversion formula up to a smoothingoperator whose kernel was explicitly obtained. The research is supported by the NSF grant DMS 1212125.
An important scenario in imaging problems is the partial (limited) data phe-nomenon. That is, the data is only collected at a subset of the observation S andthe collecting time is finite. It is quite desirable to see how the formulas work inthis situation. A natural tool is microlocal analysis. However, it seems that not allcurrently found inversion formulas can be conveniently analyzed from this point ofview.We find that the inversion formula by [Kun07] for spherical S , together with itsvariation for other geometries, comes from a simple oscillatory integral. Therefore,it is suitable to be analyzed from microlocal analytic point of view. In this articlewe consider such oscillatory integral, which depends on the observation surface S asa parameter. We show that for a general surface S , the oscillatory integral definesan operator T which can be written down in the form BPR (where B and P will bedefined later). At the same time, T is a good approximate of identity operator I .As a consequence, we obtain a good approximate inverse BP of R . We then showthat the approximation works particularly well when S has some special geometries.Our presentation has two goals:1) To understand and predict the existence of inversion formula for R undersome special geometry of the observation surface. Our approach is of micro-local analytic nature, so it does not provide the proof that a formula exactlyinverts R . However, when a formula behaves micro-locally very much likethe identity, it is reasonable to expect that it may give the exact inversion.Using this idea, we predict that BP is the exact inversion formula for allconvex quadratic of surfaces. The proof of this result is the topic of an upcoming paper.2) To understand how the formula works with the limited data problem. Weemphasize that our goal is NOT to study the general problem of what canand cannot be reconstructed from the spherical Radon transform, whichwas very deeply analyzed in [SU12]. Our goal, instead, is to see how ourparticular inversion formula works microlocally under the influence of thegeometry of the observation surface S . This might help to understand theability and limitation of this inversion formula. One of our conclusions isthe inversion formula works best with the planar observation surface, interms of constructing the singularities.The article is organized as follows. In Section 2, we consider S to be the boundaryof a convex bounded domain Ω. We show that the above-mentioned operator T (which is of the form BPR ) satisfies T = I + K , where K is pseudo-differentialoperator of order at most −
1. We then go further to obtain an asymptotic expansionof T . When Ω is an elliptical domain, we show that K is an infinitely smoothingoperator. We also show that the same result holds for a parabolic domain Ω. Thisis a good indication that T provides the inversion formula for parabolic domain. InSection 3, we consider the partial data problem. Applying the same approach as forthe case of full data, we arrive to the analog T p of the operator T . We show that T p is a pseudo-differential operator and derive a simple formula for its principal symbol.As a consequence, we deduce that T p reconstructs all the “visible” singularities of f . We then show that the full symbol of T p is equal to 1 in a conic set. Therefore, T reconstructs all the singularities of f in that conic set without any distortion.Our approach work especially well when S is a hyperplane. We consider this specialcase in Section 4. We shows that T p reconstructs all the “visible” singularities of SCILLATORY INTEGRAL AND SPHERICAL RADON TRANSFORM 3 f without any distortions. This suggests that when working with partial data, theplanar observation surface is an optimal geometry.1.1. Notations and background knowledge.
For the later convenience, we nowfix some notations. Let O be a domain; we denote by D ( O ) = C ∞ ( O ) the space ofall smooth functions compactly supported inside O . The space D ′ ( O ) is the dualof D ( O ), i.e., the space of all distributions on O . Also, E ′ ( O ) is the subspace of D ′ ( O ) that contains all distributions compactly supported in O .We denote by S m ( O ) (and Ψ m ( O )) the class of amplitudes (and pseudo-differentialoperators) on O of order at most m . The space Ψ −∞ ( O ) = T m Ψ m ( O ) contains allinfinitely smoothing operators. The reader is referred to [Tr`e80, Shu01, H¨or71] fordefinitions and basic properties of amplitudes and pseudo-differential operators. Wewill use the notation WF( f ) for the wave front set of a function/distribution f . Werefer the reader to [LQ00] for an exposition in wave front set and its connectionsto spherical Radon transform.2. Convex hypersurface S Let S be the boundary of a convex bounded domain Ω. We assume that f iscompactly supported in Ω. Then, the spherical Radon transform R is a FourierIntegral Operator (FIO) of order − n (see, e.g., [Pal10]). Therefore, R extends toa bounded operator from E ′ (Ω) to D ′ ( S × R + ). From now on, we use the notation R for this extended operator.Let us introduce the operator T , which is defined by the following oscillatoryintegral: T ( f )( x ) = 12 π n Z S Z R Z R n e i ( | y − z | −| x − z | ) λ | λ | n − h z − x, ν z i f ( y ) dy dλ dσ ( z ) , x ∈ Ω . We now decompose T . For simplicity, we first assume that f ∈ C ∞ (Ω). Then, T ( f )( x ) = 12 π n Z S h z − x, ν z i Z R Z R + e i ( r −| x − z | ) λ | λ | n − R ( f )( z, r ) dr dλ dσ ( z ) . Let P : C ∞ ( R + ) → C ∞ ( R + ) be the pseudo-differential operator defined by(2) P ( h )( r ) = Z R Z R + e i ( τ − r ) λ | λ | n − h ( τ ) dτ dλ. We obtain T ( f )( x ) = 12 π n Z S h z − x, ν z i PR ( f )( z, . ) (cid:12)(cid:12) r = | x − z | dσ ( z ) . Let B : C ∞ ( S × R + ) → C ∞ (Ω) be the back-projection type operator B ( g ) = 12 π n Z S h z − x, ν z i g ( z, | x − z | ) dσ ( z ) . We arrive to the following decomposition:(3) T ( f ) = BPR ( f ) . LINH V. NGUYEN
Since R , B are FIOs (e.g., [Pal10]) and P ∈ Ψ n − ( R + ), they all extend continuouslyto the corresponding spaces of distributions. We will see later that T is a pseudo-differential operator. Hence, it also extends continuously to E ′ (Ω) and the aboveidentity holds for all f ∈ E ′ (Ω).Working out the formula of operator P more explicitly, we obtain the followingformula of T : T ( f )( x ) = ( − n − π n − R S h z − x, ν z i (cid:2) D n − r − R ( f )( z, r ) (cid:3) (cid:12)(cid:12) r = | z − x | dσ ( z ) , for odd n, ( − n − π n R S h z − x, ν z i R R + D n − [ τ − R ( f )( z,τ ) ] | x − z | − τ τ dτ dσ ( z ) , for even n, where D = 12 r ddr . When S is a sphere or ellipse, T is the inversion formula obtained Kunyansky[Kun07] (for spherical S ), Natterer [Nat12] and Haltmeier [Hal12b] (for elliptical S ). The the factor h z − x, ν z i in the formula of T is very useful. It simplifies thesymbol calculus of T as shown in the following theorem: Theorem 2.1.
Assume that S is the boundary of a convex bounded domain Ω .Then, T is a pseudo-differential operator whose principal symbol is equal to .Proof. Let us recall the formula: T ( f )( x ) = 12 π n Z S Z R Z R n e i [ | y − z | −| x − z | ] λ | λ | n − h z − x, ν z i f ( y ) dy dλ dσ ( z ) . We observe the simple identity:(4) | y − z | − | x − z | = 2 h x − y, z − x i + | x − y | . Therefore, the Schwartz kernel of T is: K ( x, y ) = 12 π n Z S Z R + e i [ h x − y, z − x ] λ i + | x − y | λ ] | λ | n − h z − x, ν z i dλ dσ ( z )+ 12 π n Z S Z R − e i [ h x − y, z − x ] λ i + | x − y | λ ] | λ | n − h z − x, ν z i dλ dσ ( z )= X ± K ± . We first consider K + . Let us introduce the change of variables( z, λ ) ∈ S × R + −→ ξ = 2[ z − x ] λ ∈ R n . (It resembles the change from the polar to cartesian coordinates.) Straight forwardcalculations show that: dξ = 2 n | λ n − h z − x, ν z i | dλ dσ ( z ) = 2 n | λ | n − h z − x, ν z i dλ dσ ( z ) . The last equality holds since h z − x, ν z i > K + ( x, y ) = 12(2 π ) n Z R n e i h h x − y,ξ i + | x − y | | ξ | | x − z +( x,ξ ) | i dξ. SCILLATORY INTEGRAL AND SPHERICAL RADON TRANSFORM 5
Here, z + ( x, ξ ) is the intersection of S with the ray { x + tξ : t > } . The point z + ( x, ξ ) is uniquely determined since Ω is convex and x ∈ Ω.Similarly, we obtain: K − ( x, y ) = 12(2 π ) n Z R n e i h h x − y,ξ i−| x − y | | ξ | | x − z − ( x,ξ ) | i dξ, where z − ( x, ξ ) is the intersection of S with the ray { x + tξ : t < } .Due to standard theory of FIO (e.g., [Sog93, Theorem 3.2.1]), K ± are kernels ofpseudo-differential operators whose principal symbol is . Hence, K is the kernelof a pseudo-differential operator whose principal symbol is 1. This concludes theproof of the theorem. (cid:3) Corollary 2.2. T extends continuously to E ′ (Ω) → D ′ (Ω) . Moreover, T = I + K ,where K is pseudo-differential operator of order − . Since the operators B , P , R extend continuously to corresponding spaces of dis-tributions, the identity T ( f ) = BPR ( f ) holds true for all f ∈ E ′ (Ω). Corollary 2.2simply says that BP is a parametrix of R . In imaging applications, it is reasonableto consider T ( f ) = BPR ( f ) as the approximate of the image f . Although the error K ( f ) = BPR ( f ) − f might not be small, it is smoother than f . Therefore, T ( f )preserves the main part (top order) of singularities of f . For example, let us assumethat f has some jump singularities. Then, computing BPR ( f ), one recovers thesame jumps at the same locations. The Schwartz kernel of K = T − I was explicitlyfound in [Nat12] and [Hal12b]. However, our approach is more convenient frommicro-local analytic point of view. Indeed, we can compute the full symbol of T : Theorem 2.3.
Let z ± ( x, ξ ) be defined as in the proof of Theorem 2.1 and J k ( x, ξ ) = | ξ | k | x − z + ( x, ξ ) | k + ( − k | ξ | k | x − z − ( x, ξ ) | k . The the full symbol σ ( x, ξ ) of T is given by: σ ( x, ξ ) ∼ ∞ X k =1 ( − i ) k k k ! ∆ kξ J k ( x, ξ ) , ∀ ( x, ξ ) ∈ T ∗ Ω \ . Here, ∆ ξ is the Laplacian applying to the variable ξ and ∆ kξ is the k -th power of ∆ ξ . The above theorem follows from the following asymptotic behavior of T : Lemma 2.4.
Let T k be defined by the formula: (5) T k ( f )( x ) = 12(2 π ) n Z R n Z R n e i h x − y,ξ i ∆ kξ J k ( x, ξ ) f ( y ) dy dξ, x ∈ Ω . Then, for all N ∈ N : T − " I + N X k =1 ( − i ) k k k ! T k ∈ Ψ − ( N +1) (Ω) . The above lemma can be restated as:(6)
T ≡ I + ∞ X k =1 ( − i ) k k k ! T k (up to infinitely smoothing factor) . LINH V. NGUYEN
It is similar to [Bey84, Theorem 4], which was stated for the generalized Radontransform.
Proof of Lemma 2.4.
From the proof of Theorem 2.1, we obtain that T = T + + T − ,where: T + f ( x ) = 12(2 π ) n Z R n Z R n e i h h x − y,ξ i + | x − y | | ξ | | x − z +( x,ξ ) | i f ( y ) dy dξ, and T − f ( x ) = 12(2 π ) n Z R n Z R n e i h h x − y,ξ i−| x − y | | ξ | | x − z − ( x,ξ ) | i f ( y ) dy dξ. We now analyze T + . Using Taylor’s formula, we obtain: e i | x − y | | ξ | | x − z +( x,ξ ) | = N X k =0 i k k k ! | x − y | k | ξ | k | x − z + ( x, ξ ) | k + i N +1 N +1 N ! | x − y | N +1) | ξ | N +1 | x − z + ( x, ξ ) | N +1 1 Z (1 − t ) N e it | x − y | | ξ | | x − z +( x,ξ ) | dt. We arrive to T + f ( x ) = 12(2 π ) n N X k =0 i k k k ! Z R n Z R n | x − y | k e i h x − y,ξ i | ξ | k | x − z + ( x, ξ ) | k f ( y ) dy dξ + A ( f )( x ) . Here, A ( f )( x ) = i N +1 π ) n N +1 N ! Z (1 − t ) N A t ( f )( x ) dt, where A t ( f )( x ) is equal to: Z R n Z R n e i (cid:16) h x − y,ξ i + t | x − y | | ξ | | x − z +( x,ξ ) | (cid:17) | x − y | N +1) | ξ | N +1 | x − z + ( x, ξ ) | N +1 f ( y ) dy dξ. Since | x − y | k e i h x − y,ξ i = ( − k ∆ kξ e i h x − y,ξ i , we obtain T + f ( x ) = 12(2 π ) n N X k =0 ( − i ) k k k ! Z R n Z R n ∆ kξ (cid:2) e i h x − y,ξ i (cid:3) | ξ | k | x − z + ( x, ξ ) | k f ( y ) dy dξ + A ( f )( x ) . Taking integration by parts with respect to ξ , we obtain: T + f ( x ) = 12(2 π ) n N X k =0 ( − i ) k k k ! Z R n Z R n e i h x − y,ξ i ∆ kξ (cid:20) | ξ | k | x − z + ( x, ξ ) | k (cid:21) f ( y ) dy dξ + A ( f )( x ) . For any t ∈ R , due to the special form of its phase function, A t is a pseudo-differential operator (see, e.g., [Shu01]). Moreover, the amplitude function of A t is of class S N +1 (Ω × Ω) and vanishes up to order 2( N + 1) on the diagonal ∆ = SCILLATORY INTEGRAL AND SPHERICAL RADON TRANSFORM 7 { ( x, x ) : x ∈ Ω } . Due to the standard theory of FIO (e.g., [H¨or03, Proposition1.2.5]), we obtain A t ∈ Ψ − ( N +1) (Ω). Therefore, A ∈ Ψ − ( N +1) (Ω).Similarly, we obtain: T − f ( x ) = 12(2 π ) n N X k =0 i k k k ! Z R n Z R n e i h x − y,ξ i ∆ kξ (cid:20) | ξ | k | x − z − ( x, ξ ) | k (cid:21) f ( y ) dy dξ + B ( f )( x ) , where B ∈ Ψ − ( N +1) (Ω).Adding up the formulas of T + and T − , we obtain T = I + N X k =1 ( − i ) k k k ! T k + A + B. This finishes the proof of the theorem. (cid:3)
For some special geometries of S , T k is very easy to deal with. Then, some niceproperties of T can be drawn from Lemma 2.4, or equivalently the identity (8). Asexamples, we now consider the cases S is an ellipse and elliptical parabola.2.1. Elliptical domain.Theorem 2.5.
Assume that S is an ellipse. Then, K = T − I is an infinitelysmoothing operator. That is, T ( f ) − f ∈ C ∞ (Ω) for all f ∈ E ′ (Ω) .Proof of Theorem 2.5. From Lemma 2.4, it suffices to prove that T k ≡ k ≥
1. Without loss of generality, we can assume that S is defined by: S = { z ∈ R n : n X i =1 ω i z i = 1 } , where ω ′ i s are some fixed positive numbers. For two vectors x, y ∈ R n , we definethe inner product h x, y i = n X i =1 ω i x i , and the (scaled) norm: k x k = p h x, x i . We now analyze J k ( x, ξ ) = | ξ | k | x − z + ( x, ξ ) | k + ( − k | ξ | k | x − z − ( x, ξ ) | k . To this end, let us first compute the distances | x − z + ( x, ξ ) | and | x − z − ( x, ξ ) | . Werecall that z ± ( x, ξ ) ∈ S and z ± ( x, ξ ) = x + t ± ξ , for some t + > t − <
0. Tofind t ± , we solve the equation (of the ellipse S ): k ( x + tξ ) k = 1 , or k ξ k t + 2 h x, ξ i t + ( k x k −
1) = 0 . We obtain: t ± = − h x, ξ i ± √ ∆ ′ k ξ k , LINH V. NGUYEN where ∆ ′ = (1 − k x k ) k ξ k + h x, ξ i is a (homogeneous) polynomial of degree 2 in ξ . Therefore, | ξ || x − z ± ( x, ξ ) | = | ξ || t ± ξ | = 1 | t ± | = ± h x, ξ i + √ ∆ ′ (1 − k x k ) . We arrive to the formula J k ( x, ξ ) = | ξ | k | x − z + ( x, ξ ) | k + ( − k | ξ | k | x − z − ( x, ξ ) | k = (cid:2) h x, ξ i + √ ∆ ′ (cid:3) k + ( − k (cid:2) − h x, ξ i + √ ∆ ′ (cid:3) k (1 − k x k ) k . It is easy to see that J k ( x, ξ ) is a polynomial in ξ of degree at most k . Therefore,∆ kξ J k ( x, ξ ) = 0. From (5), we obtain T k = 0 for all k ≥
1. Hence, the asymptoticformula (8) then finishes the proof. (cid:3)
The above result is weaker than the identity T ( f ) = f , which was proved in[Kun07, Nat12, Hal12b]. However, our approach has some advantage when dealingwith partial data problem, which we will consider later.2.2. Parabolic domain.
An elliptical parabola is defined (up to translation), bythe equation: S = n z ∈ R n : n − X i =1 ω i z i = z n o , where ω i > i = 0 , .., n −
1. Although it is not a closed convex surface, it almostencloses the convex domain:Ω = n z ∈ R n : n − X i =1 ω i z i < z n o , in the following sense: given x ∈ Ω, for all directions ξ ∈ R n , except for the verticalones, the line through x along direction ξ intersect S at exactly two points, on theopposite sides of x .We can still define T as the a pseudo-differential operator from E ′ (Ω) to D ′ (Ω)and the general framework presented above also applies. We now prove a similarresult to Theorem 2.5: Theorem 2.6.
Assume that S is an elliptical paraboloid. Then, K = T − I is aninfinitely smoothing operator. That is, T ( f ) − f ∈ C ∞ (Ω) for all f ∈ E ′ (Ω) .Proof of Theorem 2.5. From Lemma 2.4, it suffices to prove that T k = 0 for all k ∈ N .To simplify the writing, we introduce some notations. Let x ′ = ( x , .., x n − ) , y ′ =( y , .., y n − ) ∈ R n − , we define the inner-product h x ′ , y ′ i = n − X i =1 ω i x i y i , and the corresponding norm: k x ′ k = p h x ′ , x ′ i . SCILLATORY INTEGRAL AND SPHERICAL RADON TRANSFORM 9
Then, the equation of S reads as: S = n z : k z ′ k = z n o , where z = ( z ′ , z n ).To analyze J k ( x, ξ ) = | ξ | k | x − z + ( x, ξ ) | k + ( − k | ξ | k | x − z − ( x, ξ ) | k , we first compute the distances | x − z + ( x, ξ ) | and | x − z − ( x, ξ ) | . We recall that z ± ( x, ξ ) = x + t ± ξ ∈ S , where t + > t − <
0. Therefore, t ± are the solutionsof the equation: k x ′ + tξ ′ k = x n + tξ n . or k ξ ′ k t + (2 h x ′ , ξ ′ i − ξ n ) t + ( k x ′ k − x n ) = 0 . Since ( k x ′ k − x n ) <
0, the above equation always has two solution for all ξ ′ = 0: t ± = − (2 h x ′ , ξ ′ i − ξ n ) ± √ ∆2 | ξ ′ | , where ∆ = (2 h x ′ , ξ ′ i − ξ n ) − k ξ ′ k ( x n − k x ′ k )is a (homogeneous) polynomial of degree 2 in ξ .Noting that | x − z ± ( x, ξ ) | = | t ± | | ξ | , we obtain: | ξ || x − z ± ( x, ξ ) | = 1 | t ± | = √ ∆ ± (2 h x ′ , ξ ′ i − ξ n )( x n − k x ′ k ) . We arrive to the formula J k ( x, ξ ) = | ξ | k | x − z + ( x, ξ ) | k + ( − k | ξ | k | x − z − ( x, ξ ) | k = (cid:2) √ ∆ + (2 h x ′ , ξ ′ i − ξ n ) (cid:3) k + ( − k (cid:2) √ ∆ − (2 h x ′ , ξ ′ i − ξ n ) (cid:3) k ( x n − k x ′ k ) k . It is straight forward to see that J k ( x, . ) is a polynomial of degree at most k .Therefore, ∆ kξ J k ( x, ξ ) = 0. From (5), we obtain T k = 0 for all k ≥
1. This finishesour proof. (cid:3)
Following the proof of Theorems 2.5 and 2.6, one can easily shows that thesame result holds for any surface S defined, up to translation and rotation, by thequadratic equation:(7) m X j =1 ω j x j = n X j = m +1 ω j x j + ω n +1 . for any fixed m (1 ≤ m ≤ n ) and ω i ≥ ω ′ = ( ω , .., ω m ) = (0 , ..,
0) and ω ′′ = ( ω m +1 , .., ω n ) = (0 , .., T ( f ) = f for any surface defined by (7). However, westill do not have a proof for this stronger result. It will be the subject of study foran up coming paper. Partial data problem
In practical applications, one can only measure data on a proper subset of S andin a finite time period. Therefore, the data can be modeled as g ( z, r ) = χ ( z ) η ( r ) R ( f )( z, r ) . Here, χ, η are the spatial and time cut-off functions, respectively. That is, there arebounded subsets Γ and Γ satisfying Γ ⊂ ¯Γ ⊂ Γ ⊂ S such that χ ( z ) = ( , z ∈ Γ , , z ∈ S \ Γ , and there are R > ε > η ( r ) = ( , r ≤ R, , r ≥ R + ε, The above conditions say that the (correct) data is only available on the domainΓ and time interval [0 , R ].Let us assume that S is the boundary of a convex domain Ω. We now considerthe modified operator T p ( f ) given by the formula:12 π n Z S Z R Z R n e i ( | y − z | −| x − z | ) λ h z − x, ν z i | λ | n − η ( | y − z | ) χ ( z ) f ( y ) dy dλ dσ ( z ) . Following the derivation of the identity (3) in Section 2, we obtain the decomposi-tion: T p ( f ) = BP ( g ) . We mention that, when S is an ellipse, the effect of T p is nothing but applying theinversion formula by [Kun07, Nat12, Hal12b] to partial data g . It is reasonable toconsider T p as an approximate inversion formula for limited data case. We analyzethe effect of this inversion procedure, from micro-local analytic point of view.Following the argument in the proof of Theorem 2.1, we obtain that the Schwartzkernel of T p : K p = K + p + K − p , where K + p ( x, y ) = 12(2 π ) n Z R n e i h h x − y,ξ i + | x − y | | ξ | | x − z +( x,ξ ) | i η ( | y − z + ( x, ξ ) | ) a + ( x, ξ ) dξ, and K − p ( x, y ) = 12(2 π ) n Z R n e i h h x − y,ξ i−| x − y | | ξ | | x − z − ( x,ξ ) | i η ( | y − z − ( x, ξ ) | ) a − ( x, ξ ) dξ. Here, as in the previous section, z ± ( x, ξ ) are the intersections of S with the rays { x + tξ : t > } and { x + tξ : t < } , respectively. The functions a ± ( x, ξ ) are definedby a ± ( x, ξ ) = χ ( z ± ( x, ξ )) . Using the same argument as that for Theorem 2.1, we obtain:
Theorem 3.1.
Let S be the boundary of a convex domain Ω . Then, T p : f → BP ( g ) is a pseudo-differential operator whose principal symbol is: σ ( x, ξ ) = 12 h η ( | x − z + ( x, ξ ) | ) a + ( x, ξ ) + η ( | x − z − ( x, ξ ) | ) a − ( x, ξ ) i . SCILLATORY INTEGRAL AND SPHERICAL RADON TRANSFORM 11
Let us introduce the conic sets: A + = { ( x, ξ ) ∈ T ∗ Ω \ z + ∈ Γ and | x − z + ( x, ξ ) | ≤ R }A − = { ( x, ξ ) ∈ T ∗ Ω \ z + ∈ Γ and | x − z − ( x, ξ ) | ≤ R } . Here, T ∗ Ω \ × ( R n \
0) is the cotangent bundle of Ω excluding the zero section.Let A = A + ∪ A − , then we can write A = (T ∗ Ω \ \ [ ( z,r ) ∈ Γ × (0 ,T ) N ∗ ( S r ( z )) , where N ∗ ( S r ( z )) is the conormal bundle of the sphere S r ( z ). Following [Pal00], wecall A the visible (or audible ) zone. It is described in [LQ00] that all the singularitiesof f in A are “visible” in the data g . Therefore, they should be constructed stably.We observe that a ± ( z, ξ ) ≥ x, ξ ) ∈ T ∗ Ω \
0. Moreover, for any ( x, ξ ) ∈ A ,either a + ( x, ξ ) = 1 or a − ( x, ξ ) = 1. Therefore, σ ( x, ξ ) > x, ξ ) ∈ A .Lemma 3.1 shows that T p reconstructs the singularities of f in A stably. That isany singularities of f at ( x, ξ ) ∈ A produces a corresponding singularity in T ( f )at the same location and direction ( x, ξ ) (however, the resulted singularity maybe different from the original one in terms of magnitude). The operator T p was,indeed, used in [XWAK04] to reconstruct the singularities in the visible zone A .Lemma 3.1, thus, provides a rigorous justification for their method.If ( x, ξ ) ∈ A − ∩ A + , then σ ( x, ξ ) = 1. Therefore, the main part (top order) ofsingularities of f at any ( x, ξ ) ∈ A + ∩ A − can be reconstructed exactly. Moreover,if the geometry of S is special, we can prove more: Theorem 3.2.
Let S be an ellipse or an elliptical parabola and σ ( x, ξ ) be the fullsymbol of T p . Then σ ( x, ξ ) = 1 for all ( x, ξ ) ∈ A + ∩ A − .Proof. Similar to Lemma 2.4, we obtain the asymptotics expansion:(8) T p ≡ T ,p + ∞ X k =1 ( − i ) k k k ! T k,p (up to infinitely smoothing factor) , Here, T k,p ( f )( x ) = 12(2 π ) n Z R n Z R n e i h x − y,ξ i ∆ kξ (cid:2) J k,p ( x, y, ξ ) (cid:3) f ( y ) dy dξ, x ∈ Ω , where J k,p ( x, y, ξ ) = η ( | y − z + ( x, ξ ) | ) a + ( x, ξ ) | ξ | k | x − z + ( x, ξ ) | k + ( − k η ( | y − z − ( x, ξ ) | ) a − ( x, ξ ) | ξ | k | x − z − ( x, ξ ) | k . For ( x, ξ ) ∈ A + ∩A − and y close to x , we have a ± ( x, ξ ) = 1 and η ( | y − z ± ( x, ξ ) | ) = 1;hence, J k,p ( x, y, ξ ) = J k ( x, ξ ). As proven in Theorems 2.5 and 2.6, J k ( x, ξ ) is apolynomial in ξ of degree at most k . Therefore, for all k ≥
1, the symbol of T k iszero in A + ∩ A − . Moreover, the symbol of T is equal to 1 in A + ∩ A − . Therefore,the symbol of T is equal to 1 in A + ∩ A − . (cid:3) As a consequence of the above result, we obtain T p ( f ) − f ∈ C ∞ (Ω) for any f ∈ E ′ (Ω) satisfying WF( f ) ⊂ A + ∩ A − . Therefore, T p ( f ) reconstructs all the singularities of f in A + ∩ A − without any distortions. We notice that ( x, ξ ) ∈ A + means the singularity at ( x, ξ ) is observed in the direction ξ ; and ( x, ξ ) ∈ A − means the singularity at ( x, ξ ) is observe in the opposite direction − ξ . We, hence,conclude that: if the singularity at ( x, ξ ) is observed at both directions ± ξ , it isreconstructed perfectly. This is an interesting property that only exists for somespecial geometries of S (in the above theorem, we prove for elliptical and parabolicobservation surface S ). It would be interesting to find all such geometries.4. Hyperplane
We now consider the case S is a hyperplane. Without loss of generality, weassume that S = { x ∈ R n : x n = 0 } . Let us denote Ω = { x ∈ R n : x n > } . Of course, S does not encloses (or ”almost enclose”) Ω. Therefore, the generalframework in Sections 2 and 3 does not directly apply. However, S has the followingproperty: given x in Ω, for most directions ξ , except for the horizontal ones, theline ℓ ( x, ξ ), which pass through x along direction ξ , intersects S at exactly onepoint. This, in some sense, means that S ”half-encloses” Ω. We, hence, modify theformula of T by multiplying it by 2. That is, T ( f )( x ) = 1 π n Z S Z R Z R n e i ( | y − z | −| x − z | ) λ h z − x, ν z i | λ | n − f ( y ) dy dλ dσ ( z ) . Repeating the arguments in Section 2, we obtain the following decomposition of T : T ( f )( x ) = x n π n R ∗ PR ( f )( x ) , (9)where R ∗ is the formal adjoint of R . More explcitly, T ( f ) can be written down interms of spherical Radon transform R as follows: T ( f )( x ) = ( − n − x n π n − R S D n − (cid:2) r − R ( f )( z, r ) (cid:3) (cid:12)(cid:12) r = | x − z | dσ ( z ) , if n is odd , − n − x n π n R S R R + D n − [ τ − R ( f )( z,τ ) ] | x − z | − τ τ dτ dσ ( z ) , if n is even . It is well known that T ( f ) = f (see, e.g., [BK78, NR10, XW05, Bel09]). Thatis the above formula is an exact inversion of R . However, in practical applicationssuch as SONAR (see, e.g. [QRS11]), the partial data problem is of more importance.We assume that the data is now modeled as g ( z, r ) = χ ( z ) η ( r ) R ( f )( z, r ), where η and χ are defined as in Section 3. We then consider the corresponding operator T p ( f )( x ), which is given by:1 π n Z S Z R Z R n e i ( | y − z | −| x − z | ) λ h z − x, ν z i | λ | n − χ ( z ) η ( | y − z | ) f ( y ) dy dλ dσ ( z ) . We obtain: T p ( f )( x ) = x n π n ( R ∗ P ( g ))( x ) . This is just the application of the above inversion formula to the partial data g . SCILLATORY INTEGRAL AND SPHERICAL RADON TRANSFORM 13
We now analyze the micro-local properties of T p . For any ( x, ξ ) ∈ T ∗ Ω \ ≡ R n + × ( R n \
0) such that ξ is not parallel to the plane S , ℓ ( x, ξ ) intersects with S atexactly one point z ( x, ξ ). Let us define a ( x, ξ ) = χ ( z ( x, ξ )) . The function a ( x, ξ ) can be extended smoothly to T ∗ R n + by zero. Theorem 4.1.
The operator T p : f → x n π n P ∗ R ( g ) is a pseudo-differential operatorwhose principal symbol is σ ( x, ξ ) = a ( x, ξ ) η ( | x − z ( x, ξ ) | ) . Proof.
We can rewrite the Schwartz kernel of T p as: K p ( x, y ) = 1 π n Z S Z R e i [ h x − y, z − x ] λ i + | x − y | λ ] h z − x, ν z i | λ | n − η ( | y − z | ) χ ( z ) dλ dσ ( z ) . Let us make the following change of variables( z, λ ) ∈ S × R −→ ξ ( z, λ ) = 2[ z − x ] λ ∈ R n . Straight forward calculations show that the Jacobian of the change is 2 n h z − x, ν z i | λ | n − .Noting that λ = − ξ n x n , we arrive to K p ( x, y ) = 1(2 π ) n Z R n e i h x − y,ξ i e − i | x − y | ξn xn a ( x, ξ ) η ( | y − z ( x, ξ ) | ) dξ. (10)Standard theory of FIO (see, e.g., [Sog93, Theorem 3.2.1]) concludes the proof ofthe theorem. (cid:3) As a corollary, T p extends continuously to E ′ (Ω) → D ′ (Ω). Moreover, its princi-pal symbol σ ( x, ξ ) is equal to 1 in the visible zone: A = { ( x, ξ ) ∈ T ∗ Ω \ z ( x, ξ ) ∈ Γ , and | x − z ( x, ξ ) | < R } = (T ∗ Ω \ \ [ ( z,r ) ∈ Γ × (0 ,R ) N ∗ ( S r ( z )) . Therefore, T p ( f ) = x n π n R ∗ P ( g ) reconstructs ALL the visible singularities of f with-out distorting the top order of singularities. For example, if the function f hasa jump singularity at ( x, ξ ) ∈ A then T p ( f ) = x n π n R ∗ P ( g ) also has a jump withthe same magnitude at the same location. We actually can prove even a strongerresult: Theorem 4.2.
Let f ∈ E ′ ( R n + ) such that WF ( f ) ⊂ A , then T p ( f )( x ) − f ( x ) ∈ C ∞ (Ω) . Proof.
From (10): T p ( f )( x ) = 1(2 π ) n Z R n Z R n e i h x − y,ξ i e − i | x − y | ξn xn a ( x, ξ ) η ( | y − z ( x, ξ ) | ) f ( y ) dy dξ. Using the Taylor’s expansion for e − i | x − y | ξn xn and arguing as in Lemma 2.4, wededuce the asymptotic expansion of T p :(11) T p − ∞ X k =0 i k (2 | x n | ) k k ! T k,p ∈ Ψ −∞ (Ω) . Here, T k,p ( f )( x ) = ( − k (2 π ) n Z R n Z R n | x − y | k e i h x − y,ξ i ξ kn a ( x, ξ ) η ( | y − z ( x, ξ ) | ) f ( y ) dy dξ = 1(2 π ) n Z R n Z R n ∆ kξ h e i h x − y,ξ i i η ( | y − z ( x, ξ ) | ) a ( x, ξ ) ξ kn f ( y ) dy dξ. Taking integration by parts with respect to ξ , we obtain: T k,p ( f )( x ) = 1(2 π ) n Z R n Z R n e i h x − y,ξ i ∆ kξ (cid:2) a ( x, ξ ) η ( | y − z ( x, ξ ) | ) ξ kn (cid:3) f ( y ) dy dξ. We notice that if ( x, ξ ) ∈ A and y is close to x , then a ( x, ξ ) = 1 and η ( | y − z ( x, ξ ) | ) =1. Hence, for ( x, ξ ) ∈ A the full symbol if T k is σ k ( x, ξ ) = ∆ kξ ξ kn = ( , k = 0 , , k ≥ . From the asymptotic expansion (11), we obtain the full symbol σ ( x, ξ ) of T satisfies σ ( x, ξ ) = 1 for all ( x, ξ ) ∈ A . Therefore, T ( f ) − f ∈ D ∞ (Ω) if WF( f ) ∈ A . (cid:3) The above theorem shows that x n π n R ∗ P ( g ) reconstructs ALL the “visible” sin-gularities of f up to infinite order (i.e, perfectly without any distortion). This isquite a unique property of planar observation surface S . It suggests that the pla-nar observation surface is probably the best in terms of reconstructing “visible”singularities. Acknowledgment
The author is thankful to Professor P. Stefanov for his critical comments to thepreliminary version of the article.
References [AGZL10] Gaik Ambartsoumian, Rim Gouia-Zarrad, and Matthew A. Lewis. Inversion of thecircular Radon transform on an annulus.
Inverse Problems , 26(10):105015, 11, 2010.[And88] Lars-Erik Andersson. On the determination of a function from spherical averages.
SIAM J. Math. Anal. , 19(1):214–232, 1988.[Bel09] A. Beltukov. Inversion of the Spherical Mean Transform with Sources on a Hyperplane.
ArXiv e-prints , October 2009.[Bey84] Gregory Beylkin. The inversion problem and applications of the generalized Radontransform.
Comm. Pure Appl. Math. , 37(5):579–599, 1984.[BK78] A. L. Buhge˘ım and V. B. Kardakov. Solution of an inverse problem for an elasticwave equation by the method of spherical means.
Sibirsk. Mat. ˇZ. , 19(4):749–758,953, 1978.[Che01] Margaret Cheney. Tomography problems arising in synthetic aperture radar. In
Radontransforms and tomography (South Hadley, MA, 2000) , volume 278 of
Contemp.Math. , pages 15–27. Amer. Math. Soc., Providence, RI, 2001.[Faw85] John A. Fawcett. Inversion of n -dimensional spherical averages. SIAM J. Appl. Math. ,45(2):336–341, 1985.[FHR07] David Finch, Markus Haltmeier, and Rakesh. Inversion of spherical means and thewave equation in even dimensions.
SIAM J. Appl. Math. , 68(2):392–412, 2007.[FPR04] David Finch, Sarah K. Patch, and Rakesh. Determining a function from its meanvalues over a family of spheres.
SIAM J. Math. Anal. , 35(5):1213–1240 (electronic),2004.
SCILLATORY INTEGRAL AND SPHERICAL RADON TRANSFORM 15 [Hal12a] M. Haltmeier. Inversion of circular means and the wave equation on convex planardomains.
ArXiv e-prints , June 2012.[Hal12b] M. Haltmeier. Universal inversion formulas for recovering a function from sphericalmeans.
ArXiv e-prints , June 2012.[HKN08] Yulia Hristova, Peter Kuchment, and Linh Nguyen. Reconstruction and time rever-sal in thermoacoustic tomography in acoustically homogeneous and inhomogeneousmedia.
Inverse Problems , 24(5):055006, 25, 2008.[H¨or71] Lars H¨ormander. Fourier integral operators. I.
Acta Math. , 127(1-2):79–183, 1971.[H¨or03] Lars H¨ormander.
The analysis of linear partial differential operators. I . Classics inMathematics. Springer-Verlag, Berlin, 2003. Distribution theory and Fourier analysis,Reprint of the second (1990) edition [Springer, Berlin; MR1065993 (91m:35001a)].[KK08] Peter Kuchment and Leonid Kunyansky. Mathematics of thermoacoustic tomography.
European J. Appl. Math. , 19(2):191–224, 2008.[KK10] Peter Kuchment and Leonid Kunyansky.
Mathematics of thermoacoustic and pho-toacoustic tomography , volume 2 of
Handbook of Mathematical Methods in Imaging ,chapter 19, pages 817 – 866. Springer Verlag, 2010.[Kle03] J. Klein. Inverting the spherical Radon transform for physically meaningful functions.
ArXiv Mathematics e-prints , July 2003.[Kun07] Leonid A. Kunyansky. Explicit inversion formulae for the spherical mean Radon trans-form.
Inverse Problems , 23(1):373–383, 2007.[Kun11] Leonid Kunyansky. Reconstruction of a function from its spherical (circular) meanswith the centers lying on the surface of certain polygons and polyhedra.
Inverse Prob-lems , 27(2):025012, 22, 2011.[LQ00] Alfred K. Louis and Eric Todd Quinto. Local tomographic methods in sonar. In
Sur-veys on solution methods for inverse problems , pages 147–154. Springer, Vienna, 2000.[Nat12] F. Natterer. Photo-acoustic inversion in convex domains.
Inverse Problems Imaging ,2012.[NC04] Clifford J. Nolan and Margaret Cheney. Microlocal analysis of synthetic aperture radarimaging.
J. Fourier Anal. Appl. , 10(2):133–148, 2004.[Ngu09] Linh V. Nguyen. A family of inversion formulas in thermoacoustic tomography.
InverseProbl. Imaging , 3(4):649–675, 2009.[NL79] S.J. Norton and M. Linzer. Ultrasonic reflectivity tomography: reconstruction withcircular transducer arrays.
Ultrasonic Imaging , 1(2):154–184, 1979.[NL81] S.J. Norton and M. Linzer. Ultrasonic reflectivity imaging in three dimensions: exactinverse scattering solutions for plane, cylindrical, and spherical apertures.
BiomedicalEngineering, IEEE Transactions on , (2):202–220, 1981.[Nor80a] S.J. Norton. Reconstruction of a reflectivity field from line integrals over circularpaths.
The Journal of the Acoustical Society of America , 67(3):853–863, 1980.[Nor80b] S.J. Norton. Reconstruction of a two-dimensional reflecting medium over a circulardomain: Exact solution.
The Journal of the Acoustical Society of America , 67:1266,1980.[NR10] E. K. Narayanan and Rakesh. Spherical means with centers on a hyperplane in evendimensions.
Inverse Problems , 26(3):035014, 12, 2010.[NRT95] M. M. Nessibi, L. T. Rachdi, and K. Trimeche. Ranges and inversion formulas forspherical mean operator and its dual.
J. Math. Anal. Appl. , 196(3):861–884, 1995.[Pal00] V. P. Palamodov. Reconstruction from limited data of arc means.
J. Fourier Anal.Appl. , 6(1):25–42, 2000.[Pal10] Victor Palamodov. Remarks on the general Funk transform and thermoacoustic to-mography.
Inverse Probl. Imaging , 4(4):693–702, 2010.[Pal11] V. P. Palamodov. A uniform reconstruction formula in integral geometry.
ArXiv e-prints , November 2011.[PS02] DA Popov and DV Sushko. A parametrix for the problem of optical-acoustic tomogra-phy. In
Doklady. Mathematics , volume 65, pages 19–21. MAIK Nauka/Interperiodica,2002.[QRS11] Eric Todd Quinto, Andreas Rieder, and Thomas Schuster. Local inversion of the sonartransform regularized by the approximate inverse.
Inverse Problems , 27(3):035006, 18,2011. [Sal12] Y. Salman. An inversion formula for the spherical mean transform with data on anellipsoid in two and three dimensions.
ArXiv e-prints , August 2012.[Shu01] M. A. Shubin.
Pseudodifferential operators and spectral theory . Springer-Verlag,Berlin, second edition, 2001. Translated from the 1978 Russian original by Stig I.Andersson.[Sog93] Christopher D. Sogge.
Fourier integrals in classical analysis , volume 105 of
CambridgeTracts in Mathematics . Cambridge University Press, Cambridge, 1993.[Ste09] D. Steinhauer. A Reconstruction Procedure for Thermoacoustic Tomography in theCase of Limited Boundary Data.
ArXiv e-prints , May 2009.[SU09] Plamen Stefanov and Gunther Uhlmann. Thermoacoustic tomography with variablesound speed.
Inverse Problems , 25(7):075011, 16, 2009.[SU12] P. Stefanov and G. Uhlmann. Is a curved flight path in SAR better than a straightone?
ArXiv e-prints , May 2012.[Tr`e80] Fran¸cois Tr`eves.
Introduction to pseudodifferential and Fourier integral operators.Vol. 1 . Plenum Press, New York, 1980. Pseudodifferential operators, The UniversitySeries in Mathematics.[XW05] M. Xu and L. V. Wang. Universal back-projection algorithm for photoacoustic com-puted tomography.
Physical Review E , 71, 2005.[XWAK04] Y. Xu, L.V. Wang, G. Ambartsoumian, and P. Kuchment. Reconstructions in limited-view thermoacoustic tomography.
Medical Physics , 31:724, 2004.
Department of Mathematics, University of Idaho, Moscow, Idaho 83844, USA
E-mail address ::