Leon Ehrenpreis

Print publication date: 2003

Print ISBN-13: 9780198509783

Published to Oxford Scholarship Online: September 2007

DOI: 10.1093/acprof:oso/9780198509783.001.0001

Show Summary Details
Page of

PRINTED FROM OXFORD SCHOLARSHIP ONLINE (www.oxfordscholarship.com). (c) Copyright Oxford University Press, 2017. All Rights Reserved. Under the terms of the licence agreement, an individual user may print out a PDF of a single chapter of a monograph in OSO for personal use (for details see http://www.oxfordscholarship.com/page/privacy-policy). Subscriber: null; date: 26 February 2017

(p.681) Appendix SOME PROBLEMS OF INTEGRAL GEOMETRY ARISING IN TOMOGRAPHY

Source:
The Universality of the Radon Transform
Publisher:
Oxford University Press

Peter Kuchment

Department of Mathematics, Texas A&M University, College Station, TX 77843-3368, USA

and

Eric Todd Quinto 1

Department of Mathematics, Tufts University, Medford, MA 02155, USA

This appendix is designed to coordinate the abstract nature of this book with more concrete and practical aspects of the Radon transform.

Beautiful mathematical theorems in integral geometry form the foundation of tomography. The goal of tomography is to recover the interior structure of a nontransparent body using external measurements. In this appendix, we will outline a few tomographic problems that draw on the integral geometry in this book, including support theorems, inversion formulas, range theorems, and microlocal analysis. We will try to show how these theorems are used in tomography. We do not intend to present here any concise survey of tomography, rather a collection of examples chosen according to their connection to the topics considered in the book and to the authors’ interest. The bibliography we provide is also far from being complete, and we refer the reader to books, surveys, and collections of articles [ 14 , 4852 , 67 , 72 , 73 , 80 , 124 , 126 , 138 , 139 , 150 , 161 , 170 ] for details and further references on tomography and related problems of integral geometry.

X-ray computed tomography (CT) is the most basic type of tomography; the model is the classical Radon transform that integrates functions over lines. However, the applied mathematics is deep and the properties of the transform provide insight into all parts of integral geometry. In this section we will discuss some of the basic ideas in X-ray CT and then talk about various types of limited (p.682) data tomography. Our goal is to show some of the ways integral geometry and microlocal analysis can help one understand these problems.

The goal of X-ray CT is to get a picture of the internal structure of an object by X-raying the object from many different directions. Let f : R3 → R be the density function of the object; that is, f(x) is the density of the object at the point x ∈ R3. Mathematically, the goal of X-ray CT is to recover the density f from these measurements. As X-rays travel on a line L from the X-ray source through the object to an X-ray detector, they are attenuated by the material on the line they travel along (scatter can be neglected in most X-ray CT problems). If the X-rays are monochromatic, then the attenuation at xL is proportional to the density f(x) (with proportionality constant c > 0). If I(x) is the number of X-ray photons in the beam when it arrives at x, then the intensity in a small segment of length Δx is decreased by a multiplicative factor of cf(x)Δx, so

(A.1)
$Display mathematics$
By integrating (( A.1 )), we get:
$Display mathematics$
This integral transform R is exactly the classical Radon transform of f on the line L [( 71 ), ( 124 )]. Often, X-rays are taken in parallel slices through the object, and the problem becomes a planar problem. We will primarily consider this planar case.

Complete tomographic data is X-ray data over all lines. In practice, this means data is collected on a fairly evenly distributed set of lines throughout the object (the concept of complete data can be made precise using sampling theory as in [( 124 ), ( 126 )]). In this case, the commonly used reconstruction algorithm is filtered backprojection, which we will now discuss. We start with some standard notation. Let ω = (cosθ sinθ) be a unit vector on the plane and p be a real number. We denote by L(p, ω) = {x ∈ R2|x·ω = p} the line perpendicular to ω and p directed units from the origin. We will also often use ω = (− sin θ, cos θ), the unit vector π/2 units counterclockwise from ω. We should mention that there are two standard ways of parametrizing lines, the one we use in this section (“hyperplane” or nonparametric notation) and the one used in Section ( A.3 ) (“parallel beam” notation). The theorems and ideas in each section are easier to express using the corresponding notation.

The nonparametric Radon transform (( A.2 )) of a function f in the plane can be naturally interpreted as a function of (p, ω) (see Chapter ( 2 )):

(A.2)
$Display mathematics$

(p.683) We define now the backprojection operator, the dual Radon transform, as

$Display mathematics$
This is the integral of g over all lines through x, since L(x·ω, ω) is the line through x and perpendicular to ω. Let g$C 0 ∞$ (R × S 1); then the Fourier transform of g in the p variable is
$Display mathematics$
For f$C 0 ∞$ (R 2) we define the two-dimensional Fourier transform by
$Display mathematics$
This allows us to define the Riesz potential I −1 as the operator with Fourier multiplier |σ|:
$Display mathematics$
A variation of the usual inversion formula for the nonparametric Radon transform (Section 2.1 ) is

Theorem A.1 [( 124 ), ( 151 ), ( 160 )] Let $f ∈ C 0 ∞ ( R 2 )$. Then $f = 1 4 π R * ( I − 1 R f ) ( x ) .$

The theorem is applied in practice by truncating and smoothing the multiplier |σ| in I−1 and writing this truncated multiplier as a convolution operator in p [( 124 ), ( 126 ), ( 151 )]. The resulting approximate inversion algorithm becomes fR*(Φ*p Rf) where Φ is the inverse Fourier transform of the truncated multiplier and *p denotes convolution in the p variable.

Proof The proof uses some of the key elementary formulas for the Radon transform. It begins with the projection-slice theorem (which is discussed in Chapter ( 2 ) (( 2.3 ))ff.):

(A.3)
$Display mathematics$
Equation ( A.3 ) is proven using Fubini's theorem. It is an elementary application of Radon transforms on spreads: to prove ( A.3 ), we integrate Rf over the spread of lines perpendicular to ω with respect to exp(− ipσ) dp and then set p = x · ω for xL(p,ω).

(p.684) To finish the proof, one writes the two-dimensional Fourier inversion formula in polar coordinates and uses ( A.3 ) to get

$Display mathematics$

In Section ( A.2.2 ) we will use the properties of R as an elliptic Fourier integral operator (FIO) to develop subtle results about singularity detection in tomography. Now, we will use the projection-slice theorem ( A.3 ) to show R is, in fact, an elliptic FIO. Let J :C (R 2) → C (R×S 1) be defined by Jf(p,ω) = f(pω). We take the one-dimensional inverse Fourier transform in ( A.3 ):

(A.4)
$Display mathematics$
Of course, this integral does not converge absolutely, although the operation $ℱ σ → p − 1 ∘ J ∘ ℱ x → ξ$ is defined for $f ∈ C 0 ∞ ( R 2 )$. To make the integrals converge, one does integrations by parts as with FIOs in general to show ( A.4 ) can be defined on distributions of compact support [( 76 ), ( 168 )]. Equation ( A.4 ) shows that R is an elliptic FIO with phase function φ(x, p, ω, σ) = (p − ω ·x)σ and amplitude 1/2π.

A similar exercise using the fact that R* Rf(x) = (2/|x|)* f shows that

(A.5)
$Display mathematics$
This shows R* R is a classical elliptic pseudo-differential operator with symbol 1/(π|ξ|). The interplay between microlocal analysis and integral geometry is rich, and its history, beginning with Guillemin's seminal work, will be discussed in detail after Theorem A.3.

Wave front sets provide a subtle and elegant classification of singularities of functions and distributions, and microlocal analysis allows us to understand how FIOs transform wave front sets (compare Section ( 5.3 )).

In particular, we will use this methodology to learn about how the FIO R detects singularities and to predict how well reconstruction algorithms will perform when they have limited tomographic data. The wave front set of f ∈ 𝒟′(R n ) is defined in Chapter 5 . It is useful because it classifies not only singularities of f using points x, 0 at which f is not smooth, but also directions above each point in which f is not smooth. We also need the concept of Sobolev wave front set. The distribution f is in H 9(R2) if and only if its Fourier transform, = ℱ x→ξ f, (p.685) is in L2(R2, (1 +|ξ|2) s ). This relates global smoothness of f to integrability of its Fourier transform. A local version of this at a point x 0 ∈ R n is obtained by multiplying f by a smooth cut-off function $ψ ∈ C 0 ∞ ( R n )$ (with ψ(x 0) ≠ 0) and seeing if $( ψ f ) ˜$ is in this weighted L 2 space. However, this localized Fourier transform $( ψ f ) ˜$ gives even more specific information—microlocal information— namely, the directions near which $( ψ f ) ˜$ is in L2(R2,(1 +|ξ|2) s ). The precise definition is (see [ 135 ], p. 259)

Definition A.2 A distribution f is in the Sobolev space H s locally near x 0 ∈ R n if and only if there is a cut-off function ψ ∈ D(R n ) with ψ(x0) ≠ 0 such that the Fourier transform $( ψ f ) ˜$(ξ) ∈ L2(R n ,(1 +|ξ|2) s ). Let ξ0 ∈ R n \ 0. The distribution f is in H s microlocally near (x00) if and only if there is a cut- off function ψ ∈ D(R n ) with ψ(x0) ≠ 0 and function u(ξ) homogeneous of degree 0 zero and smooth on R n \ 0 and with u(ξ0) ≠ 0 such that the product $u ( ξ ) ( ψf ) ˜ ( ξ ) ∈ L 2 ( R n , ( 1 + | ξ | 2 ) s )$. The H s wave front set of f, WF S (f), is the complement of the set of (x0, ξ0) near which f is microlocally in H s .

It follows from this definition that, if (x0, ξ0) ∉ WF(f), then for all s, f is H s near (x0, ξ0). The wave front set and microlocal Sobolev smoothness are usually defined on T*(R n ) \ 0, the cotangent space of R n with its zero section removed because such a definition can be extended invariantly to manifolds using local coordinates. To this end, let x0 ∈ R2. If y = (y1, y 2) ∈ R2, then we let ydx = y 1dx1 +y2/2dx2 be the cotangent vector corresponding to y in $T x 0 * R 2$. Now let (p,ω) ∈ R × S l . We let dp and dω be the standard basis of $T ( p , ω ) * ( R × S 1 )$. Here, θ maps to ω = (cosθ, sinθ) provides coordinates on S 1 and the basis covector dω is the dual covector to ∂/∂θ. The wave front set is extended to distributions on R × S l using these local coordinates, and it is a subset of T*(R × S1).

The fundamental theorem that gives the relation between Sobolev wave front of a function and its Radon transform is the following.

Theorem A.3 [( 146 )] Let $f ∈ ε ′ ( R 2 ) , x 0 ∈ L ( p 0 , ω 0 ) , η 0 = dp − ( x 0 · ω 0 ⊥ ) dω ,$ and $a ∈ R n \ 0.$ The correspondence between WF(f) and WF(Rf) is

(A.6)
$Display mathematics$
Given (p 0, ω0; aη0), (x 0; aω0 dx) is uniquely determined by ( A.6 ). Moreover, f ∈ H s microlocally near (x0; aω0 dx) if and only if Rf ∈ H s+1/2 microlocally near (p0, ω0; aη0). Singularities of Rf above (p 0, ω0) give no information about singularities of f above points not on L(p 0, ω0) or at points on this line in directions not conormal to the line.

This theorem is a corollary of fundamental results of Guilleman [( 59 )], [( 61 )] on the microlocal analysis of R. The general results will first be described and then the relation to Theorem A.3 will be given and finally the proof will be outlined.

Guillemin first developed the microlocal analysis of the Radon transform. For general Radon transforms defined by double fibrations on manifolds with smooth (p.686) nowhere-zero measures, he proved that R is an elliptic FIO. He also proved that R*R is an elliptic pseudo-differential operator under the Bolker assumption on the Lagrangian manifold associated to the FIO R [( 60 )]. At that time, he developed a new description of FIO using push forward and pull backs generalizing those used to define Radon transforms using double fibrations [( 59 )]. This work was included in [( 61 )] (see pp 364ff.).

Quinto [( 140 )] described the symbol of R*R for all generalized Radon transforms satisfying the Bolker assumption in terms of the measures involved into the transform. He also proved more results for the hyperplane transform, including a simple expression for the symbol and a proof of invertibility for translation invariant transforms. Beylkin [( 19 )] considered Radon transforms satisfying the Bolker assumption integrating over surfaces in R n and analyzed the asymptotics of the symbol of R†KR where R† is a backprojection related to Rη and K is a type of filter. We will discuss the generalized Radon transform further in Section ( A.2.3.3 ) and Section ( A.3. )

The microlocal correspondence between WF(f) and WF(Rf) in Theorem A.3 follows from Guillemin's seminal work [( 60 ), ( 61 )]. Because R is an FIO corresponding to a specific Lagrangian manifold, ( A.6 ) holds. Sobolev continuity is a basic property for FIOs; any FIO, such as R, of order ½ will map functions in H s of fixed compact support continuously to functions in $H loc s + 1 / 2$. If the operator is elliptic, then the function must be ½ order less smooth than its image. Theorem A.3 is a refinement of this observation. Not only does R smooth of order ½ but it maps functions that are in H s microlocally near a given covector to functions that are in H s+½ near the covector given by the correspondence ( A.6 ).

The proof outline is as follows. For this classical Radon transform, the key to the correspondence ( A.6 ) is ( A.4 ), that R is an elliptic FIO with phase function φ(x, p, ω, ω) = (px · ω)ω and the explicit Lagrangian manifold associated to this operator [( 60 ), ( 61 )] (see [( 140 )] for details). The fact that (p 0, ω0; aη0) in ( A.6 ) determines (x 0, aω0dx) is a simple exercise. The assertion about H s is given in [( 146 )], and, as discussed above, is proven using general Sobolev continuity results for FIO [ 76 , 168 ].

Theorem A.3 allows one to understand what R does to singularities in a precise and rigorous way and it provides an application of microlocal analysis to tomography. It gives an exact correspondence between singularities of f and those of Rf. Moreover, it states that the singularities of Rf that are detected from the data are of Sobolev order ½ smoother than the corresponding singularities of f. For typical singularities of f (jump singularities in H ½∈) one can realistically expect the corresponding singularities of Rf not to be masked by noise.

The theorem has the corollary

Corollary A.4 Data Rf for (p, ω) arbitrarily near (p 0, ω0) detects singularities of f perpendicular to the line L(p 0, ω0) but not in other directions. If Rf is in H s+½ locally near (p 0), then f is microlocally in H s at all points on (p.687) L(p 0, ω0) at all directions conormal to L(p 0, ω0). Conversely if a direction above (p 0, ω0) is in WFs+½(Rf) then f has H s wave front set conormal to L(p 0, ω0) at the point given by ( A.6 ).

Corollary A.4 follows from the correspondence ( A.6 ) for WF and WF S and the fact that data WFs+½Rf above (p 0, ω0) provides information only about singularities of f conormal to L(p 0, ω0). Palamodov stated a closely related idea in [( 133 )]. The “tangent casting” effects of [( 156 )] are related to Corollary A.4.

A simple illustration will give an intuitive feeling for the corollary. Let f : R 2R be equal to one on the unit disk and zero outside. Then, $R f ( p , ω ) = 2 1 − p 2$ for |p| ≤ 1 and Rf(p, ω) = 0 for |p| > 1. The only lines where Rf is not smooth are those with |p| = 1 and these lines are tangent lines to the unit circle. Since WF(f) is the set of vectors normal to the unit circle, singularities of Rf precisely capture singularities of f.

Limited tomographic data is tomographic data given on some proper open subset of R × S 1. Theorem A.3 and Corollary A.4 provide a paradigm to decide which singularities of f are stably visible from limited tomographic data, and we will examine what this predicts for the three common types of limited data: limited angle CT (Section ( A.2.3.1 )), exterior CT (Section ( A.2.3.2 )), and the interior problem or local CT (Section ( A.2.3.3 )).

The basic idea is as follows. If Rf is in H s+½ near (p 0, ω0) ∈ R × S 1, then, by Corollary A.4, f must be in H s at all points of L(p 0, ω0) in directions conormal to L(p 0, ω0). By Theorem A.3, if Rf is not in H s+½ microlocally in some direction at (p 0, ω0) then f cannot be in H s conormal to L(p0, ω0) at a specific point on this line that is given by ( A.6 ). So, singularities of Rf above (p0, ω0) determine specific singularities of f conormal to L(p 0, ω0). However, data Rf arbitrarily close to (p0, ω0) does not tell anything in a stable way about singularities of f away from L(p 0, ω0) or at points of L(p 0, ω0) in directions not conormal to this line. By ( A.6 ), these singularities are visible from other data.

A typical tomographic density f is often a piecewise smooth function that is smooth on open sets with well-behaved boundaries. So, singularities of f occur at the boundaries, and the singularities are in H ½ − ∈ for ∈ > 0. By Theorem A.3, singularities of the Radon data Rf will be in H l − ∈. A limitation of this analysis is that any discrete data Rf can be considered smooth. However, standard singularities of Rf are going to have large norm in H 1 and so should be visible. Moreover, the paradigm of the preceding paragraph is observed in all typical CT reconstructions from limited data.

Ramm and Zaslavsky [( 152 )] have analyzed how the Radon transform itself behaves on functions that are smooth except at smooth boundary surfaces. They give precise asymptotics of Rf at lines tangent to boundaries depending on the curvature of the boundaries, and they have proposed a singularity detection (p.688) method using this information on the raw data [( 152 )]. This method has been tested on simulated data [( 84 )]. A more general method has been proposed [( 146 )].

We now apply Theorem A.3 and Corollary A.4 as well as some classical theorems to three common types of limited data tomography problems in the plane. We will use the microlocal results to determine which singularities of f should be visible from limited tomographic data. One can also understand stability of limited data tomography using singular value decompositions [( 35 ), ( 107 ), ( 110 ), ( 113 ), ( 114 ), ( 116 )], and they reflect the principle predicted by the microlocal analysis: singular functions associated with large singular values (which are easy to reconstruct) oscillate in directions in which the wave front is easily detectable, and vice versa.

A.2.3.1 Limited angle tomography

Let U be a proper open subset of S 1 such that U = −U. Limited angle data is tomographic data given on R × U only, and the limited angle Radon transform is the Radon transform $R : L 0 1 ( R 2 ) → L 0 1 ( R × U ) where L 0 1$ is the set of L l functions of compact support. The goal is to recover f using data Rf(p, ω) for (p, ω) ∈ R × U. This problem comes up in electron microscopy [( 33 )] and in X-ray airport inspections of luggage traveling on a belt. In both cases, one cannot acquire X-ray CT data all around the object but only in a limited angular range. Alternatively, in electron microscopy, one has a picture that shows many molecules with the same crystalline structure but oriented in random directions. This is like regular CT with complete data, but with no information about the directions [( 126 )].

Using the projection-slice theorem ( A.3 ), one can prove this limited data transform is injective for functions of compact support: knowing Rf on R × U gives the Fourier transform f̃ on the cone in R 2, V = {tω|t ∈ R, ω ∈ U}. Since f has compact support, f̃ is real analytic and this Fourier data determines f. Tuy [ 169 ] developed a reconstruction algorithm for this problem by extending this Fourier data f̃ on the cone V to R 2 by making f̃ equal to zero off of V, then taking the inverse Fourier transform. This method destroys compact support since the Fourier transforms of compactly supported functions are real analytic. His inversion is really a pseudo-differential operator that is elliptic in directions from V so standard results about wave front sets [( 76 ), ( 168 )] show that all wave fronts of f in the directions in U will be preserved.

Other appealing and successful inversion methods involve extending the Radon data from R × U to R × S 1 using the range conditions on the Radon transform [( 103 ), ( 107 )]. The basic idea is to project the data onto the range of the Radon transform with complete data (on R × S1) and then invert this consistently completed data using filtered backprojection. Here are the details. Let B be a ball of radius M > 0. To implement this, one uses the singular value decomposition for R : L 2 (B)L 2 ([− M, M] × S 1) in appropriate measures [( 30 ), ( 124 )]. One projects the incomplete data onto this basis of singular functions on L 2 ([− M, M] × S 1). (p.689) The projection is unstable because it requires an extrapolation: the singular basis is not orthogonal on [− M, M] × U but on the larger set [− M, M] θ S 1. The condition number of the projection matrix increases exponentially as the set Ubecomes smaller [( 35 ), ( 107 )]. Alternatively, if one were to construct an orthonormal basis on L 2 [− M, M] × U then the singular functions would have large norm outside of U.

This instability is intrinsic to the problem, and it is reflected in the fact that some singularities of f are not stably visible from the limited angle data. By Corollary A.4, the only singularities of f that can be detected in a stable way are those with directions in U. To see this, choose xR 2 and ω ∈ U. Any wave front of f at (x; ωdx) is detected by limited angle data because the line L(x·ω, ω) is in this data set. For the same reason, the wave front of f at (x; ωdx) for ω U will not be stably detected by this limited angle data.

Moreover, an argument of Finch related to this singularity analysis shows that inversion of limited angle or exterior data is ill-posed in any range of Sobolev norms [( 40 )]. The argument for limited angle CT is as follows. We construct a nonsmooth function f that is smooth at all points in R 2 in all directions in U. Then, by ( A.6 ), Rf is smooth on R × U. So, Rf is in every local Sobolev space above R × U but f is not in all Sobolev spaces on R 2 since f is not smooth. Therefore, inversion of the limited angle Radon transform is not continuous in any range of Sobolev norms.

Here is one way to construct such an f. Let ω ∉ U, let L be the line through the origin normal to ω and let H be a half-plane with boundary L. We take a smooth nonzero radial function of compact support, ψ and let f be the result of setting ψ to zero on H. Then, f is smooth at all points in R 2 in all directions in U but f is not in C . In fact, WF(f) is a subset of points in L with direction ωdω, i.e. conormals to L.

This instability is also illustrated by the singular functions for the limited angle problem [( 107 )]. Those corresponding to large singular values (easy to reconstruct) oscillate generally in directions in U and those corresponding to small singular values (hard to reconstruct) oscillate generally in directions outside of U.

A.2.3.2 The exterior problem

Let M > 1 and assume supp f ⊂ {x ∈ R 2||x| ≤ M}. In the exterior problem, one has data Rf(p, ω) for all ω but only for |p| > 1. By the support theorem for the Radon transform ([( 67 )], see Chapter 2 ), one can reconstruct f(x) for |x| > 1. This problem comes up in studies of the solar corona [( 6 )] in which data is total intensities of the corona of the Sun along lines exterior to the core of the Sun that go from the solar corona to the observer on Earth. Exterior data also occurs in industrial tomography of very large objects, for which X-ray data through their centres is too highly attenuated to be usable [( 148 )].

Let |x| > 1 and ω ∈ S l . Then the only singularities of f at x that are reconstructed in a stable manner are those for ω with L(x · ω, ω) in the data set, i.e. for|x · ω| > 1. Other singularities of f are not stably detected. Because (p.690) some wave front directions are not stably detectable by exterior data, inversion for the exterior problem is highly ill-posed; a translate of Finch's example in Section ( A.2.3.1 ) can be used to show inversion of the exterior transform is discontinuous in any range of Sobolev norms. Lissianoi [( 105 )] has extended Finch's Sobolev discontinuity result to show that in the exterior problem even recovery of the function in a smaller ring than where the data is given does not help to improve stability. However, logarithmic stability has been proven by Isakov [( 78 )].

Lewitt and Bates [( 103 )], Louis [( 108 )], and Natterer [( 120 )] have developed good reconstruction algorithms that use exterior data. Lewitt and Bates’ algorithm completes the exterior data by projecting it on the range of the complete Radon transform as discussed in Section ( A.2.3.1. ) The projection step is unstable because the singular functions are not orthogonal on the annulus but on a disk. Natterer's algorithm is an effective regularization method. Quinto has developed an exterior reconstruction algorithm which employs a singular value decomposition [ 134 ] for the Radon transform on domain L 2({x ∈ R 2||x| ≥ 1}) and a priori information about the shape of the object to be reconstructed. Reconstructions for “medical” phantoms are in [( 142 )] and those for industrial phantoms are in [( 143 ), ( 144 )] and real industrial data in [( 148 )]. Exactly those singularities that are supposed to be stably reconstructed are clearly imaged in the reconstruction. In Quinto's reconstruction algorithm, singularities that are not “visible” are blurred.

A.2.3.3 The interior problem and local tomography

Let M > 1 and assume supp f ⊂ {x ∈ R 2||x| ≤ M}. Interior tomographic data is data Rf(p, ω) for all ω but only for |p| < 1. Data is missing over lines outside the unit sphere, even though supp f can meet the annulus {x ∈ R 2 | l ≤ |x| ≤ M}. The goal of interior CT is to reconstruct information about f(x) for |x| ≤ 1. This problem comes up whenever scientists want information only about some portion of the object to be reconstructed, not the whole object, or in problems, such as high-resolution tomography of very small parts of objects or electron microscopy, for which it is difficult or impossible to get complete high-resolution CT data.

Simple examples (derived using the range theorem for the Radon transform) show the interior transform, the Radon transform with this limited data, is not injective. However, according to Corollary A.4, one can detect all singularities of f in |x| < 1. To see this, choose a point x inside the unit disk and choose a direction ω ∈ S 1. Then the line through x and normal to ω is in the data set for interior tomography. Therefore, by Theorem A.3, any singularity of f at (x; ωdx) is stably detected by interior data.

Maaβ [( 114 )] has developed a singular value decomposition for this problem. See also [( 110 )]. The authors of [( 103 )] have developed a reconstruction method that projects the interior data onto the range of the Radon transform with complete data and then inverts this completed data. Because of the nonuniqueness of the problem, this projection step is not unique. However, Maaβ has shown that singular functions associated to small singular values are fairly constant inside (p.691) the region of interest, the unit ball [( 114 )]. These singular functions corresponding to small singular values are difficult to reconstruct, but they do not add much detail inside the region since they are relatively constant there. This reflects the fact that all Sobolev singularities inside the region of interest are stably reconstructed.

Even though the interior transform is not invertible, we have shown that all singularities in the region can be stably reconstructed. This explains why singularity detection methods work so well for this problem. Lambda tomography [( 38 ), ( 39 ), ( 160 ), ( 171 )] is a well-developed algorithm that finds singularities of a function using interior data over a very small region. One starts with the filtered backprojection inversion formula, f = 1/4πR* I −1 Rf, and replaces I −1 by I −2 = I −1 ο I −l = −d 2 / dp 2. The beauty of this idea is that the formula

(A.7)
$Display mathematics$
is a local reconstruction formula. This is true because one needs only data Rf(p, ω) near (x · ω, ω) to calculate (d 2 / dp 2 Rf)(x · ω, ω) and then to calculate R*(d 2 / dp 2) Rf(x). Moreover, √−Δ is an elliptic pseudo-differential operator and so WF(f) = WF(√−Δf) and WF S (f) = WF s (√−Δf). This reconstruction formula takes singularities of f and makes them more pronounced; any singularity of f in H s becomes a singularity in H s−1 because λ is an elliptic pseudo-differential operator of order 1. At the same time, there is a cupping effect at the boundaries. Because of this, the developers of lambda CT chose to add a multiple of R* Rf = λ−1f to the reconstruction to smooth it out; moreover, with a good multiple, the cupping effects are decreased [( 38 ), ( 39 )].

For the X-ray transform with sources on a curve in R 3, Louis and Maafi [( 109 )] have developed a very promising generalization of lambda CT. Greenleaf and Uhlmann [( 56 ), ( 57 )] completely analyzed the microlocal properties of R* R for admissible transforms on geodesies, including this example. The microlocal analysis of this three-dimensional lambda CT operator has been investigated by Finch [( 42 )], Katsevich [( 83 )], and Lan [( 102 )]. The operator R*λ adds singularities to f because R is not well enough behaved; for the classical line transform in the plane, R*λ R is an elliptic pseudo-differential operator and therefore preserves singularities. The theorem corresponding to Theorem A.3 for the X-ray transform on lines through a curve is given in [( 146 )], and it explains which singularities are stably reconstructed by this transform. The prediction is observed in the reconstructions in [( 109 )]. The general analysis in [( 62 )] and [( 56 )] provides the microlocal results needed to prove this theorem. The microlocal properties of this line transform were also studied in [( 22 )] as a basis for support and uniqueness theorems for the X-ray transform.

Limited angle and exterior versions of lambda CT have been developed, and they are promising on simulations [( 92 )] and tests on industrial and electron microscope data (E.T. Quinto unpublished).

(p.692) Another way to deal with interior data is to truncate the data to zero by multiplying by a cut-off function which is equal to one in the region of interest. If one now applies a standard Radon inversion formula to the truncated data, one can check that the result is a pseudo-differential operator with the symbol equal to one in the local region of interest. This means that although the exact values of the function are not recoverable from the local data, its singularities are, by simple usage of truncated data in inversion formulas. One sees, for instance, that the values of jumps of the function are recovered correctly [( 86 ), ( 90 )]. This idea was further developed by Katsevich and Ramm [( 85 )] under the name of pseudo-local tomography. The authors of [( 38 )] have provided a different method to calculate the jumps of f at boundaries using lambda CT.

Madych [( 115 )] used wavelet analysis to show the strong relationship between lambda CT and pseudo-local CT and regular filtered backprojection. Authors of the papers [( 36 ), ( 153 )] have also used wavelet techniques for local tomography. They use Radon data to calculate wavelet coefficients of the density to be reconstructed. Although they need some data slightly outside the region of interest, their methods are fairly local.

Candès and Donoho [( 27 )] have defined ridgelets, wavelets that are not radially symmetric and are more sensitive to singularities in specified directions. They have developed a tomographic reconstruction algorithm using ridgelets that provides high-quality reconstructions in practice. Their ideas reflect the singularity detection predictions of Theorem A.3.

Technicians currently use the sinogram, the graph of Rf(p, ω) in rectangular coordinates, to find boundaries, but this method is subjective, and some industrial data is too homogeneous for this to work [( 149 )].

Local tomography has been developed for generalized Radon transforms [( 92 )], such as the attenuated Radon transform discussed in Section ( A.3 ) and in particular in Section ( A.3.1. ) If ν = ν(x,ω) is a nowhere-zero smooth weight, we define the generalized Radon transform of $f ∈ C 0 ∞ ( R 2 )$ to be

(A.8)
$Display mathematics$
the integral of f on the line L(p, ω) in weight ν (note that the weight ν is determined by x and ω because x · ω) = p when xL(p, ω)).

An important generalized Radon transform is the attenuated transform T μ defined in Section ( A.3 ), ( A.11 ). The weight for this transform is

$Display mathematics$
where μ(y) is the attenuation coefficient of the body at the point y.

The generalized dual transform of R ν is

$Display mathematics$
(p.693) for gC (R × S1). Let μ(x,ω) and ν(x,ω) be nowhere-zero smooth functions on R 2 × S 1.

We use this to develop a generalized lambda tomography. The calculus of FIOs shows that the generalized lambda operator (compare with ( A.7 ))

(A.9)
$Display mathematics$
is an elliptic pseudo-differential operator of order 1 with principal (top-order) symbol [( 92 )]
$Display mathematics$
The analogous theorem was proven about $R ν * R ν$ in [( 140 )]. The operator defined by ( A.9 ) is a local operator for the same reasons the lambda operator ( A.7 ) is. Because $R ν * d 2 / dp 2 R μ$ is an elliptic pseudo-differential operator, it preserves wavefront set and
$Display mathematics$
Good reconstructions of simulations using this operator from limited angle, exterior, and interior data are given in [( 92 )]. In addition, Katsevich [( 81 )] uses properties of pseudo-differential operators to provide an algorithm to detect the values of density jumps at boundaries of regions using this operator ( A.9 ).

In general, the attenuated Radon transform will have a measure that is not smooth on boundaries of regions at which the attenuation changes. This gives rise to pseudo-differential operators $R ν * d 2 / dp 2 Rν$ that have nonsmooth symbols. Katsevich studied these operators, and he proved that one can detect jumps of f that are not masked by singularities of the measure using this information, if one knows the measure [( 82 )]. This work is important even if one does not know the measure ν exactly. Because the operator is an elliptic pseudo-differential operator (albeit with a nonsmooth symbol), the singularities of f that are not masked by singularities of ν will still be visible from the data, even if the measure is unknown. One practical issue would be whether the large amount of noise in SPECT data might create problems for any singularity detection algorithm that “overdifferentiates” the data.

Some problems of medical tomography, radiation therapy, and industrial nondestructive testing naturally lead to consideration of a special type of weighted X-ray transforms. We will briefly describe one such model, referring the reader to publications [( 24 ), ( 31 ), ( 32 ), ( 95 ), ( 124 ), ( 138 ), ( 139 ), ( 158 )] for further details, examples, and references.

The so-called single photon emission computed tomography (SPECT) deals with the situation when a nontransparent body contains a distributed radiation source. The goal is to reconstruct the interior intensity distribution f(x) of the (p.694) radiation sources by conducting exterior measurements of the intensity of outgoing radiation. Let us denote by μ(x) the linear attenuation coefficient (or absorption) of the body at the point x. This means that due to the absorption a beam passing through the point x suffers at a small distance Δx the relative intensity loss equal to μ(x)Δx. Then assuming that effects of scatter are small and hence can be neglected (which is in fact not always true), the stationary one-velocity transport equation can be written as

(A.10)
$Display mathematics$
Here u(x, φ) is the density of particles at x moving in the direction of the vector νφ = (cosφ,sinφ). In other words,
$Display mathematics$
Exterior detectors are placed on the boundary. They are collimated in order to detect only the particles following a specific oriented line L. In fact, perfect collimation is impossible, so the detector is affected by all rays from a solid angle around the line L, but we will disregard this circumstance. Assuming that there are no particles entering the region from outside, one can easily solve the transport equation for the intensities measured at the detectors. This gives as the total detected intensity of the beam L the expression
(A.11)
$Display mathematics$
Here L x is the segment of the line L between the point x and the detector and dy denotes the standard linear measure on L.

The operator T μ is said to be the attenuated X-ray transform of the function f(x) and μ(x) is called the attenuation. The formula above is a nonparametric version of this transform. One can also introduce a parametric one. As mentioned in the last section, we will use a slightly different parametrization of lines in this section than we did in Section ( A.2 ). The theorems and ideas in each section are easier to express using the respective notation. Namely, let ω be the unit vector parallel to the direction of the oriented line L, ω be its orthogonal hyperplane, and y = L∩ω. Then the oriented lines L are in one-to-one correspondence with pairs (ω, y), where y ∈ ω. (Now, this notation differs from that of the previous chapters; in Chapter 6 we studied attenuations which do not depend on y.) In the present notation we write

$Display mathematics$

One notices that in contrast with the standard X-ray transform f → ∫L f(x)dx, ( A.2 ), the resulting function depends on the orientation of the line L. (p.695) In practice one often averages over the two orientations, thus arriving at a function of nonoriented lines.

In many algorithms one assumes that the attenuation is constant inside the body and zero outside. Let us also assume that the body has a known convex shape. In this case, according to [( 117 )], the attenuated transform can be reduced to a simpler one. Indeed, the integral

$Display mathematics$
is the so-called divergent beam transform of μ (the integral of μ over the ray in direction ω starting at y +tω), and it can be evaluated explicitly when μ is constant:
$Display mathematics$
Here s(ω, y) is the value r ≥ t for which y +τω belongs to the boundary of the body. Since the shape of the body is known, s(ω, y) is a known function. Then the transform T μ can be rewritten as follows:
$Display mathematics$
Here the function exp(−μs(ω, y)) is known and
$Display mathematics$
is the so-called exponential X-ray transform of the function f. (In order to avoid confusion, the reader should notice that we use notations for these transforms different from the ones in [( 124 )]:) It is necessary to mention that one needs to require that the function f(x) decays at infinity exponentially and sufficiently fast to offset the effect of the exponential weight in the integral.

Among the natural questions to ask about the attenuated and exponential X-ray transforms are:

Injectivity. Can a function of a natural class be reconstructed from its attenuated or exponential transforms? In other words, are these operators injective?

Inversion formulas in cases when injectivity is established.

Stability of inversion.

Range. Are these operators surjective in natural functional classes? Judging by the precedent of the standard Radon transform, the reader certainly should feel that the answer is probably negative, and some nontrivial range conditions should be satisfied. Then the question arises of describing the ranges of these operators. (p.696) Simultaneous reconstruction of the sources density f and attenuation μ. In most cases not only the value of the intensity f(x) distribution, but the attenuation coefficient μ(x) as well is unknown. Is it possible to extract any information about both functions from the values of T μ f or R f?

These problems will be briefly addressed below. We would also like to mention papers [( 37 ), ( 95 ), ( 131 ), ( 132 ), ( 157 )], where some interesting relations between the properties of the exponential transform and complex analysis are discovered that are beyond the scope of this text.

Uniqueness of reconstruction is obviously the first question one should ask. Non-uniqueness would mean that one is unable to recover the interior structure of the object of interest from the exterior information. When we discuss uniqueness we will assume that the function f to be reconstructed has compact support.

The problem of uniqueness for the attenuated X-ray transform in dimension 2 has proven to be hard. The first results on uniqueness were the local ones. Since for smooth μ the operator T μ (after averaging over opposite directions) is an elliptic FIO [( 19 ), ( 60 ), ( 61 ), ( 140 )], one concludes that for each point x there is a neighborhood U of x such that uniqueness holds for functions f with support in U. The idea is that locally the attenuated transform looks almost like the standard X-ray transform. It was shown in [( 118 )] that C2 smoothness of the attenuation μ is sufficient. The paper [( 66 )] contains the proof of the following statement. Let us require that the supports of all functions f under consideration belong to a fixed bounded domain Ω Assume, for instance, that Ω is the unit ball centreed at the origin. Then the transformation T μ is a semi-Fredholm operator between the Sobolev spaces $H 0 s ( Ω )$ and H s+½(S 1 × [−1, 1]) for any s (see also [( 124 )]). This means that the range of this operator is closed, and its kernel Ker T μ is finite dimensional. In the class of such operators property Ker A = 0 is stable under perturbations of small operator norm. This observation immediately leads to a local uniqueness result: for small Ω, operator Tμ is norm close to the X-ray transform, and hence has zero kernel. The next step was made by Finch [( 41 )], who managed to obtain a “semi-local” uniqueness results. Namely, he proved uniqueness under the condition ||μ|| diamΩ < 5.37. This result is sufficient for many practical situations. For instance, in medical applications it restricts the diameter of an object to 35.8 cm. The proof is nontrivial and involves energy estimates. A breakthrough came recently in brilliant works [( 8 ), ( 128 ), ( 129 )], where a positive solution of the uniqueness problem under some mild smoothness condition on the attenuation was obtained [( 8 )] and an explicit inversion formula found [( 128 ), ( 129 )] (see also [( 26 ), ( 43 ), ( 58 ), ( 96 ), ( 98 ), ( 125 )] for different derivations and implementations). We will discuss the inversion formulas in a little more detail in the next section.

It was noticed in [( 41 )] that a similar uniqueness problem in dimensions 3 and higher is trivial. The reason is that the three-dimensional attenuated X-ray (p.697) transform is the collection of the corresponding two-dimensional transforms in all affine planes. If there is a compactly supported function f(x) annihilated by the transform, then we can choose a plane which just barely touches the support of f. In this plane one can use the local uniqueness theorem to conclude that the function is in fact zero there. Proceeding further in this way, one can “eat the support away” and finally conclude that it is empty.

The problem of uniqueness was also considered for transforms with more general positive weights ν the generalized Radon transform (A.8). Here there are local uniqueness results for fairly arbitrary weights and global uniqueness theorems for some special types of weights (translation invariant, rotation invariant) [( 100 ), ( 140 ), ( 141 )]. In the case of an analytic weight uniqueness in the class of compactly supported functions follows from analytic ellipticity of the corresponding FIO [( 21 )]. Analytic ellipticity has been used to prove support theorems and uniqueness in many settings such as for Radon transforms on geodesic spheres in realanalytic manifolds [( 145 )] and for circular transforms related to the wave equation (e.g. [( 2 )]) (see also [( 1 )] and [( 172 )] for related results). However, as the famous counterexample by Boman [( 20 )] shows, the condition of infinite smoothness of the weight function ω alone does not guarantee uniqueness.

Let us move now to the much simpler case of the exponential transform Rμ. One of the most important (albeit simple) properties of the X-ray transform is the Fourier-slice (also called projection-slice) theorem (see Chapter 2 and (A.3)). There is an analog of this property for the exponential transform [( 121 ), ( 124 )], Namely, let f be a compactly supported function on the plane. A straightforward calculation shows that

(A.12)
$Display mathematics$
where the hat on the left side denotes the one-dimensional Fourier transform with respect to the linear variable s, while the tilde on the right side denotes the two-dimensional Fourier transform. This means that the one-dimensional Fourier transform of the projection data R μ f provides the values of the Fourier transform of the function f on a surface in C 2. One can see this as an indication of the possibility to use methods of the theory of functions of several complex variables. Results of [( 3 ), ( 44 ), ( 37 ), ( 89 ), ( 93 ), ( 95 ), ( 117 ), ( 121 ), ( 131 ), ( 132 )] show that the relation between two theories is indeed very deep. We will try to show this in our exposition.

Let us turn now to the problem of uniqueness of reconstruction, i.e. injectivity of the operator R μ. We will follow here considerations of [( 117 )] and [( 74 )]. Let us assume that function f is either compactly supported or decays as exp(−(μ +∈) |x|) for a positive ∈. Consider the surface

$Display mathematics$
It is straightforward to check that, except for the points where σ = 0, this surface is a totally real, smooth two-dimensional submanifold of C 2. This means, in (p.698) particular, that it is a uniqueness set for analytic functions. The formula (A.12) shows that knowing R μ f one can recover the values of the Fourier transform f̃ of the function f on the surface S μ. Due to our decay assumption, f̃ is analytic in a tubular neighborhood of R 2 containing Sμ. This implies that Rμf determines f̃ and hence f uniquely.

As soon as injectivity is established, it is natural to look for inversion formulas. By now many inversion formulas have been established for the attenuated Radon transform (see [( 8 ), ( 128 ), ( 129 )] and also [( 26 ), ( 43 ), ( 58 ), ( 91 ), ( 96 ), ( 98 ), ( 125 )] for further discussion and numerical implementation). We will address these below in this section. We would like to start with the much simpler, but still very rich, example of the exponential transform.

A.3.2.1 Exponential transform

The first explicit inversion formula for the exponential X-ray transform in the plane was provided by Tretiak and Metz in [( 167 )] (see also its discussion in [( 124 )]). An inversion procedure was also provided in [( 12 )]. Let us introduce the dual exponential X-ray transform (or exponential backprojection) $R μ #$ as follows: applied to a function g(s, ω) it produces a function on the plane according to

$Display mathematics$
One can verify the formula
$Display mathematics$
Hence, in order to reconstruct the function f from $R − μ # R μ f$ one needs to perform a deconvolution. This can be done by using appropriate generalized Riesz-type potentials on the line. Namely, let
$Display mathematics$
and $I μ − 1$ (a generalized Riesz potential) be the Fourier multiplier by ζμ(σ).

Theorem A.5 [( 167 )] Let $f ∈ C 0 ∞ ( R 2 )$. Then

(A.13)
$Display mathematics$

The original proof of this filtered backprojection-type formula in [( 167 )] and the proof provided in [( 124 )] are interesting and instructive, although rather (p.699) technical. We will now present a different simple approach based on complex analysis [( 89 ), ( 95 )].

Let us go back to the projection-slice formula (A.12). As we have already discussed, it implies that given R μ f one can obtain the values of the Fourier transform of the unknown function f on the surface S μ in C 2. If instead of S μ we dealt with R 2, we would just use the Fourier inversion formula to find f. So, the question is whether one can develop a Fourier inversion formula that uses the data on S μ instead of R 2. The simple idea is to try to use the Cauchy theorem in order to deform the surface of integration from R 2 to S μ. Let us look at the standard Fourier inversion formula in R 2:

$Display mathematics$
(an additional multiplicative constant may arise depending on normalization of the Fourier transform). Assume that f is compactly supported or decays sufficiently rapidly exponentially, so f̃(ξ) is analytic in all areas of interest. We consider the holomorphic differential form Φx = f̃(z)exp(ix · z) dz 1dz 2 in C 2. Then the Fourier inversion formula reads
$Display mathematics$
Now it is clear that using the Cauchy theorem is a good idea. There are, however, two complications. The first one is minor and can be easily dealt with. The surface S μ when projected to R 2 by zRe(z) covers R 2\{0} twice (since a point σω can be also represented as (−σ)(−ω)). To treat this problem we can restrict the values of σ in the definition of S μ to nonnegative numbers only. The second complication is much more serious and interesting. The surface S μ is not homological to R 2, in particular since it is not simply connected. That is, the points where σ = 0 constitute the boundary of a circular hole in the surface. This hole is the disk of radius μ in the imaginary plane. The idea to paste this disk to the surface does not work, since the values of f̃ on this disk are not known. The hope is to find a different disk that, on being pasted to the surface, would fix this problem. However, no matter what kind of disk we use, the values of f̃ on this disk will still be unknown. What difference then does the choice of disk make? Assume that we found an analytic disk; then the holomorphic (2,0) form Φx would vanish on this disk. Thus, for an analytic disk the integral of the form Φx over this disk vanishes independently on the values of f̃. So, we are led to a classical problem of complex analysis about pasting analytic disks to a totally real surface in C2. A little experimentation shows that such a disk Γ1 does exist and can be described as follows:
(A.14)
$Display mathematics$

(p.700) It is contained in the complex line z 2 = iz1 and hence is analytic. We now define the following part of the surface S μ:

$Display mathematics$
Defining Γ = Γ1 ∪ Γ2, applying the Cauchy theorem, and using the equality ∫Γ 1Φx = 0, we get an inversion formula
$Display mathematics$
Writing the integral ∫Γ 2 Φx explicitly using the definition of the form and coordinates σ and ω on Sμ, one arrives at the inversion formula (A. 13). The details are easily workable and can be found in [( 89 ), ( 95 ), ( 158 )].

The idea of using the Cauchy theorem and Cauchy formula for inverting the exponential X-ray transform was initially used in papers [( 44 ), ( 117 ), ( 121 )], and developed in the presented form in [( 89 ), ( 95 )]. This approach leads to a wide variety of inversion formulas [( 95 ), ( 158 )], including cases when the attenuation μ depends on the direction ω [( 95 ), ( 158 )] and when imperfect detector collimation is taken into account [( 89 )]. Invertibility of the exponential X-ray transform was also studied analytically and numerically by different methods in [( 63 )−( 65 )]. One might also be interested in related inversion formulas available for the cases of half-view [( 127 )], three-dimensional limited data acquisition [( 97 )], and imperfect collimation [( 89 )].

Another interesting relation with complex analysis is the following. The possibility of pasting disks to totally real surfaces in R 2 is a popular topic of complex analysis (e.g. [( 55 )]). In particular, the question of rigidity of such a pasted disk has been studied. Rigidity means absence of parametric deformation of the pasted analytic disks. The natural way the disk (A. 14) pasted to the surface S φ appears, suggests that one might expect its rigidity. This happens to be true [( 112 )]. Moreover, a general totally real surface in C 2 in a neighborhood of the boundary of a pasted analytic disk can be reduced to a canonical form resembling the definition of S μ, which in turn leads to the possibility of counting the number of parameters in deformations of this disk [( 112 )]. This provides a different approach to the known results [( 55 )] concerning this problem.

A.3.2.2 Attenuated transform

We follow in this section the technique introduced in [( 8 ), ( 25 )] with an additional development from [( 91 )]. Let as before f(x) be the source's intensity in a planar domain ω, u(x, φ) the density of particles at x moving in the direction of νφ = (cosφ, sinφ) and μ(x) the linear attenuation coefficient. When scattering is absent, the transport equation (A.10) can be rewritten as

$Display mathematics$

(p.701) Using the complex coordinate z = x1 − ix2 one gets

(A.15)
$Display mathematics$
or when the attenuation is zero,
(A.16)
$Display mathematics$

The boundary data u(x, φ)|x∈∂ω known. Let us consider the spaces L 2(S) and H 2L 2(S) (the Hardy space) on the circle and the orthogonal projector P : L 2(S) → H 2. We also denote by u k = u k (x) the Fourier coefficients of u(x, φ). Let R be the right shift in H 2:(Ru)k = u k−1 and L = R* be the left shift. Then we easily obtain from (A. 16) that

$Display mathematics$
or
(A.17)
$Display mathematics$
and
(A.18)
$Display mathematics$

Let υ = P υ = υe ⊕ υ0, where υe = (υ0, υ2,…) and υ0 = (υ1, υ3,…). Then ( A.18 ) can be rewritten as

(A.19)
$Display mathematics$
which splits into independent equations for υe and υ0:
$Display mathematics$

The boundary values of both vectors υe and υ0 at ∂Ω are known. According to ( A.17 ), only the odd part υ0 plays a role in the recovery of f(x).

Let us introduce a generalized ∂¯ operator:

$Display mathematics$

The equation ( A.19 ) is a generalized ∂¯ equation for the function υ If one could prove uniqueness of its solution under given boundary conditions, this would immediately imply uniqueness of recovery of f, and hence uniqueness for the attenuated Radon transform. This would follow, for instance, from a Cauchy-type formula for such “analytic” vector-valued functions (if A = L 2, such functions are called A-analytic in [( 8 ),( 25 )]). This is exactly what will be outlined in the next subsection.

(p.702) A.3.2.3 Cauchy-type formulas and inversion of Radon transform

Let ω be a planar domain and u be an H 2-valued function. Then analogs of the standard considerations lead to a Cauchy–Green-type formula:

Theorem A.6 If uC l (ω¯,H2) and z 0 ∈ ω then

(A.20)
$Display mathematics$

Corollary A.7 If Du = 0 (i.e. u is A-analytic), then

(A.21)
$Display mathematics$
This gives an alternative reconstruction procedure for the standard Radon transform:
1. (1) Find the boundary values for u e and u 0.

2. (2) Find u e and u 0 using ( A.21 ).

3. (3) Find f from υ0 using ( A.17 ).

Remark 1. As one can read in [( 8 ), ( 25 )], the resolvent (zz 0L 2 (z¯ − z¯0)) involved in the above formulas is formally not defined, since the point λ = (zz 0) / (z¯ − z¯0) always belongs to the spectrum of the operator L 2. It can, however, be extended from the exterior of the disk |λ| ≤ 1 as a strong limit and hence defined as an unbounded operator.

Remark 2. It is interesting to recognize the standard parts of the Radon inversion formula (Hilbert transform, differentiation, and back projection) in the procedure above. The continuation of the resolvent up to the spectrum gives the Hilbert transform, the integration in the Cauchy formula ( A.21 ) corresponds to the back projection, and the differentiation is performed in ( A.17 ).

A.3.2.4 Inversion of the attenuated transform

We have dealt so far with the case of zero attenuation. Let us outline now how one can resolve the situation with the attenuated transform, i.e. when the transport equation has the form

$Display mathematics$
Availability of a Cauchy-type formula would have the same effect of producing the uniqueness result and inversion formulas as above. On a formal level, if we had a “function” h such that
(A.22)
$Display mathematics$
(p.703) then Leibnitz's rule would imply that D(hu) = hf. This would mean that the attenuation is essentially eliminated, and the above Cauchy-type procedure is applicable to recover hf and hence f. It is easy to find one such function:
$Display mathematics$
This, however, is not as simple as it sounds. Indeed, h(x, ω) for a fixed x acts as the Toeplitz operator φ → P(hφ)) in the Hardy space, and hence we are dealing with an operator-valued function of x. On the other hand, the differential operator D = ∂¯ + L 2∂ also has operator coefficients. In particular, D and h do not commute in general. This means that Leibnitz's rule does not work. If, however, one could find a solution of ( A.22 ) that for each x has only nonpositive frequencies in its Fourier series with respect to ω then h and D would commute and the procedure would work. This was implemented (in a not entirely explicit form) in [( 8 )], which produced the long-awaited proof of uniqueness for the attenuated Radon transform. However, applying more effort one can show that the procedure can be made explicit [( 91 )], which leads to a variety of inversion formulas, including the one obtained in a different way in [( 128 )] (see also a different derivation in [( 125 )]):
(A.23)
$Display mathematics$
where
$Display mathematics$
This formula was implemented numerically in [( 58 ), ( 96 ), ( 98 ), ( 125 )].

We start with the case of the exponential transform, which is simpler than the one of the general attenuated transform and also leads to interesting analysis.

A.3.3.1 Exponential transform

Probably the first appearance of the range conditions for this transform was in the papers [( 12 ), ( 167 )] devoted to its inversion. We consider a compactly supported smooth function f(x) in the plane and its exponential X-ray (Radon) transform with attenuation μ:

$Display mathematics$
Here ω = (cosφ, sinφ) is a unit vector and ω = (− sinφ, cosφ). Now we expand g(s, ω, μ), into the Fourier series with respect to the angle φ:
$Display mathematics$

(p.704) We consider the Fourier transform ∶ι l(σ,μ) of gl (s,μ) with respect to s. It was observed in [( 12 ), ( 167 )] that the function

$Display mathematics$
is even with respect to σ for Any ι ∈ Z. Papers [( 93 )] and [( 94 )] were devoted to finding tho complete set of range conditions. It was shown In particular that the set of conditions described above is complete. Namely, the following theorem holds:

Theorem A.8 A function g(s, ω) can be represented as Rμ f for some $f ∈ C 0 ∞ ( R )$ if and only if

1. (i) $f ∈ C 0 ∞ ( R × S 1 )$

2. (ii) function (σ +μ)ιĝ(σ, μ) is even with respect to σ for any ι ∈ Z.

An interesting corollary of the range condition provided in this theorem is that the function ĝι(σ, μ) has a zero of order |ι| at the point μ × sgn ι. Existence of these zeros, however, is obviously not equivalent to the whole set of conditions. One can check directly that ( A.30 ) in the case of constant attenuation leads to the above-mentioned root conditions. This shows in particular that the range conditions ( A.30 ) for the attenuated transform that will be discussed later are incomplete [( 93 )].

The range conditions of Theorem A.8 do not have the usual momentum form as they do in the case of the Radon transform. A momentum-type set of conditions was also found in [( 93 )] and [( 94 )].

Theorem A.9 A function g(s, ω) can be represented as Rμf for some $f ∈ C 0 ∞ ( R 2 )$ if and only if

1. (i) $f ∈ C 0 ∞ ( R × S 1 )$

2. (ii) the following identity is satisfied for any odd natural n:

(A.24)
$Display mathematics$
Here i is the imaginary unit, ω = (cosφ, sinφ), and ο denotes composition of differential operators.

The condition ( A.24 ) is not very intuitive. One way we can try to understand it is the following. Let us assume that g is representable as R μf and plug R μ f instead of g into ( A.24 ). In the case of the standard Radon transform, this would lead to a condition that is obviously satisfied. In other words, for the Radon transform, necessity of the momentum conditions (as soon as they are formulated) is simple to observe, and only sufficiency requires a proof. However, even necessity of ( A.24 ) is not obvious. A direct calculation shows that necessity of this condition is equivalent to the following series of identities for the function sinφ.

(p.705) For any odd natural n

(A.25)
$Display mathematics$
The reader might want to try to establish these identities directly. This was done in [( 93 )], where it was also shown that an analogous identity can be written in any commutative differential algebra. A similar identity for all natural numbers n holds for exp (x). Further discussion and reformulation of these conditions can be found in [( 131 )].

A better understanding of the meaning of the range conditions for the exponential Radon transform can be achieved using complex analysis. This was done in the series of papers [( 3 ), ( 4 ), ( 131 ), ( 132 )]. It was shown that the range description problem essentially reduces to Bernstein-Hartogs’ type of theorems on extension of separately analytic functions (see Chapter 5 and further details later on in this chapter).

Consider again formula ( A. 12 ):

$Display mathematics$
If now we have a function g(s, ω) which we suspect of being in the range of the operator R μ, then the function defined on S μ as
(A.26)
$Display mathematics$
must be extendable from S μ to the whole C 2 as an entire function of the Paley-Wiener class. If this is correct, then the extension provides the Fourier transform of a compactly supported function f such that g = R μ f. So, the question arises of extendability of the function defined on the left side of ( A.26 ) to an entire function of the Paley-Wiener class. Such an extendability condition would provide a range description for the operator R μ. In fact, such a condition can be easily found and has a simple geometric meaning.

Returning to the projection-slice formula ( A.12 ), we notice that although we initially used the real values of σ only, we can in fact use all complex values of σ provided that the function f is compactly supported. This leads to the determination of values of on a three-dimensional variety in C 2. Let us denote this variety by ∑μ:

$Display mathematics$
This variety is smooth at all points where Re σ ≠ 0. At these points it has a complex tangent line. Thus, a boundary ∂¯ equation must be satisfied if we want our function to be extendable. One can check, however, that this condition is automatically satisfied due to analyticity of ĝ(σ, ω) with respect to σ. So, where (p.706) is the extendability condition hidden? It turns out that the right place to look for it is at imaginary values of σ. If σ = ip for a real p, then the definition of the function F looks as follows:
(A.27)
$Display mathematics$
The line ppμ +μω is tangent at the point μω to the circle of radius μ in the plane. Changing ω, we rotate this tangent line, spanning the whole exterior of the disk of radius μ. Let us also notice that the function F is entire with respect to p (i.e. along each tangent line). Since there are two such tangent lines passing through each point in the exterior of this disk, the formula ( A.27 ) might provide a self-contradictory definition of the function F. Let us impose a condition that guarantees that this does not happen. That is, we assume that
(A.28)
$Display mathematics$
whenever
$Display mathematics$
It turns out that this simple geometric condition is sufficient. That is, the following theorem holds.

Theorem A. 10 Let F be a function defined in the exterior of the disk of a positive radius in the plane. If this function is entire along each tangent line to the disk, it is extendable to an entire function of two variables.

This theorem was proven (although not formulated explicitly in the present form) in [( 3 )]. The proof uses an analytic trick, whose role has not yet been well understood. For instance, the question remains whether one can replace tangent lines to a circle by tangent lines to a more general convex algebraic curve. This is certainly a separate analyticity theorem of Bernstein's kind [( 5 ), ( 18 )], since in a neighborhood of any point in the exterior of the disk the pairs of tangent lines provide a local coordinate system, for which we have separate analyticity. A proof that uses purely complex analytic tools was provided in [( 132 )]. It was shown in [( 136 )] and [( 131 )] that the strange moment conditions ( A.24 ) can be understood as the infinitesimal form of ( A.28 ) at the boundary of the disk.

The last theorem leads to the following range description [( 3 )]:

Theorem A.11 Let g(s, ω) be a function on R × S 1, and μ be a positive real. Then g = R μ f for some function $f ∈ C 0 ∞ ( R 2 )$ if and only if $g ∈ C 0 ∞ ( R × S 1 ) ,$ and the Fourier transform ĝ(ξ,φ) with respect to the s variable of g(s, ω(φ)) satisfies the following condition for all real σ:

(A.29)
$Display mathematics$

Here ω(φ) = (cosφ, sinφ).

(p.707) The condition ( A.29 ) is a simple restatement of ( A.28 ).

A generalization of this result to the multidimensional case is provided in [( 4 )]. A nice discussion can be found in [( 131 )] and [( 132 )].

Let us mention briefly some applications of these range conditions. It was shown in [( 137 )] how the conditions can be effectively used in SPECT in order to detect and correct some data errors arising from hardware imperfection. The paper [( 111 )] uses the range conditions in order to treat numerically incomplete data problems of SPECT. It is suggested in paper [( 32 )] that the range conditions might be used in some problems of radiation treatment planning (see also [( 31 )] and [( 88 )]). We will not go into details of these results.

A.3.3.2 Attenuated transform

The question about the range descriptions for the attenuated Radon transform turns out to be a tricky one (see also Section A.2 ). In fact, one could think that there is no chance of getting any explicit range conditions for the attenuated X-ray transform. Despite this natural belief, Natterer [( 124 )] found the following infinite set of such range conditions. Let us introduce some notation first. We denote by H the Hilbert transform on the line

$Display mathematics$
where υ. p. denotes the principal value integral (the limit as ∈ → 0 of the integral over R \ (−∈, ∈)).

Theorem A. 12 [( 124 )] Let f and μ belong to the Schwartz space S(R 2). Then, for k > m ≥ 0 integers, we have

(A.30)
$Display mathematics$
where ω = (cosφ,sinφ), I is the identity operator, H is the Hilbert transform, andis the Radon transform of μ.

As has already been explained, the set of conditions provided in this theorem is incomplete. Nevertheless, it is of practical importance, since it can be used for the simultaneous recovery of the sources f(x) and attenuation μ (x) (see [( 121 )−( 123 )] and further discussion below).

The question had remained on how to write down the missing conditions, until the very recent publication [( 130 )], where a complete set of conditions was obtained.

It can be shown [( 91 )] that the approach of [( 8 )] provides an explanation and natural derivation of ( A.30 ). In order to explain this, let us start with the case of zero attenuation (i.e. standard Radon transform) and look at it in terms of .A-analytic functions (see Section ( A.3.2.2 )). In the notation of that section, the role of the Radon data is played by a function u on ∂ω with values in the Hardy (p.708) space. It belongs to the range iff it is 4-analytically extendable to ω. Due to an analog of Sokhotsky's formula for such functions, it is easy to show that this happens iff

(A.31)
$Display mathematics$
for any z 0 ∉ Ω¯. These conditions, however, do not resemble the standard momentum ones. How do we get those? In standard complex analysis, the boundary values of an analytic function as a measure on ∂Ω must be orthogonal to the boundary values of any function analytic in Ω and continuous up to the boundary. Any base of such functions would do. For instance, powers of z give
$Display mathematics$
Similarly for A-analytic functions: ( A.31 ) means orthogonality on the boundary to the operator functions
$Display mathematics$
Is there any other convenient basis? What are the analogs of the powers z n here? The operator-valued function H(z) = Rz − L¯z satisfies DH = 0. It is an A-analytic analog of z. Due to the operator nature of H and D, the product Hu is not necessarily A-analytic. Namely,
$Display mathematics$
This expression is not necessarily equal to zero. However, the commutator [L 2, R] = −P 1 L, where P l is the projector onto the first l terms of Fourier expansion. In particular, L[L 2, R] = 0. This implies that
$Display mathematics$
and analogously
$Display mathematics$
Theorem A.13 [( 91 )] These conditions coincide (modulo evenness) with the standard momentum conditions.

How can one now include attenuation in this scheme? Using an appropriate “correction” function h as in Section ( A.3.2.2 ) and then applying the above conditions to hu one obtains the range conditions. In particular, one can obtain ( A.30 ).

(p.709) Theorem A.14 [( 91 )] Orthogonality to the suitably understood functions H k h (i.e. analogs of orthogonality to z k ) gives Natterer's range conditions ( A.30 ).

We would like to mention again that a full set of range conditions was obtained in [( 130 )].

As we have discussed in the beginning of this section, simultaneous recovery of the source's density f(x) and of the attenuation μ (x) is an important applied issue. At the first glance this problem might look hopeless. Indeed, in the plane case the data g = Tμf is a single function of two variables, while we are trying to recover two functions f(x) and μ (x) of two variables. However, this counting of functions and variables would be persuasive only if the operator Tμ were close to an invertible one (for instance, if it were a Fredholm operator). In fact, the study of the range conditions shows that the operator Tμ has an infinite dimensional cokernel. This means that when μ changes, the range could in principle rotate in such a way that for two different values of the attenuation the two ranges would have zero (or a “very small”) intersection. If this were true, then both f and μ would be recoverable or “almost recoverable.” For instance, some additional a priori information could help. We will describe now some results that show that something like this does happen, although the understanding of the situation is still very superficial.

Let us start with the simplest case of the exponential X-ray transform. The first attempt to recover both f and μ from Rμf was made in [( 75 )]. This problem was resolved in [( 162 )] (see also [( 163 )]). Consider the range conditions in the form of evenness of the function

$Display mathematics$
This implies the equality
$Display mathematics$
whenever ĝl(σ, μ) is not identically zero. This means that if there is a nonzero harmonic ĝl(σ, μ) for some l ≠ 0, then the rational function
$Display mathematics$
is uniquely determined. This, in turn, determines μ uniquely as the pole of this function. The only exception arises when ĝ i(σ, μ) = 0 for all l ≠ 0. In this case, however, the function g = Rμ f, and hence f itself, is radial. For radial functions no information about the value of μ can be extracted from Rμf. Geometrically this means that the intersection of ranges of R μ for different values of μ consists of radial functions. In other words, we established the following theorem.

(p.710) Theorem A.15 [( 162 ), ( 163 )] Simultaneous recovery of a compactly supported function f and constant attenuation μ from the values of R μ f is possible if and only if f is not radial.

In fact, in the non-radial case one can provide an algorithm of finding the attenuation μ (see [( 162 )] and [( 163 )]).

An alternative approach to recovery of attenuation was described in [( 7 )].

Recovery of a variable attenuation is definitely much more difficult. As in the exponential case, the range theorems are used to this end. Theorem A.12 has been used under some additional restrictions for the sources distribution in order to accomplish this task [( 121 )−( 123 )]. Papers [( 8 ), ( 128 )−( 130 )] contain some additional indications on how their techniques and results might be employed for that purpose.

Electrical impedance tomography (EIT) is a very promising, inexpensive, and technologically simple (at least in comparison with X-ray, SPECT, PET, and MRI scans) method of medical diagnostics and of industrial nondestructive testing (see, for instance, [( 9 ), ( 10 ), ( 16 ), ( 17 ), ( 23 ), ( 28 ), ( 29 ), ( 47 ), ( 54 ), ( 77 ), ( 87 ), ( 104 ), ( 150 ), ( 154 ), ( 155 ), ( 164 ), ( 165 )] and references therein). Here is the idea of EIT: one places electrodes around a body, creates some known currents through the electrodes, and then measures the corresponding boundary voltage drops. All of the data is collected on the boundary. After many such measurements one wants to recover the distribution of the electric conductivity inside the body. The information about the electric conductivity is very important for medical diagnostics; it is also vital for some electrical procedures, such as defibrillation; it could also provide a cheap and reliable nondestructive evaluation technology. There are several strong groups in industry and academe that work on the practical and theoretical aspects of this problem using various techniques. We will consider here only one approach, which is related to the integral geometry problems discussed in the book. The reader can refer to the papers quoted above that address other methods that do not use the Radon transform.

Let us describe first the setup of the inverse conductivity problem, which is the mathematical formulation of EIT. Let U ⊂ R n be a domain with sufficiently smooth boundary γ. An unknown function β(x) (the conductivity) must be recovered from the following data. Given a function ψ on γ (the current) one solves the Neumann boundary value problem

$Display mathematics$
where ν is the unit outer normal vector on γ. Then one measures the boundary value φ = u|γ (the potential). All the pairs (ψ, φ) are assumed to be accessible.

(p.711) In other words, the so-called Dirichlet-to-Neumann operator λβ : φ → ψ is known. One needs to recover the conductivity β from this data. The problems is obviously nonlinear.

The inverse conductivity problem is a hard nut to crack both analytically and numerically. The main questions are about uniqueness of determination of the conductivity, stability of the reconstruction, and inversion algorithms. Although the problem of uniqueness can be considered as principally resolved (see [( 78 ), ( 119 ), ( 166 ), ( 170 )], and references therein), the other ones are still under thorough investigation. The general understanding is that the problem is highly unstable, so there is no hope of achieving the quality of reconstruction known to other common tomographic techniques, at least without involving additional information about the image to be reconstructed.

The first practical algorithm of Barber and Brown [( 9 )−( 11 )] deals with the linearized problem, which is still highly unstable. A thorough investigation of this algorithm started by Santosa and Vogelius [( 155 )] lead in the works of Berenstein and Casadio Tarabusi [( 16 ), ( 17 )] to the understanding that the linearized two-dimensional problem can be treated by means of hyperbolic geometry (see Section ( 7.3 )). Consider the two-dimensional case when U is the unit disk. It is well known (see for instance [( 13 ), ( 69 )]) that the unit disk serves as the Poincare model of the hyperbolic plane H 2. Here are the indications why the hyperbolic geometry might play some role in the (at least linearized) inverse conductivity problem. First of all, the Laplace operator that arises in the linearized problem is invariant with respect to the group of Mobius transformations. Another indication is that if one creates a dipole current through a point on the boundary of U, then the equipotential lines and the current lines form families of geodesies and horocycles in H 2. Following the analysis done in [( 155 )] of the algorithm suggested in [( 9 ), ( 10 )], Berenstein and Casadio Tarabusi [( 16 ), ( 17 )] discovered that in fact for n = 2 the linearized inverse conductivity problem can be reduced to the following integral geometry problem on H 2: the available data enables one to find the function

$Display mathematics$
where R G is the geodesic Radon transform on H 2, A is an explicitly described radial function on H 2
$Display mathematics$
and the asterisk * denotes the (non-euclidean) convolution on H 2. Now one can hope to use the known methods of harmonic analysis on H 2 to recover β. That is, one needs to be able to invert the geodesic Radon transform and to deconvolve. Necessary details of harmonic analysis on H2 can be found in [( 68 ), ( 70 )]. Several different inversion formulas for the geodesic Radon transform are developed in [( 15 ), ( 69 ), ( 106 )]. The paper [( 52 )] provides a unified approach to such formulas. The reader might also be interested in the related papers [( 99 ), ( 101 )]. The formula obtained in [( 106 )] was numerically implemented in [( 45 )]. Another part of inversion would be deconvolution. Its numerical implementation can be done by using the (p.712) Fourier transform on the hyperbolic plane described in [( 68 )] and its inversion. A different approach to the Fourier transform is given in Chapter ( 7 ). The Fourier transform is defined as follows:
$Display mathematics$
where b ∈ ∂H 2, λ ∈ R, 〈z,b〉 is the (signed) hyperbolic distance from the origin to the horocycle passing through the points z and b, and dm(z) is the invariant measure on H 2. The inverse Fourier transform is
$Display mathematics$
These operators were numerically implemented in [( 45 )]. (By an editorial error, all pictures have been omitted in [( 45 )]. They can be found at the URL http://www.math.tamu.edu/kuchment/hypnum.pdf.) Although the problem of numerical implementation of all these transforms looks familiar and similar to the euclidean case, it is in fact much more complex. In particular, even the problem of sampling is not a trivial one. One natural way of sampling would be the following. One chooses a discontinuous group of motions of the hyperbolic disk and uses the orbit of the origin as the grid. Choosing a small mesh, one samples with appropriate accuracy. Unfortunately, this idea does not work, due to the known rigidity of such lattices. Namely, the size of the pixel (i.e. the fundamental domain) for any discrete group acting on H 2 cannot be made arbitrarily small [( 13 )]. The way of overcoming this difficulty chosen in [( 45 )] was to create a combination of euclidean and hyperbolic grids. A “superlattice” was created using a discrete group of hyperbolic motions, and then a finer euclidean grid was chosen in each of the fundamental domains.

It seemed at first that the hyperbolic geometry approach would not work in dimensions higher than two, due to lack of hyperbolic invariance of the governing equations. It was shown, however, in [( 46 )] that a combination of euclidean and hyperbolic integral geometries does the trick.

Notes:

(1) The second author was partially supported by NSF grants 9877155 and 0200788.