7 Mathematics of thermoacoustic tomography
Peter Kuchment∗and Leonid Kunyansky†
November 2, 2018
Abstract
The paper presents a survey of mathematical problems, techniques,
and challenges arising in the Thermoacoustic (also called Photoacous-
tic or Optoacoustic) Tomography.
1 Introduction
Computerized tomography has had a huge impact on medical diagnostics.
Numerous methods of tomographic medical imaging have been developed
and are being developed (e.g., the “standard” X-ray, single-photon emission,
positron emission, ultrasound, magnetic resonance, electrical impedance, op-
tical) [55, 59, 75, 76, 77]. The designers of these modalities strive to in-
crease the image resolution and contrast, and at the same time to reduce
the costs and negative health effects of these techniques. However, these
goals are usually rather contradictory. For instance, some cheap and safe
methods with good contrast (like optical or electrical impedance tomogra-
phy) suffer from low resolution, while some high resolution methods (such
as ultrasound imaging) often do not provide good contrast. Recently re-
searchers have been developing novel hybrid methods that combine different
physical types of signals, in hope to alleviate the deficiencies of each of the
types, while taking advantage of their strengths. The most successful exam-
ple of such a combination is the Thermoacoustic Tomography (TAT)
∗Mathematics Department, Texas A& M University, College Station, TX 77843-3368,
USA. kuchment@math.tamu.edu
†Mathematics Department, University of Arizona, Tucson, AZ 77843-3368, USA.
leonk@math.arizona.edu
http://arxiv.org/abs/0704.0286v2
1[62]. Albeit not being a common feature in clinics yet, TAT scanners are
actively researched, developed and already manufactured, for instance by
OptoSonics, Inc. (http://www.optosonics.com/), founded by the pioneer of
TAT R. Kruger.
After a substantial effort, major breakthroughs have been achieved in the
last couple of years in the mathematical modeling of TAT. The aim of this
article is to survey this recent progress and to describe the relevant models,
mathematical problems, and reconstruction procedures arising in TAT, and
to provide references to numerous research publications on this topic.
The main thrust of this text is toward mathematical methods; consid-
erations of the text length, as well as authors’ background do not let us
discuss in any detail industrial and physical set-ups and parameters of the
TAT technique, and limitations of the corresponding mathematical models.
Fortunately, the excellent recent surveys by M. Xu and L.-H. V. Wang [117]
and by A. A. Oraevsky and A. A. Karabutov [87, 88] accomplish all of these
tasks, and thus the reader is advised to consult with them for all such is-
sues (see also the recent textbook [113]). On the other hand, in spite of
the significant recent progress in mathematics of TAT, there is no compre-
hensive survey text addressing in details the relevant mathematical issues,
although the surveys [88, 117] do mention some mathematical reconstruction
techniques.
The structure of the paper is a follows: Section 2 contains a brief descrip-
tion of the TAT procedure. The next Section 3 provides the mathematical
formulation of the TAT problem. In general, it is formulated as an inverse
problem for the wave equation. However, in the case of the constant sound
speed, it can be also described in terms of a spherical mean operator (a
spherical analog of the Radon transform). The section also contains the list
of natural questions to be addressed concerning this model. These issues are
addressed then one by one in the following sections. In particular, Section
4 discusses uniqueness of reconstruction, i.e. the question of whether the
data collected in TAT is sufficient for recovery of the information of interest.
Albeit, for all practical purposes this issue is resolved in Corollary 2, we pro-
vide an additional discussion of unresolved uniqueness problems, which are
probably of more academic interest. Section 5 addresses inversion formulas
and algorithms. In Section 6 effects of having only partial data are discussed.
1TAT is also called Photoacoustic (PAT) or Optoacoustic (OAT) Tomography and is
sometimes abbreviated as TCT, which stands for Thermoacoustic Computed Tomography
http://www.optosonics.com/
Section 7 contains results concerning the so called range conditions, i.e. the
conditions that all ideal data must satisfy. Section 8 provides additional
remarks and discussions of the issues raised in the previous sections. The
paper ends with an Acknowledgments section and bibliography. Concerning
the latter, we need to mention that the engineering and biomedical literature
on TAT is rather vast and no attempt has been made in this text to create
a comprehensive bibliography of the topic from the engineering prospective.
The references in [87, 88, 109, 112, 117] to a large extent fill this gap. The
authors, however, have tried to present a sufficiently complete review of the
existing literature on mathematics of TAT.
2 Thermoacoustic tomography
In TAT, a short duration EM pulse is sent through a biological object (e.g.,
woman’s breast in mammography) with the aim of triggering a thermoa-
coustic response in the tissue. As it is explained in [117], the radiofrequency
(RF) and the visible light frequency ranges are currently considered to be the
most suitable for this purpose. Since mathematics works exactly the same
way in both of these frequency ranges, we will not make such distinction and
will be talking about just “an EM pulse”. E.g., in Figure 1 a microwave
pulse is assumed. In most cases the pulse is spatially wide, so that the
Figure 1: The TAT procedure.
whole object is more or less uniformly irradiated. Some part of EM energy
is absorbed throughout the object. The amount of energy absorbed at a
location strongly depends on local biological properties of the cells. Oxygen
saturation, concentration of hemoglobin, density of the microvascular net-
work (angiogenesis), ionic conductivity, and water content are among the
parameters that influence the absorption strongly [117]. Thus, if the energy
absorption distribution function f(x) were known, it would provide a great
diagnostic tool. For instance, it could be useful for detecting cancerous cells
that absorb several times more energy in the RF range than the healthy
ones [62, 88, 115, 117]. However, as an imaging tool neither RF waves, nor
visual light alone would provide acceptable resolution. In the RF case, this
is due to the long wave length. One can use shorter microwaves, but this
will be at the expense of the penetration depth. In the optical region, the
problem is with the multiple scattering of light. So, a different mechanism,
the so called Photoacoustic Effect [46, 107, 113, 117], is used to image f(x).
Namely, the EM energy absorption results in thermoelastic expansion and
thus in a pressure wave p(x, t) (an ultrasound signal) that can be measured
by transducers placed around the object. Now one can attempt to recover the
function f(x) (the image) from the measured data p(x, t). Such a measuring
scheme, utilizing two types of waves, brings about the high resolution of the
ultrasound diagnostics and the high contrast of EM waves. It overcomes the
adverse effect of the low contrast of ultrasound with respect to soft tissue.
In fact, such a low contrast is a good thing here, allowing one to assume in
the first approximation that the sound speed is constant. This often used
approximation is not always appropriate, but it is the most studied case at
the moment. Later on in this text we will describe some initial considerations
of the variable sound speed case, following [4, 57].
For this TAT method (and in particular, for the mathematical model
described below) to work, several conditions must be met. For instance, the
time duration of the EM pulse must be shorter than the time it takes the
sound wave to traverse the smallest feature that needs to be reconstructed.
The ultrasound detector must be able to resolve the time scale of the duration
of the EM pulse. On the other hand, the transducer must be also able to
detect much lower frequencies. Thus, one needs to have extra-wide-band
transducers, and these are currently available. One can find the technical
discussion of all these issues, for instance, in [88, 117]. In this text we will
assume that all these conditions are met and thus the mathematical models
described are applicable.
In the next section we present a mathematical description of the relation
between f(x) and p(x, t) (similar mathematical problems arise in sonar [73]
and radar [81] imaging, as well as in geophysics [27]).
3 Mathematical model of TAT:
wave equation and the spherical mean trans-
3.1 The wave equation model
We assume that the ultrasound speed at location x is equal to c(x). Then,
modulo some constant coefficients that we will assume all to be equal to 1,
the pressure wave p(x, t) satisfies the following problem for the standard wave
equation [28, 107, 115]:
ptt = c
2(x)∆xp, t ≥ 0, x ∈ R
p(x, 0) = f(x),
pt(x, 0) = 0
The goal is to find, using the data measured by transducers, the initial value
f(x) at t = 0 of the solution p(x, t).
In order to formalize what data is in fact measured, one needs to specify
what kind of transducers is used, as well as the geometry of the measurement.
By the geometry of the measurement we mean the distribution of locations
of transducers used to collect the data.
We briefly describe here the commonly considered measurement proce-
dure, which uses point detectors. Line and planar detectors have also been
suggested (see Section 8.1.1). It is too early to judge which one of them
will become most successful, but the one using point transducers has been
more thoroughly studied mathematically and experimentally, and thus will
be mostly addressed in this article. In this case, the transducers are assumed
to be point-like, i.e. of sufficiently small dimension. A transducer at time
t measures the average pressure over its surface at this time, which for the
small size of the transducer can be assumed to be just the value of p(y, t) at
the location y of the transducer. Dimension count shows immediately that in
order to have enough data for reconstruction of the function f(x), one needs
to collect data from the transducers’ locations y running over a surface S in
3. Thus, the data at the experimentalist’s disposal is the function g(y, t)
that coincides with the restriction of p(x, t) to the set of points y ∈ S.
Taking into account that the measurements produce the values g(y, t) of
the pressure p(x, t) of (1) on S × R+, the set of equations (1) extends to
become
ptt = c
2(x)∆xp, t ≥ 0, x ∈ R
p(x, 0) = f(x),
pt(x, 0) = 0
p(y, t) = g(y, t), y ∈ S × R+
Figure 2: An illustration to (2).
The problem now reduces to finding the initial value f(x) in (2) from the
knowledge of the lateral data g(x, t) (see Figure 3.1). A person familiar with
PDEs might suspect first that there is something wrong with this problem,
since we seem to have insufficient data for the recovery of the solution of the
wave equation in a cylinder from the lateral values alone. This, however,
is an illusion, since in fact there is a significant additional restriction: the
solution holds in the whole space, not just inside the cylinder S × R+. We
will see soon that in most cases, the data is sufficient for recovery of f(x).
3.2 Spherical mean model
We now introduce an alternative formulation of the problem that works in
the constant speed case only. We will assume that the units are chosen in
such a way that c(x) = 1. The known Poisson-Kirchhoff formula [25, Ch.
VI, Section 13.2, Formula (15)] for the solution of (1) gives
p(x, t) = c
(t(Rf)(x, t)) , (3)
where
(Rf)(x, r) =
|y|=1
f(x+ ry)dA(y) (4)
is the spherical mean operator applied to the function f(x), and dA is the
normalized area element on the unit sphere in R3. Hence, knowledge of the
function g(x, t) for x ∈ S and all t ≥ 0 essentially means knowledge of the
spherical mean Rf(x, t) at all points (x, t) ∈ S × R+. One thus is lead to
studying the spherical mean operator R : f → Rf and in particular its
restriction RS to the points x ∈ S only (these are the points where we place
transducers):
RSf(x, t) =
|y|=1
f(x+ ty)dA(y), x ∈ S, t ≥ 0. (5)
This explains why, in many works on TAT, the spherical mean operator has
been the model of choice. Albeit the (unrestricted) spherical mean operator
has been studied rather intensively and for a long time (e.g., [17, 25, 58]),
its version RS with the centers restricted to a subset S appears to have been
studied since early 1990s only [1]-[14], [16, 26, 30, 31, 33, 32, 34, 35, 36, 39,
41, 63, 64, 68, 69, 71, 72, 73, 77, 80, 82, 83, 89, 90, 91, 95, 96, 97, 104, 121]
and offers quite a few new and often hard questions.
In what follows, we will alternate between these two (PDE and integral
geometry) interpretations of the TAT model, since each of them has its own
advantages.
3.3 Main mathematical problems of TAT
We now formulate the typical list of problems one would like to address in
order to implement the TAT reconstruction.
1. For which sets S ∈ R3 is the data collected by transducers placed along
S sufficient for unique reconstruction of f? In terms of the spheri-
cal mean operator, the question is whether RS has zero kernel on an
appropriate class of functions, say continuous with compact supports.
2. If the data collected from S is sufficient, what are inversion formulas
and algorithms?
3. How stable is the inversion?
4. What happens if the data is “incomplete”?
5. What is the space of all possible “ideal” data g(t, y) collected on a sur-
face S? Mathematically (and in the constant sound speed case) it is
the question of describing the range of the operator RS in appropriate
function spaces. This question might seem to be unusual (for instance,
to people used to partial differential equations), but in tomography im-
portance of knowing the range of Radon type transforms is well known.
Such information is used to improve inversion algorithms, complete in-
complete data, discover and compensate for certain data errors, etc.
(e.g., [30, 38, 39, 40, 53, 54, 55, 76, 77, 90]).
4 Uniqueness of reconstruction
Many of the problems of interest to TAT can be formulated in any dimension
d, albeit the practical dimensions are only d = 3 and d = 2. We will consider
an arbitrary dimension d whenever we see this suitable.
Let S ⊂ Rd be the set of locations of transducers and f be a compactly
supported function (one can show that for purposes of uniqueness of recon-
struction problem, one can always assume that f is smooth [7]). Does the
absence of the signal on the transducers, i.e. g(t, y) = 0 for all t and y in
S, imply that f = 0? If the answer is a “yes,” we call S - a uniqueness
set, otherwise a non-uniqueness set. In other words, in terms of TAT, the
uniqueness sets are those that distributing transducers along them provides
enough data for unique reconstruction of the function f(x).
In terms of the wave equation, uniqueness sets are the sets of complete
observability, i.e. such that observing the motion on this set only, one gets
enough information to reconstruct the whole oscillation. In terms of the
spherical mean operator, the question is of whether the equality RSf = 0
implies that f = 0.
We will address this problem for the constant sound speed case first.
4.1 Constant speed case
As it has been discussed, the dimension count makes it clear that S must
be (d − 1)-dimensional, i.e. a surface in 3D or a curve in 2D. We will
also see that most of such surfaces are “good”, i.e. are uniqueness ones (or,
in other words, provide enough information for reconstruction). Thus, we
should rather discuss the problem of describing the “bad”, non-uniqueness
sets. The following simple statement is very important and not immediately
obvious.
Lemma 1 [7, 71, 72, 121] Any non-uniqueness set S is a set of zeros of a
(non-trivial) harmonic polynomial. In particular,
1. If there is no non-zero polynomial vanishing on S, then S is a unique-
ness set.
2. If there is no non-zero harmonic function vanishing on S, then S is a
uniqueness set.
The proof of this lemma is very simple. It works under the assumption
of exponential decay of the function f(x), not necessarily of compactness of
its support. It also introduces some polynomials that play significant role in
the whole analysis of the spherical mean operator RS.
Let k ≥ 0 be an integer. Consider the convolution
Qk(x) = |x|
2k ∗ f(x) =
|x− y|2kf(y)dy. (6)
This is clearly a polynomial of degree at most 2k. Rewriting the integral in
polar coordinates centered at x and using radiality of |x − y|, one sees that
Qk(x) is determined if we know the values Rf(x, t) of the spherical mean of
f centered at x:
Qk(x) = cd
t2k+d−1Rf(x, t)dt.
In particular, If RSf ≡ 0, then each polynomial Qk vanishes on S.
Another observation that is easy to justify is that if the function f is
exponentially decaying (e.g., is compactly supported), then if all polynomials
Qk vanish identically, the function itself must be equal to zero. (This is not
necessarily true anymore if f and its derivatives decay only faster than any
power, rather than exponentially.)
Thus, we conclude that if f is not identically equal to zero, then there is at
least one non-zero polynomial Qk. Since, as we discussed, equality RSf = 0
implies that Qk|S = 0, we conclude that S must be algebraic.
Now notice the following simple to verify equality (with a non-zero con-
stant ck):
∆Qk = ckQk−1, (7)
where ∆ is the Laplace operator. This implies that the lowest k non-zero
polynomial Qk is harmonic. Since Qk|S = 0, this proves the lemma.
Consider now the case when S is a closed (hyper-)surface (i.e., the bound-
ary of a bounded domain). Since, as it is well known, there is no non-zero
harmonic function in the domain that would vanish at the boundary (the
spectrum of the Dirichlet Laplace operator is strictly positive), we conclude
that such S is a uniqueness set for harmonic polynomials. Thus, we get the
following important
Corollary 2 [7, 63] Any closed surface is uniqueness set for the spherical
mean Radon transform.
An older alternative proof of this corollary provides an additional insight
into the problem. We thus sketch it here. Let us assume for simplicity
that the dimension d ≥ 3 is odd (even dimensions require a little bit more
work). Suppose that the closed surface S remains stationary (nodal) for the
oscillation described by (1). Since the oscillation is unconstrained and the
initial perturbation is compactly supported, after a finite time, the interior
of S will become stationary. On the other hand, we can think that S is
fixed (since it is not moving anyway). Then, the energy inside S must stay
constant. This is the contradiction that proves the statement of Corollary 2.
We will see in the next Section that the same method works in some cases
of variable sound speed, providing the needed uniqueness of reconstruction
result.
This corollary resolves the uniqueness problems for most practically used
geometries. It fails, however, if f does not decay sufficiently fast (see [3],
where it is shown in which Lp(Rd) classes of functions f(x) closed surfaces
remain uniqueness sets).
It also provides uniqueness for some “limited data” problems. For in-
stance, if S is an open (even tiny) piece of an analytic closed surface Σ, it
suffices. Indeed, if it did not, then it would be a part of an algebraic non-
uniqueness surface. Uniqueness of analytic continuation would show then
that the whole Σ is a non-uniqueness set, which we know to be incorrect.
This result, however, does not say that it would be practical to reconstruct
using observations from a tiny S. We will see later that this would not lead
to a satisfactory reconstructions, due to instabilities.
A geometry sometimes used is the planar one, i.e. detectors are placed
along a plane S (line in the 2D). In this case, there is no uniqueness of
reconstruction when the sound speed is constant. Indeed, if f(x) is odd with
respect to S, then clearly all measured data g(t, y) will vanish. However, it is
well known [25, 58] that functions even with respect to S can be recovered.
What saves the day in TAT is that the object to be imaged is located on one
side of S. Then, extending f(x) as an even function with respect to S, one
can still recover it from the data.
Although, for all practical purposes the uniqueness of reconstruction prob-
lem is essentially resolved by the Corollary 2, the complete understanding of
uniqueness problem has not been achieved yet. Thus, we include below some
known theoretical results and open problems.
4.1.1 Non-uniqueness sets in R2.
In this Section, we follow the results and exposition of [7, 71, 72] in discussing
uniqueness sets in 2D. What are simple examples of non-uniqueness sets? As
we have already mentioned, any line S (or a hyperplane in higher dimensions)
is a non-uniqueness set, since any function f odd with respect to S will clearly
produce no signal: RSf = 0. Analogously, consider a Coxeter system ΣN of
N lines passing through a point and forming equal angles (see Fig. 3).
Figure 3: Coxeter cross ΣN .
Choosing the intersection point as the pole and expanding functions into
Fourier series with respect to the polar angle, it is easy to discover existence
of an infinite dimensional space of functions that are odd with respect to
each of the N lines. Thus, such a cross ΣN is also a non-uniqueness set. Less
obviously, one can use the infinite dimensional freedom just mentioned to
add any finite set Φ of points still preserving non-uniqueness. The following
major and very non-trivial result was conjectured in [71, 72] and proven in
[7]. It shows that there are no other bad sets S besides the ones we have just
discovered:
Theorem 3 A set S ⊂ R2 is a non-uniqueness set for the spherical mean
transform in the space of compactly supported functions, if and only if
S ⊂ ωΣN ∪ Φ,
where ΣN is a Coxeter system of lines, ω is a rigid motion of the plane, and
Φ is a finite set.
A sketch of a rather intricate proof of this result is provided in Section
4.1.2 Higher dimensions
Here we present a believable conjecture of how the result should look like in
higher dimensions.
Conjecture 4 [7]A set S ⊂ Rd is a non-uniqueness set if and only if S ⊂
ωΣ ∪ Φ, where Σ is the surface of zeros of a homogeneous harmonic polyno-
mial, ω is a rigid motion of Rd, and Φ is an algebraic surface of dimension
at most d− 2.
Figure 4: A picture of a 3-dimensional non-uniqueness set.
The progress towards proving this conjecture has been slow, albeit some
partial cases have been treated ([1]-[12]). E.g., in some cases one can prove
that S is a ruled surface (i.e., consists of lines), but proving that these lines
(rules) pass through a common point remains a challenge. It is known,
though, that both the zero sets of homogeneous harmonic polynomials and
algebraic subsets of dimension at most d − 2 are non-uniqueness sets [2, 7],
and thus one should avoid using them as placements of transducers for TAT.
4.1.3 Relations to other areas of analysis
The problem of injectivity of RS has relations to a wide variety of areas of
analysis (see [1, 7] for many examples). In particular, the following interpre-
tation is important:
Theorem 5 [7, 63] The following statements are equivalent:
1. S ⊂ Rd is a non-uniqueness set for the spherical mean operator.
2. S is a nodal set for the wave equation, i.e. there exists a non-zero
compactly supported f such that the solution of the wave propagation
problem 
= ∆u,
u(x, 0) = 0,
ut(x, 0) = f(x)
vanishes on S for any moment of time.
3. S is a nodal set for the heat equation, i.e. there exists a non-zero
compactly supported f such that the solution of the problem
= ∆u,
u(x, 0) = f(x)
vanishes on S for any moment of time.
The interpretation in terms of the wave equation provides important PDE
tools and insights, which have lead to a recent progress [33, 12] (albeit it
has not lead yet to a complete alternative proof of Theorem 3). The rough
idea, originally introduced in [33], is that if S is a nodal set, then it might be
considered as the fixed boundary. In this case, the signals must go around S.
However, in fact, there is no obstacle, so signals can propagate along straight
lines. Thus, in order to avoid discrepancies in arrival times, S must be very
special. One can find details in [33] and in [12].
4.2 Uniqueness in the case of a variable sound speed
It is shown in [35, Theorem 4] that uniqueness of reconstruction also holds in
the case of a smoothly varying (strictly positive) sound speed, if the source
function f(x) is completely surrounded by the observation surface S (in other
words, if there is no US signal coming from outside of S). The proof uses
the celebrated unique continuation result by D. Tataru [108].
One can also establish uniqueness of reconstruction in the case of the
source not necessarily completely surrounded by S. However, here we need
to impose an additional non-trapping condition on the sound speed. We
assume that the sound speed is strictly positive c(x) > c > 0 and such that
c(x)− 1 has compact support, i.e. c(x) = 1 for large x.
Consider the Hamiltonian system in R2nx,ξ with the Hamiltonian H =
c2(x)
|ξ|2: 
x′t =
= c2(x)ξ
ξ′t = −
∇ (c2(x)) |ξ|2
x|t=0 = x0, ξ|t=0 = ξ0.
The solutions of this system are called bicharacteristics and their projections
into Rnx are rays.
We will assume that the following non-trapping condition holds:
all rays (with ξ0 6= 0) tend to infinity when t→ ∞.
Theorem 6 [4] Under the non-trapping conditions formulated above, com-
pactly supported function f(x) is uniquely determined by the data g measured
on S for all times. (No assumption of f being supported inside S is imposed.)
One should mention that ray trapping can occur for some sound speed
profiles. For instance, if c(x) = |x| for some range r1 < |x| < r2, then there
are rays trapped in this spherical shell. We are not sure what happens in
this case to the uniqueness of reconstruction statement of Theorem 6 and
inversion formula of Theorem 7.
5 Reconstruction: formulas and examples
Here we will address the procedures of actual reconstruction of the source
f(x) from the data g(t, y) measured by transducers.
5.1 Constant sound speed
We assume here that the sound speed is constant and normalized to be equal
to 1.
5.1.1 Inversion formulas
Before we move to our case of interest, which is spheres centered on a closed
surface S surrounding the object to be imaged, we briefly refer to related
but somewhat different work. Namely, the problem of recovering functions
from integrals over spheres centered on a (hyper)plane S has attracted a
lot of attention over the years. Albeit, as it has been mentioned before,
there is no uniqueness in this case (functions odd with respect to S are
annihilated), even functions can be recovered. Thus also functions supported
on one side of the plane can be as well, by means of their even extension.
Many explicit inversion formulas and procedures have been obtained for this
situation [16, 26, 31, 39, 41, 60, 77, 80, 89, 90, 101, 102, 103]. We will not
provide any details here, since this acquisition geometry is not very useful.
In particular, this is due to “invisibility” of some parts of the interfaces,
see Section 6, which arises from truncating the plane. The same problem
is encountered with some other unbounded acquisition surfaces, such as a
surface of an “infinitely” long cylinder.
Thus, it is more practical to place transducers along a closed surface
surrounding the object. The simplest surface of this type is a sphere.
5.1.2 Fourier expansion methods
Let us assume that S is the unit sphere in Rn. We would like to reconstruct
a function f(x) supported inside S from the known values of its spherical
integrals g(y, r) with the centers on S:
g(y, r) =
f(y + rω)rn−1dω, y ∈ S.
The first inversion procedures for the case of spherical acquisition were de-
scribed in [82] in 2D and in [83] in 3D. These solutions were obtained by
harmonic decomposition of the measured data and the sought function, and
by equating coefficients of the corresponding Fourier series.
In particular, the 2-D algorithm of [82] is based on the Fourier decompo-
sition of f and g in angular variables:
f(x) =
fk(ρ)e
ikϕ, x = (ρ cos(ϕ), ρ sin(ϕ)) (9)
g(y(θ), r) =
gm(r)e
ikθ, y = (R cos(θ), R sin(θ)).
Following [82] we consider the Hankel transform ĝm,J(λ) of the Fourier coef-
ficients gm(r) (divided by 2πr)
ĝm,J(λ) =
gm(r)J0(λr)dr = H0
gm(r)
. (10)
To simplify the presentation we introduce the convolution GJ(λ, y) of the
sought function with the Bessel function J0(λ|x− y|).
GJ(λ, y) =
f(x)J0(λ|x− y|)dx, (11)
One can notice that ĝm,J(λ) are the Fourier coefficients of GJ(λ, y) in θ:
ĝm,J(λ) =
GJ(λ, y)e
−imθdθ. (12)
Now coefficients fm(ρ) can be recovered from gm(r) by application of the
addition theorem for the Bessel function J0(λ|x− y|):
J0(λ|x− y|) =
Jm(λ|x|)Jm(λ|y|)e
−im(ϕ−θ). (13)
Indeed, by substituting equations (9) and (13) into (11), and (11) into (12)
one obtains
ĝm,J(λ) = 2πJm(λ|R|)
fm(ρ)Jm(λρ)ρdρ = Hm(fm(ρ)),
where Hm is the m-th order Hankel transform. Since the latter transform
is self-invertible, the coefficients fm(ρ) can be recovered by the following
formula
fm(ρ) = Hm
ĝm,J(λ)
Jm(λ|R|)
Jm(λ|R|)
gm(r)
, (14)
which is the main result of [82]. Function f(x) can now be reconstructed by
summing series (9).
Note that the above method requires a division of the Hankel transform
of the measured data by Bessel functions Jm that have infinitely many zeros.
Theoretically, there is no problem; the Hankel transform H0
gm(r)
has to
have zeros that would cancel those in the denominator. However, since the
measured data always contain some error, the exact cancelation is not likely
to happen, and one needs a sophisticated regularization scheme to keep the
total error bounded.
This problem can be avoided by replacing in (10) Bessel function J0 by
Hankel function H
ĝm,H(λ) =
gm(r)H
0 (λr)dr.
The addition theorem for H
0 (λ|x− y|) takes form
0 (λ|x− y|) =
Jm(λ|x|)H
m (λ|y|)e
−im(ϕ−θ),
and by proceeding as before one can obtain the following formula for fm(ρ):
fm(ρ) = Hm
ĝm,H(λ)
m (λ|R|)
m (λ|R|)
gm(r)H
0 (λr)dr
Unlike Jm, Hankel functions H
m (t) do not have zeros for all real values of t
and therefore problems with division by zeros do not arise in this amended
version of the method [82].
This derivation can be repeated in 3-D, with the exponentials eikθ replaced
by the spherical harmonics, and with cylindrical Bessel functions replaced by
their spherical counterparts. By doing this one will arrive at the Fourier
series method of [83]. Our use of Hankel function H
0 above is similar to
the way the authors of [83] utilized spherical Hankel function h
0 to avoid
the divisions by zero.
5.1.3 Filtered backprojection methods
The favorite way of inverting Radon transform for tomography purposes is
by using filtered backprojection type formulas, which involve filtration in
Fourier domain followed (or preceded) by a backprojection. In the case of
the set of spheres centered on a closed surface (e.g., sphere) S, one expects
such a formula to involve a filtration with respect to the radius variable and
then some integration over the set of spheres passing through the point of
interest. For quite a while, no such type formula had been discovered. This
did not prevent practitioners from reconstructions, since good approximate
inversion formulas (parametrices) could be developed, followed by an iterative
improvement of the reconstruction, see e.g. reconstruction procedures in
[114, 115, 118, 119, 120], and especially [96, 97].
The first set of exact inversion formulas of the filtered backprojection type
was discovered in [33]. These formulas were obtained only in odd dimensions.
Several different variations of such formulas (different in terms of exact order
of the filtration and backprojection steps) were developed. Let us denote by
g(p, r) = r2RSf the spherical integral, rather than the average, of f . Then
various versions of the 3D inversion formulas that reconstruct a function f(x)
supported inside S from its the spherical mean data RSf , read:
f(x) = − 1
g(y, |y − x|)dA(y),
f(x) = − 1
g(y, t)
) ∣∣∣∣∣
t=|y−x|
dA(y),
f(x) = − 1
g(y,t)
)) ∣∣∣∣∣
t=|y−x|
dA(y).
Recently, analogous formulas were obtained for even dimensions in [32]. De-
noting by g, as before the spherical integrals (rather than averages) of f , the
formulas of [32] in 2D look as follows:
f(x) =
g(y, t) log(t2 − |x− y|2) dt dl(y), (16)
f(x) =
g(y, t)
log(t2 − |x− y|2) dt dl(y), (17)
A different set of explicit inversion formulas that work in arbitrary dimensions
was presented in [69].
f(x) =
4(2π)n−1
n(y)h(y, |x− y|)dA(y). (18)
h(y, t) =
Y (λt)
J(λt′)g(y, t′)dt′
−J(λt)
Y (λt′)g(y, t′)dt′
λ2n−3dλ, (19)
J(t) =
Jn/2−1(t)
tn/2−1
, Y (t) =
Yn/2−1(t)
tn/2−1
Jn/2−1(t) and Yn/2−1(t) are respectively the Bessel and Neumann functions
of order n/2− 1, and n(y) is the vector of exterior normal to ∂B.
In 2-D equations (18), (19) can be simplified to yield the following recon-
struction formula
f(x) = −
g(y, t′)
|x− y|2 − t′2
 dl(y).
A similar simplification is also possible in 3D resulting in the formula
f(x) =
g(y, t)
) ∣∣∣∣∣
t=|y−x|
dA(y). (20)
Equation (20) is equivalent to one of the formulas derived in [116] for the 3D
case. It is interesting to notice that the “universal” formula of [116] holds for
all geometries when the backprojection type formulas are known: spherical,
cylindrical, and planar. It is not very likely that such explicit formulas would
be available for any closed surfaces S different from spheres (see a related
discussion in [15, 27]).
5.1.4 Series solutions for arbitrary geometries
Although, as we have just mentioned, we do not expect such explicit formulas
to be derived for non-spherical closed surfaces S, there is, however, a different
approach [70] that theoretically works for any closed S and that is practically
useful in some non-spherical geometries.
Let λ2m and um(x) be the eigenvalues and normalized eigenfunctions of
the Dirichlet Laplacian −∆ on the interior Ω of a closed surface S:
∆um(x) + λ
mum(x) = 0, x ∈ Ω, Ω ⊆ R
n, (21)
um(x) = 0, x ∈ S,
||um||
|um(x)|
2dx = 1.
As before, we would like to reconstruct a compactly supported function f(x)
from the known values of its spherical integrals g(y, r) with the centers on S:
g(y, r) =
f(y + rω)rn−1dω, y ∈ S.
We notice that um(x) is the solution of the Dirichlet problem for the Helmholtz
equation with zero boundary conditions and the wave number λm, and thus
it admits the Helmholtz representation
um(x) =
Φλm(|x− y|)
um(y)ds(y) x ∈ Ω, (22)
where Φλm(|x− y|) is a free-space rotationally invariant Green’s function of
the Helmholtz equation (21).
The eigenfunctions {um(x)}
0 form an orthonormal basis in L2(Ω). There-
fore, f(x) can be represented by the series
f(x) =
αmum(x) (23)
um(x)f(x)dx.
Since f(x) is C10 , series (23) converges pointwise. A reconstruction formula of
αm, and thus of f(x), will result if we substitute representation (22) into (23)
and interchange the order of integrations. Indeed, after a brief calculation
we will get
um(x)f(x)dx =
I(y, λm)
um(y)dA(x), (24)
where
I(y, λ) =
Φλ(|x− y|)f(x)dx. (25)
Certainly, the need to know the spectrum and eigenfunctions of the
Dirichlet Laplacian imposes a severe constraint on the surface S. However,
there are simple cases when the eigenfunctions are well known, and fast sum-
mation formulas for the corresponding series are available. Such is the case
of a cubic measuring surface S (see [70]); the eigenfunctions um are products
of sine functions
um(x) =
πm1x1
πm2x2
πm3x3
, (26)
where m = (m1, m2, m3), m1, m2, m3 ∈ N, and the eigenvalues are easily
found as well
λm = π
2|m|2/R2. (27)
Sum (23) is just a regular 3-D Fourier sine series easily computable by ap-
plication of the Fast Sine Fourier transform algorithm. The algorithmic
trick that allows one to calculate fast the coefficients αm consists in com-
puting first integrals (25) on a uniform mesh in λ. This is easily done by
a one-dimensional Fast Cosine Fourier transform algorithm, with Φλ(t) =
cos(λt)/t. The normal derivatives of um(x) are also products of sine func-
tions, this time two-dimensional ones. This, in turn, permits rapid evaluation
of integrals
I(y, λ) ∂
um(y)dA(x) for each mesh value of λ, and for each
one of the six faces ∂Ωi, i = 1, ..., 6 of the cube. Finally, the computation of
αm using equation (24) reduces to the interpolation in the spectral param-
eter λ, since the integrals in the right hand side of this equation have been
computed for the mesh values of this parameter (not for λm). Due to oscil-
latory nature of the integrals (25) a low order interpolation here would lead
to inaccurate reconstructions. Luckily, however, these integrals are analytic
functions of parameter λ (due to the finite support of g). Hence, high order
polynomial interpolation is applicable, and numerics yields very good results.
The algorithm we just described requires O(m3 logm) floating point op-
erations if the reconstruction is to be performed on an m×m×m Cartesian
grid, from comparably discretized data measured on a cubic surface. In prac-
tical terms, it yields reconstructions in the matter of several seconds on grids
with total number of nods exceeding a million [70].
5.1.5 Time reversal (backpropagation) methods
In the constant speed case, the following approach is possible in 3D: due
to the validity of the Huygens’ principle (i.e., the signal escapes from any
bounded domain in finite time), the pressure p(t, x) inside S will become
equal to zero for any time T larger than the time required to cross the domain
(i.e., time that it takes the sound to move along the diameter of S, which
for c = 1 equals the diameter). Thus, one can impose the zero conditions
on p(t, x) for t = T and solve the wave equation (2) back in time, using the
measured data g as the boundary values. The solution of this well posed
problem at t = 0 gives the desired source function f(x). Such methods have
been successfully implemented [22].
Although in 2D or in presence of sound speed variations, Huygens’ princi-
ple does not hold anymore, and thus the signal theoretically will stay forever,
one can find good approximate solutions using a similar approach [4, 18], see
discussions of the variable speed case below.
5.1.6 Examples of reconstructions and additional remarks about
the inversion formulas
• It is well known that different analytic inversion formulas in tomogra-
phy can behave differently in numerical implementation (e.g., in terms
of their stability), However, numerical implementation seems to show
that the analytic (backprojection type) formulas (15)-(20), in spite of
some of them being not equivalent, work equally well. See, for example
the results of an analytic formula reconstruction in 3D shown in Fig.
• It is worth noting that although formulas (15)-(16) and (18)-(20) will
yield identical results when applied to functions that can be represented
as the spherical mean Radon transform of a function supported inside
S, they are in general not equivalent when applied to functions with
larger supports. Simple examples (e.g., of f being the characteristic
Figure 5: A mathematical phantom in 3D (left) and its reconstruction using
an analytic inversion formula.
set of a large ball containing S) show that these two types of formulas
provide different reconstructions.
• An interesting observation is that backprojection formulas (15)-(20) do
not reconstruct the function f correctly inside the surface S, if f has
support reaching outside S. For instance, applying the reconstruction
formulas to the function RS(χ|x|≤3) leads to an incorrect reconstruction
of the value of f = χ|x|≤3 inside S = {|x| ≤ 1}. (Here by χV we denote
the characteristic function of the set V , i.e. it takes the value 1 in V
and zero outside. So, χ|x|≤3 is the characteristic function of the ball of
radius 3 centered at the origin.)
An another example: if one adds to the phantom shown in Fig. 5
two balls to the right of the surrounding sphere S, this leads to strong
artifacts, as seen on Fig. 6.
What is the reason for such a distortion? If one does not know in ad-
vance that f has support inside S, the backprojection formulas shown
before use insufficient information to recover a function with a larger
support, and thus uniqueness of reconstruction is lost. Then the for-
mulas misinterpret the data, wrongly assuming that they came form a
function supported inside S and thus reconstructing the function in-
correctly.
Notice that the series reconstruction of the preceding Section is free of
such problem. E.g., the reconstruction shown in Fig. 7 confirms this.
Figure 6: A perturbed reconstruction, due to presence of two additional balls
outside S (not shown on the picture).
Figure 7: In the phantom shown on the left, most disks are located outside
the square acquisition surface S indicated by the dotted line. This does not
perturb the reconstruction inside S (right).
5.2 Reconstruction in the variable speed case
We will assume here that the sound speed c(x) is smooth, positive, constant
for large x, and non-trapping. Although most analytic techniques we de-
scribed above do not work in the variable speed case, some formulas can be
derived and algorithms can be designed. This work is in a beginning stage
and the results described below most surely can and will be improved.
5.2.1 “Analytic” inversions
Let us denote by Ω the interior of the observation surface S, i.e. the area
where the object to be imaged is located. Consider in Ω the operator A =
−c2(x)∆ with zero Dirichlet conditions on the boundary S = ∂Ω. This
operator is self-adjoint, if considered in the weighted space L2(Ω; c−2(x)).
We also denote by E the operator of harmonic extension, which trans-
forms a function φ on S to a harmonic function on Ω which coincides with φ
on S.
The following result provides a formula for reconstructing f from the data
Theorem 7 [4] The function f(x) in (2) can be reconstructed in Ω as fol-
lows:
f(x) = (Eg|t=0)−
2 sin (τA
2 )E(gtt)(x, τ)dτ. (28)
The validity of this result hinges upon decay estimates for the solution (so
called local energy decay [29, 110, 111]), which hold under the non-trapping
condition. These estimates guarantee a qualified decay of the solution p(t, x)
inside any bounded region, e.g. in Ω, when time t increases. In odd dimen-
sions decay is exponential, but only polynomial in even dimensions. The
decay can be used instead of Huygens’ principle to solve the wave equation
backwards, starting at the infinite time. This leads to the formula (28).
Due to functions of the operator A being involved, it is not that clear how
explicit this formula can be made. For instance, it would be interesting to
see whether one can derive from (28) a backprojection inversion formula for
the case of a constant sound speed and S being a sphere (we have already
seen that such formulas are known).
5.2.2 Backpropagation
The exponential decay at large values of time can be used as follows: for a
sufficiently large T , one can assume that the solution is practically zero at
t = T . Thus, imposing zero initial conditions at t = T and solving in reverse
time direction, one arrives at t = 0 to an approximation of f(x) [18].
5.2.3 Eigenfunction expansions
One natural way to try to use the formula (28) is to use eigenfunction expan-
sion of the operator A in Ω (assuming that such expansion is known). This
immediately leads to the following result:
Theorem 8 Under the same conditions on the sound speed as before, func-
tion f(x) can be reconstructed inside Ω from the data g in (2), as the following
L2(B)-convergent series:
f(x) =
fkψk(x), (29)
where the Fourier coefficients fk can be recovered using one of the following
formulas:
fk = λ
k gk(0)− λ
sin (λkt)g
k(t)dt,
fk = λ
k gk(0) + λ
cos (λkt)g
k(t)dt, or
fk = −λ
sin (λkt)gk(t)dt = −λ
sin (λkt)g(x, t)
(x)dxdt,
gk(t) =
g(x, t)
(x)dx.
Here ν denotes the external normal to S.
One notices that this is a generalization to the variable sound speed case
of the expansion method of [70], discussed in Section 5.1.4. An interesting
feature is that, unlike in [70], we do not need to know the whole space Green’s
function for A (which is certainly not known).
It is not clear yet how feasible numerical implementation of this approach
could be.
6 Partial data. “Visible” and “invisible” sin-
gularities
Uniqueness of reconstruction does not imply practical recoverability, since the
reconstruction procedure might be severely unstable. This is well known to
be the case, for instance, in incomplete data situations in X-ray tomography,
and even for complete data problems in some imaging modalities, such as
the electrical impedance tomography [64, 68, 76, 77].
In order to describe the results below, we need to explain the notion
of the wave front set WF (f) of a function f(x). This set carries detailed
information on singularities of f(x). It consists of pairs (x, ξ) of a point x
in space and a wave vector (Fourier domain variable) ξ 6= 0. It is easier to
say what it means that a point (x0, ξ0) is not in the wave front set WF (f).
This means that one can smoothly cut-off f to zero at a small distance from
x0 in such a way that the Fourier transform φ̂f(ξ) of the resulting function
φ(x)f(x) decays faster than any power of ξ in directions that are close to
the direction of ξ0. We remind the reader that if this Fourier transform
decays that way in all directions, then f(x) is smooth near the point x0. So,
the wave front set contains pairs (x0, ξ0) such that f is not smooth near x0,
and ξo indicates why it is not: the Fourier transform does not decay well in
this direction. For instance, if f(x) consists of two smooth pieces joined non-
smoothly across a smooth interface Σ, then WF (f) contains pairs (x, ξ) such
that x is in Σ and ξ is normal to Σ at x. One can find simple introduction
to the notions of microlocal analysis, such as the wave front set, for instance
in [106].
Analysis done in [99] for the constant speed case (equivalently, for the
spherical mean transform RS), showed which parts of the wave front (and
thus singularities) of a function f can be recovered from its partial X-ray
data. An analog of this result also holds for the spherical mean transform
RS [73] (see also [120] for a practical discussion). We formulate it below in
an imprecise form (see [73] for precise formulation).
Theorem 9 [73] A wavefront set point (x, ξ) of f is “stably recoverable”
from RSf if and only if there is a circle (sphere in higher dimensions) centered
on S, passing through x, and normal to ξ at this point.
As we have already mentioned, this result does not exactly hold the way it is
formulated and needs to include some precise conditions (see [73, Theorem
3]). The statement is, for instance, correct if S is a smooth hypersurface and
the support of f lies on one side of the tangent plane to S at the center of
the sphere mentioned in the theorem.
Talking about jump singularities only (i.e., interfaces between smooth
regions inside the object to be imaged), this result says that in order for a
piece of the interface to be stably recoverable (dubbed “visible”), one should
have for each point of this interface, a sphere centered at S and tangent to
the interface at this point. Otherwise, the interface will be blurred away
(even if there is a uniqueness of reconstruction theorem). The reason is that
if all spheres of integration are transversal to the interface, the integration
smoothes off the singularity, and therefore its recovery becomes highly unsta-
ble (numerically, one has to deal with inversion of a matrix with exponentially
fast decaying singular values). The Figure 8 below shows an example of an
incomplete data reconstruction from spherical mean data. One sees clearly
the effect of disappearance of the parts of the boundaries that are not touched
tangentially by circles centered at transducers’ locations.
Figure 8: Effect of incomplete data: the phantom (left) and its incomplete
data reconstruction. The transducers were located along a 180o circular arc
(the left half of a large circle surrounding the squares).
7 Range conditions
As it has already been mentioned, the space of functions g(t, y) that could
arise as exact data measured by transducers (i.e., the range of the data),
is very small (of infinite codimension in the spaces of all functions of t >
0, y ∈ S). Knowing this space (range) is useful for many theoretical and
practical purposes (reconstruction algorithms, error corrections, incomplete
data completion, etc.), and thus has attracted a lot of attention (e.g., [30,
38, 39, 40, 53, 54, 64, 66, 67, 68, 74, 76, 77, 78, 90, 100].
For instance, for the standard Radon transform
f(x) → g(s, ω) =
x·ω=s
f(x)dx, |ω| = 1,
the range conditions on g(s, ω) are:
1. evenness: g(−s,−ω) = g(s, ω)
2. moment conditions: for any integer k ≥ 0, the kth moment
Gk(ω) =
skg(ω, s)ds
extends from the unit circle of vectors ω to a homogeneous polynomial
of degree k in ω.
The evenness condition is obviously necessary and is kind of “trivial”. It
seems that the only non-trivial conditions are the moment ones. However,
here the standard Radon transform misleads us, as it often happens. In fact,
for more general transforms of Radon type it is often easy (or easier) to find
analogs of the moment conditions, while analogs of the evenness conditions
are often elusive (see [64, 66, 67, 76, 77, 84] devoted to the case of SPECT
(single photon emission tomography)). The same happens in TAT.
Let us deal first with the case of a constant sound speed, when one can
think of the spherical mean transform RS instead of the wave equation model.
An analog of the moment conditions was already present implicitly (without
saying that these were range conditions) in [71, 72, 7] and explicitly formu-
lated as such in [95]. Indeed, our discussion in Section 4 of the polynomials
Qk provides the following conditions of the moment type:
Moment conditions [7, 71, 72, 95] on data g(p, r) = RSf(p, r) look as
follows: for any integer k ≥ 0, the moment
Mk(ω) =
r2k+d−1g(p, r)dr
can be extended from S to a (non-homogeneous) polynomial Qk(x) of degree
at most 2k.
These conditions, however, are incomplete, and in fact infinitely many
others, which play the role of an analog of evenness, need to be added.
Complete range descriptions for RS when S is a circle in 2D were discov-
ered in [13] and then in odd dimensions in [34]. They were then extended
to any dimension and interpreted in several different ways in [6]. These
conditions happen to be intimately related to PDEs and spectral theory.
Figure 9:
In order to describe these conditions, we need to introduce some notations.
Let B be the unit ball in Rd, S - the unit sphere, and C - the cylinder B×[0, 2]
(see Fig. 9).
We introduce the spherical mean operator RS as before:
RSf(x, t) =
|y|=1
f(x+ ty)dA(y), x ∈ S.
Several different range descriptions for RS were provided in [6], out of
which we only show a few:
Theorem 10 [6] The following three statements are equivalent:
1. The function g ∈ C∞0 (S × [0, 2]) is representable as RSf for some
f ∈ C∞0 (B). (In other words, g represents an ideal (free of errors) set
of TAT data.)
2. (a) The moment conditions are satisfied.
(b) Let −λ2 be any eigenvalue of the Laplace operator in B with zero
Dirichlet conditions and ψλ be the corresponding eigenfunction.
Then the following orthogonality condition is satisfied:
S×[0,2]
g(x, t)∂νψλ(x)jn/2−1(λt)t
n−1dxdt = 0. (31)
Here jp(z) = cp
Jp(z)
is the so called spherical Bessel function.
3. (a) The moment conditions are satisfied.
(b) Let ĝ(x, λ) =
g(x, t)jn/2−1(λt)t
n−1dt. Then, for any m ∈ Z, the
mth spherical harmonic term ĝm(x, λ) of ĝ(x, λ) vanishes at all
zeros λ 6= 0 of Bessel function Jm+n/2−1(λ).
Remark 11 [6]
1. In odd dimensions, moment conditions are not necessary, and thus con-
ditions 2(b) or 3(b) suffice. (A similar earlier result was established for
a related transform in [34].)
2. The range conditions (2) of the previous Theorem are also necessary
when S is the boundary of any bounded domain, not necessarily a
sphere.
3. An analog of these conditions can be derived for a variable sound speed
(without non-trapping conditions imposed).
8 Concluding remarks
8.1 Variations of the TAT procedure
8.1.1 Planar and linear transducers
Assuming that transducers are point-like, is clearly an approximation, and
in fact a transducer measures the average pressure over its area. It has
been rightfully claimed that the point approximation for transducers should
lead to some blurring in the reconstructions. This, as well as intricacies of
reconstructions from the data obtained by point transducers, triggered recent
proposals for different types of transducers (see [20, 21], [47]-[52], [92, 93]).
In these papers, it was suggested to use either planar, or line detectors.
In the first case [47], the detectors are assumed to be large and planar,
ideally assumed to be approximations of infinite planes that are placed tan-
gentially to a sphere containing the object. Thus, the data one collects is
the integrals of the pressure over these planes, for all values of t > 0. If one
takes the standard 3D Radon transform of the pressure p(x, t) with respect
to x:
P (x, t) 7→ q(s, t, ω) =
x·ω=s
p(x, t)dA(x),
where dA is the surface measure and ω is a unit vector in R3, this is well
known to reduce the 3D Laplace operator ∆x to the second derivative ∂
2/∂s2
[30, 38, 39, 40, 53, 54], and thus the 3D wave equation to the string vibra-
tion problem. The measured data provide the boundary conditions for this
problem. The initial conditions in (1) mean evenness with respect to time,
and thus the standard d’Alambert formula leads to the immediate realiza-
tion that the measured data is just the 3D Radon transform of f(x). Thus,
the reconstruction boils down to the well known inversion formulas for the
Radon transform.
Another proposal ([20, 21], [49]-[52], [92, 93]) is to use line detectors
that provide line integrals of the pressure p(x, t). Such detectors can be
implemented optically, using either Fabry-Perot [20], or Mach-Zehnder [93]
interferometers.
Suppose that the object is surrounded by a surface that is rotation in-
variant with respect to the z-axis. It is suggested to place the line detectors
perpendicular to the z-axis and tangential to the surface. The same consid-
eration as above then shows that after the 2D Radon (or X-ray, which in
2D is the same) transform in each plane orthogonal to z-axis, the 3D wave
equation converts into the 2D one for the Radon data. The measurements
provide the boundary data. Thus, the reconstruction boils down to solving
a 2D problem similar to the one in the case of point detectors, and then
inverting the 2D Radon transform.
Due to the recent nature of these two projects, it appears to be too early
to judge which one will be superior in the end. For instance, it is not clear
beforehand, whether the approximation of infinite size (length, area) of the
linear or planar detectors works better than the zero dimension approxima-
tion for point detectors. Further developments will resolve these questions.
8.1.2 Direct imaging techniques
Some direct imaging techniques have been suggested, which might not require
mathematical reconstructions. See, for instance, [79] about an acoustic lens
system.
8.1.3 Using contrast agents
Contrast agents to improve TAT imaging have been developed (e.g., [24]).
8.1.4 Passive thermoacoustic imaging
The TAT model we have considered can be called “active thermoacoustic to-
mography,” due to the set-up when the practitioner creates the signal. There
has been some recent development of the “passive thermoacoustic tomogra-
phy,” where the thermoacoustic signal is used to image the temperature
sources present inside the body. One can find a survey of this area in [94].
8.2 Uniqueness
8.2.1 Sketch of the proof of Theorem 10
We provide here a brief outline of the rather technical proof of Theorem 10.
Suppose that f is compactly supported, not identically zero, and such
that RSf = 0. Our previous considerations show that one can assume that
S is an algebraic curve (not a straight line) that is contained in the set of
zeros of a non-trivial harmonic polynomial. Now one touches the boundary
of the support of f from outside by a circle centered on S. Then microlocal
analysis of the operator RS (which happens to be an analytic Fourier Integral
Operator, FIO [19, 42, 43, 44, 45, 65, 98]) shows that, due to the equality
RSf = 0, at the tangency point the vector co-normal to the sphere should not
belong to the analytic wave front of f (microlocal regularity of solutions of
RSf = 0). This, for instance, can be also extracted from the results of [105].
On the other hand, a theorem by Hörmander and Kashiwara [56, Theorem
8.5.6] shows that this vector must be in the analytic wave front set, since
f = 0 on one side of the sphere (a microlocal version of uniqueness of analytic
continuation). This way, one gets a contradiction. Unfortunately, the life is
not so easy, and the proof sketched above does not go through smoothly, due
to possible cancelation of wavefronts at different tangency points. Then one
has to involve the geometry of zeros of harmonic polynomials [37] to exclude
the possibility of such a cancelation.
Thus, the proof uses microlocal analysis and geometry of zeros of har-
monic polynomials. Both these tools have their limitations. For instance,
the microlocal approach (at least, in the form it is used in [7]) does not al-
low considerations of non-compactly supported functions. Thus, the validity
of the Theorem for arbitrarily fast decaying, but not compactly supported,
functions is still not established, albeit it most certainly holds. On the other
hand, the geometric part does not work that well in dimensions higher than
two. Development of new approaches is apparently needed in order to over-
come these hurdles. A much simpler PDE approach has emerged recently
[33] (see also [12] and the next Section), albeit its achievements have been
limited so far.
8.2.2 Some open problems concerning uniqueness
As it has already been mentioned, one can consider the practical problems
about uniqueness resolved. However, the mathematical understanding of the
uniqueness problem for the restricted spherical mean operators RS is still
unsatisfactory. Here are some questions that still await their resolution:
1. Describe uniqueness sets in dimensions larger than 2 (prove the Con-
jecture 4). Recent limited progress, as well as variations on this theme
can be found in [1]-[12].
2. Prove Theorem 3 without using microlocal and harmonic polynomial
tools.
3. Prove Theorem 3 on uniqueness sets S under the condition of suffi-
ciently fast decay (rather than compactness of support) of the function.
Very little is known for the case of functions without compact support.
The main known result is of [3], which describes for which values of
1 ≤ p ≤ ∞ the result of Corollary 2 still holds:
Theorem 12 [3] Let S be the boundary of a bounded domain in Rd and
f ∈ Lp(Rd) such that RSf ≡ 0. If p ≤ 2d/(d−1), then f ≡ 0 (and thus
S is injectivity set for this space). This fails for any p > 2d/(d− 1).
8.3 Inversion
Albeit closed form (backprojection type) inversion formulas are available now
for the cases of S being a plane (and object on one side from it), cylinder,
and a sphere, there is still some mystery surrounding this issue.
1. Can one write a backprojection type inversion formula in the case of
the constant sound speed for a closed surface S which is not a sphere?
We suspect that the answer to this question is negative (see also related
discussion in [15, 27]).
2. The inversion formulas for S being a sphere assume that the object to
be imaged is inside S. One can check on simplest examples that if the
support of function f(x) reaches outside S, the inversion formulas do
not reconstruct the function correctly even inside of S. See [5] for a
discussion.
3. The I. Gelfand’s school of integral geometry has developed a marvelous
machinery of the so called κ operator, which provides a general ap-
proach to inversion and range descriptions for transforms of Radon type
[38, 39]. In particular, it has been applied to the case of integration of
various collections (“complexes”) of spheres in [39, 41]. This consider-
ation seems to suggest that one should not expect explicit closed form
inversion formulas for RS when S is a sphere. We, however, know that
such formulas have been discovered recently [33, 69]. This apparent
controversy has not been resolved.
4. Can one derive any more explicit analytic formulas from (28)?
5. Can the series expansion formulas of Theorem 8 be efficiently imple-
mented?
One can also mention that in some works [15, 23] it is suggested to use in
the TAT problem not only the values of the pressure measured by transducers
on the observation surface S, but its normal derivative to S as well. If one
knows both, then taking Fourier transform in the time variable and using the
whole space Green’s function for the Helmholtz equation leads immediately
to a reconstruction formula for the solution (which seems to be much simpler
than what is proposed in [23]). The problem is that this normal derivative is
not measured by TAT devices. Under some circumstances (e.g., when there
are no sources of ultrasound outside S), one can prove the theoretical pos-
sibility of recovering the missing normal derivative. This, however, does not
seem to us to be a plausible procedure. In rare cases (planar, cylindrical,
or spherical surface S), when involvement of the normal derivative can be
eliminated (e.g., [15, 27]), this might lead to feasible inversion algorithms,
but in these cases, as explained before in this text, explicit and nicely im-
plementable analytic inversion formulas are available. So, jury is still out on
this issue as well.
8.4 Stability
Stability of inversion when S is a sphere surrounding the support of f(x) is
the same as for the standard Radon transform, as the results of [91] and sec-
ond statement of Theorem 11 show. However, if the support reaches outside,
albeit Corollary 2 still guarantees uniqueness of reconstruction, stability (at
least for the parts outside S) is gone. Indeed, Theorem 9 shows that some
parts of singularities of f outside S will not be stably “visible.”
8.5 Range
As Theorem 9 states, the range conditions 2 and 3 of Theorem 10 are neces-
sary also for non-spherical closed surfaces S and for functions with support
outside S. They, however, are not expected to be sufficient, since Theorem 9
indicates that one might expect non-closed ranges in some cases. The same
applies for non-constant sound speed case.
Acknowledgments
The work of the first author was partially supported by the NSF DMS grants
0604778 and 0648786. The second author was partially supported by the
DOE grant DE-FG02-03ER25577 and NSF DMS grant 0312292. Part of
this work was completed when the first author was at the Isaac Newton
Institute for Mathematical Sciences. The authors express their gratitude to
the NSF, DOE and INI for this support. The authors thank M. Agranovsky
for extremely useful discussions, M. Anastasio, G. Beylkin and M. Klibanov
for providing preprints and references, and the reviewers for very helpful
remarks.
References
[1] Agranovsky, M. 1997 Radon transform on polynomial level sets and
related problems. Israel Math. Conf. Proc., 11, 1-21.
[2] Agranovsky, M. 2000 On a problem of injectivity for the Radon trans-
form on a paraboloid. Analysis, geometry, number theory: the mathe-
matics of Leon Ehrenpreis (Philadelphia, PA, 1998). Contemp. Math.
251, AMS, Provodence, RI, 1–14.
[3] Agranovsky, M., Berenstein, C., & Kuchment, P. 1996 Approximation
by spherical waves in Lp-spaces, J. Geom. Anal. 6, no. 3, 365–383.
[4] Agranovsky, M. & Kuchment, P. 2007 Uniqueness of reconstruction and
an inversion procedure for thermoacoustic and photoacoustic tomogra-
phy with variable sound speed. Inverse Problems 23, 2089–2102.
[5] Agranovsky, M., Kuchment, P. & Kunyansky, L. 2007 On reconstruc-
tion formulas and algorithms for the thermoacoustic and photoacoustic
tomography. Submitted.
[6] Agranovsky, M., Kuchment, P., & Quinto, E. T. 2007 Range descriptions
for the spherical mean Radon transform. J. Funct. Anal. 248, 344–386.
[7] Agranovsky, M. & Quinto, E. T. 1996 Injectivity sets for the Radon
transform over circles and complete systems of radial functions. Journal
of Functional Analysis, 139, 383–414.
[8] Agranovsky, M. & Quinto, E. T. 2001 Geometry of stationary sets for
the wave equation in Rn:the case of finitely suported initial data. Duke
Math. J., 107, no. 1, 57–84.
[9] Agranovsky, M. & Quinto, E. T. 2003 Stationary sets for the wave equa-
tion in crystallographic domains. Trans. AMS, 355, no. 6, 2439–2451.
[10] Agranovsky, M. & Quinto, E. T. 2006 Remarks on stationary sets for the
wave equation. Integral Geometry and Tomography, Contemp. Math.
405, 1–11.
[11] Agranovsky, M., Volchkov, V. V. & Zalcman, L. 1999 Conical uniqueness
sets for the spherical Radon transform. Bull. London Math. Soc., 31, no.
4, 363–372.
[12] Ambartsoumian, G. & Kuchment, P. 2005 On the injectivity of the
circular Radon transform. Inverse Problems 21, 473–485.
[13] Ambartsoumian, G. & Kuchment, P. 2006 A range description for the
planar circular Radon transform. SIAM J. Math. Anal. 38, no. 2, 681–
[14] Ambartsoumian, G. & Patch, S. 2007 Thermoacoustic tomography: nu-
merical results. Proceedings of SPIE 6437, Photons Plus Ultrasound:
Imaging and Sensing 2007: The Eighth Conference on Biomedical Ther-
moacoustics, Optoacoustics, and Acousto-optics, Alexander A. Oraevsky,
Lihong V. Wang, Editors, 64371B.
[15] Anastasio, M. A., Zhang, J., Modgil, D., and La Rivière, P. J. 2007
Application of inverse source concepts to photoacoustic tomography,
preprint.
[16] Andersson, L.-E.. 1988 On the determination of a function from spherical
averages. SIAM J. Math. Anal. 19 no. 1, 214–232.
[17] Asgeirsson, L. 1937 Über eine Mittelwerteigenschaft von Lösungen ho-
mogener linearer partieller Differentialgleichungen zweiter Ordnung mit
konstanten Koeffizienten. Ann. Math., 113, 321–346.
[18] Bangerth, W., Georgieva-Hristova, Y. & Kuchment, P. 2007 On recon-
struction in thermoacoustic tomography with variable speed, in prepa-
ration.
[19] Beylkin, G. 1984 The inversion problem and applications of the gener-
alized Radon transform. Comm. Pure Appl. Math. 37, 579–599.
[20] Burgholzer, P., Hofer, C., Paltauf, G., Haltmeier, M., & Scherzer, O.
2005 Thermoacoustic tomography with integrating area and line detec-
tors. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency
Control 52(9), 1577–1583.
[21] Burgholzer, P.,, Hofer, C., Matt, G. J., Paltauf, G., Haltmeier, M.
& Scherzer, O. 2006 Thermoacoustic tomography using a fiber-based
Fabry-Perot interferometer as an integrating line detector. Proc. SPIE
6086, 434–442.
[22] Burgholzer, P., Matt, G., Haltmeier, M. & Patlauf, G. 2007 Exact and
approximate imaging methods for photoacoustic tomography using an
arbitrary detection surface. Phys. Rev. E 75, 046706.
[23] Clason, C. and Klibanov, M. 2007 Quasireversibility method in ther-
moacoustic tomography in heterogeneous medium, preprint.
[24] Copland, J. A. et al. 2004 Bioconjugated gold nanoparticles as a molec-
ular based contrast agent: implications for imaging of deep tumors using
optoacoustic tomography. Molecular Imaging and Biology 6, no. 5, 341–
[25] Courant, R. & Hilbert, D. 1962 Methods of Mathematical Physics, Vol-
ume II Partial Differential Equations, Interscience, New York.
[26] Denisjuk, A. 1999 Integral geometry on the family of semi-spheres. Fract.
Calc. Appl. Anal. 2, no. 1, 31–46.
[27] Devaney, A. J. and Beylkin, G. 1984 Diffraction tomography using arbi-
trary transmitter and receiver surfaces. Ultrasonic Imaging 6, 181-193.
1984.
[28] Diebold, G. J., Sun, T. & Khan, M. I. 1991 Photoacoustic monopole
radiation in one, two, and three dimensions. Phys. Rev. Lett. 67, no.
24, 3384–3387.
[29] Egorov, Yu. V. & Shubin, M. A. 1992 Partial Differential Equations I.
Encyclopaedia of Mathematical Sciences, (Springer Verlag), 30, 1–259
[30] Ehrenpreis, L. 2003 The Universality of the Radon Transform, Oxford
Univ. Press.
[31] Fawcett, J. A. 1985 Inversion of n-dimensional spherical averages. SIAM
J. Appl. Math. 45, no. 2, 336–341.
[32] Finch, D., Haltmeier, M. & Rakesh 2007 Inversion of spherical
means and the wave equation in even dimensions. Preprint arXiv
math.AP/0701426.
[33] Finch, D., Patch, S. & Rakesh 2004 Determining a function from its
mean values over a family of spheres. SIAM J. Math. Anal. 35, no. 5,
1213–1240.
[34] Finch, D. & Rakesh 2006 The range of the spherical mean value operator
for functions supported in a ball. Inverse Problems 22, 923-938.
[35] Finch, D. & Rakesh 2007 Recovering a function from its spherical mean
values in two and three dimensions. Preprint.
http://arxiv.org/abs/math/0701426
[36] Finch, D. & Rakesh 2007 The spherical mean value operator with centers
on a sphere. Preprint. To appear in Inverse Problems.
[37] Flatto, L., Newman, D. J. & Shapiro, H. S. 1966 The level curves of
harmonic functions. Trans. Amer. Math. Soc. 123, 425–436.
[38] Gelfand, I., Gindikin, S. & Graev M. 1980 Integral geometry in affine
and projective spaces. J. Sov. Math. 18, 39–167.
[39] Gelfand, I., Gindikin, S. & Graev M. 2003 Selected Topics in Integral
Geometry. Transl. Math. Monogr. v. 220, Amer. Math. Soc., Providence
[40] Gelfand, I., Graev M. & Vilenknin, N. 1965 Generalized Functions, v.
5: Integral Geometry and Representation Theory, Acad. Press.
[41] Gindikin, S. 1995 Integral geometry on real quadrics, in Lie groups and
Lie algebras: E. B. Dynkin’s Seminar, 23–31, Amer. Math. Soc. Transl.
Ser. 2, 169, Amer. Math. Soc., Providence, RI.
[42] Greenleaf, A. & Uhlmann, G. 1990 Microlocal techniques in integral
geometry. Contemporary Math. 113, 149–155.
[43] Guillemin, V. 1975 Fourier integral operators from the Radon transform
point of view. Proc. Symposia in Pure Math., 27, 297–300.
[44] Guillemin, V. 1985 On some results of Gelfand in integral geometry.
Proc. Symposia in Pure Math., 43, 149–155.
[45] Guillemin, V. & Sternberg S. 1977 Geometric Asymptotics. Amer. Math.
Soc., Providence, RI.
[46] Gusev, V. E. & Karabutov, A. A. 1993 Laser Optoacoustics. American
Inst. of Physics, NY.
[47] Haltmeier, M., Burgholzer, P., Paltauf, G. & Scherzer, O. 2004 Ther-
moacoustic computed tomography with large planar receivers. Inverse
Problems 20, 1663–1673.
[48] Haltmeier, M., Schuster, T. & O. Scherzer. 2005 Filtered backprojec-
tion for thermoacoustic computed tomography in spherical geometry.
Mathematical Methods in the Applied Sciences, 28, 1919–1937.
[49] Haltmeier, M., Paltauf, G., Burgholzer, P. & Scherzer, O. 2005 Ther-
moacoustic Tomography with integrating line detectors. Proc. SPIE
5864:586402-8.
[50] Haltmeier, M., Burgholzer, P., Hofer, C., Paltauf, G., Nuster, R. &
Scherzer, O. 2005 Thermoacoustic tomography using integrating line
detectors. Ultrasonics Symposium 1, 166–169.
[51] Haltmeier, M., Scherzer, O., Burgholzer, P. & Paltauf, G. 2005 Ther-
moacoustic Computed Tomography with large planar receivers. ECMI
Newsletter 37, pp. 31-34. http://www.it.lut.fi/mat/EcmiNL/ecmi37/
[52] Haltmeier, M., & Fidler, T. Mathematical Challenges Aris-
ing in Thermoacoustic Tomography with Line Detectors, preprint
arXiv:math.AP/0610155.
[53] Helgason, S. 1980 The Radon Transform, Birkhäuser, Basel.
[54] Helgason, S. 2000 Groups and Geometric Analysis. Amer. Math. Soc.,
Providence, R.I.
[55] Herman, G.(Ed.) 1979 Image Reconstruction from Projections . Topics
in Applied Physics, v. 32, Springer Verlag, Berlin, New York.
[56] Hörmander, L. 1983 The Analysis of Linear Partial Differential Opera-
tors, vol. 1, Springer-Verlag, New York.
[57] Jin, X. & Wang, L. V. 2006 Thermoacoustic tomography with correction
for acoustic speed variations. Physics in Medicine and Biology 51, 6437–
6448.
[58] John, F. 1971 Plane Waves and Spherical Means Applied to Partial
Differential Equations, Dover.
[59] Kak, A. C. & Slaney, M. 2001 Principles of Computerized Tomographic
Imaging. SIAM, Philadelphia.
[60] Köstli K. P., Frenz, M., Bebie, H. & Weber H. P. 2001 Temporal back-
ward projection of optoacoustic pressure transients using Fourier trans-
form methods. Phys. Med. Biol. 46, 1863–1872
http://www.it.lut.fi/mat/EcmiNL/ecmi37/
http://arxiv.org/abs/math/0610155
[61] Kruger, R. A., Kiser, W. L., Reinecke, D. R. & Kruger, G. A. 2003 Ther-
moacoustic computed tomography using a conventional linear trans-
ducer array. Med. Phys. 30, no 5, 856–860.
[62] Kruger, R. A., Liu. P., Fang, Y. R. & Appledorn, C. R. 1995 Photoa-
coustic ultrasound (PAUS)reconstruction tomography. Med. Phys. 22,
1605–1609.
[63] Kuchment, P. 1993, unpublished.
[64] Kuchment, P. 2006 Generalized Transforms of Radon Type and Their
Applications. in [85], pp. 67–91.
[65] Kuchment, P., Lancaster, K. & Mogilevskaya, L. 1995 On local tomog-
raphy. Inverse Problems, 11, 571–589.
[66] Kuchment, P. & Lvin S. 1990 Paley-Wiener theorem for the exponential
Radon transform. Acta Applicandae Mathematicae, no.18, 251–260.
[67] Kuchment, P. & Lvin S. 1991 The Range of the Exponential Radon
Transform. Soviet Math Dokl, 42, no.1, 183–184.
[68] Kuchment, P. & Quinto, E. T. 2003. Some problems of integral geometry
arising in tomography. Chapter XI in [30].
[69] Kunyansky, L. 2007 Explicit inversion formulae for the spherical mean
Radon transform. Inverse problems 23 (2007), 737-783.
[70] Kunyansky, L. 2007 A series solution and a fast algorithm for the
inversion of the spherical mean Radon transform. Preprint arXiv
math.AP/0701236.
[71] Lin, V. & Pinkus, A. 1993 Fundamentality of ridge functions. J. Approx.
Theory, 75, 295–311.
[72] Lin, V. & Pinkus, A 1994 Approximation of multivariate functions. In
Advances in Computational Mathematics, H. P. Dikshit & C. A. Mic-
chelli, Eds., World Sci. Publ., 1–9.
[73] Louis, A. K. & Quinto, E. T. 2000 Local tomographic methods in Sonar.
In Surveys on solution methods for inverse problems, Springer, Vienna,
147–154.
http://arxiv.org/abs/math/0701236
[74] Lvin S. 1994 Data correction and restoration in emission tomography,
in E.T. Quinto, M. Cheney, and P. Kuchment (Editors), Tomography,
Impedance Imaging, and Integral Geometry, Lectures in Appl. Math.,
vol. 30, AMS, Providence, RI, 149–155.
[75] Mathematics and Physics of Emerging Biomedical Imaging,
The National Academies Press 1996. Available online at
http://www.nap.edu/catalog.php?record id=5066#toc.
[76] Natterer, F. 1986 The mathematics of computerized tomography, Wiley,
New York.
[77] Natterer, F. & Wübbeling, F. 2001 Mathematical Methods in Image Re-
construction, Monographs on Mathematical Modeling and Computation
5, SIAM, Philadelphia, PA.
[78] Nessibi, M. M., Rachdi, L. T. & Trimeche, K. 1995 Ranges and inversion
formulas for spherical mean operator and its dual. J. Math. Anal. Appl.
196, no. 3, 861–884.
[79] Niederhauser, J. J., Jaeger, M., Lemor, R., Weber, P. & Frenz, M.
2005 Combined ultrasound and optoacoustic system for real-time high-
contrast vascular imaging in vivo. IEEE Transactions on medical Imag-
ing 24, 436–440.
[80] Nilsson, S. 1997 Application of fast backprojection techniques for some
inverse problems of integral geometry. Linkoeping studies in science and
technology, Dissertation 499, Dept. of Mathematics, Linkoeping univer-
sity, Linkoeping, Sweden.
[81] Nolan, C. J. & Cheney, M. 2002 Synthetic aperture inversion. Inverse
Problems 18, 221–235.
[82] Norton, S. J. 1980 Reconstruction of a two-dimensional reflecting
medium over a circular domain: exact solution. J. Acoust. Soc. Am.
67, 1266–1273.
[83] Norton, S. J. & Linzer, M. 1981 Ultrasonic reflectivity imaging in three
dimensions: exact inverse scattering solutions for plane, cylindrical, and
spherical apertures. IEEE Transactions on Biomedical Engineering, 28,
200–202.
http://www.nap.edu/catalog.php?record_id=5066#toc
[84] Novikov, R. 2002 On the range characterization for the two-dimensional
attenuated X-ray transform. Inverse Problems 18, 677–700.
[85] Olafsson, G. & Quinto, E. T. (Editors), 2006 The Radon Transform, In-
verse Problems, and Tomography. American Mathematical Society Short
Course January 3–4, 2005, Atlanta, Georgia, Proc. Symp. Appl. Math.,
v. 63, AMS, RI.
[86] Oraevsky, A. A., Esenaliev, R. O., Jacques, S. L. Tittel, F. K. 1996
Laser optoacoustic tomography for medical diagnostics principles. Proc.
SPIE 2676, 22.
[87] Oraevsky, A. A. & Karabutov, A. A. 2002 In Handbook of Optical
Biomedical Diagonstics, edited by V. V. Tuchin, SPIE, Bellingham, WA,
Chap. 10.
[88] Oraevsky A. A.& A. A. Karabutov, A. A., 2003 Optoacoustic Tomogra-
phy, Ch. 34 In Biomedical Photonics Handbook, edited by T. Vo-Dinh,
CRC, Boca Raton, FL, Chap. 34, 34-1 – 34-34.
[89] Palamodov, V. P. 2000 Reconstruction from limited data of arc means.
J. Fourier Anal. Appl. 6, no. 1, 25–42.
[90] Palamodov, V. P. 2004 Reconstructive Integral Geometry, Birkhäuser,
Basel.
[91] Palamodov, V. 2006 Remarks on the general Funk transform. Preprint,
Tel Aviv University, August.
[92] Paltauf, G., Burgholzer, P., Haltmeier, M. & O. Scherzer 2005 Ther-
moacoustic Tomography using optical Line detection. Proc. SPIE 5864,
7–14.
[93] Paltauf, G., R. Nuster, Haltmeier, M. & Burgholzer, P. 2007(?) Ther-
moacoustic Computed Tomography using a Mach-Zehnder interferome-
ter as acoustic line detector. Submitted
[94] Passechnik, V. I., Anosov, A. A. & Bograchev, K. M. 2000 Fundamentals
and prospects of passive thermoacoustic tomography. Critical reviews in
Biomed. Eng. 28, no. 3&4, 603–640.
[95] Patch, S. K. 2004 Thermoacoustic tomography - consistency conditions
and the partial scan problem. Phys. Med. Biol. 49, 1–11.
[96] Popov D. A. & Sushko, D. V. 2002 A parametrix for the problem of
optical-acoustic tomography. Dokl. Math. 65, no. 1, 19–21.
[97] Popov D. A. & Sushko, D. V. 2004 Image restoration in optical-acoustic
tomography. Problems of Information Transmission 40, no. 3, 254–278.
[98] Quinto, E. T. 1980 The dependence of the generalized Radon transform
on defining measures. Trans. Amer. Math. Soc. 257, 331–346.
[99] Quinto, E. T. 1993 Singularities of the X-ray transform and limited data
tomography in R2 and R3. SIAM J. Math. Anal. 24, 1215–1225.
[100] Quinto, E. T. 2006 An introduction to X-ray tomography and Radon
transforms. In [85], 1–23.
[101] Ramm, A. G. 1985 Inversion of the backscattering data and a problem
of integral geometry. Phys. Lett. A 113, no. 4, 172–176.
[102] Ramm, A. G. 2002 Injectivity of the spherical means operator. C. R.
Math. Acad. Sci. Paris 335, no. 12, 1033–1038.
[103] Romanov, V. G. 1967 Reconstructing functions from integrals over a
family of curves. Sib. Mat. Zh. 7, 1206–1208.
[104] Schuster, T. & Quinto, E. T. 2005 On a regularization scheme for linear
operators in distribution spaces with an application to the spherical
Radon transform. SIAM J. Appl. Math. 65 (4), 1369–1387.
[105] Stefanov, P. & Uhlmann, G. Integral geometry of tensor fields on a class
of non-simple Riemannian manifolds. Preprint arXiv:math/0601178.
[106] Strichartz, Robert S. 2003 A Guide to Distribution Theory and Fourier
Transforms. World. Sci.
[107] Tam, A. C. 1986 Applications of photoacoustic sensing techniques. Rev.
Mod. Phys. 58, no. 2, 381–431.
[108] Tataru, D. 1995 Unique continuation for solutions to PDEs; between
Hörmander’s theorem and Holmgren’s theorem. Comm. PDE 20, 814–
http://arxiv.org/abs/math/0601178
[109] Tuchin, V. V. (Editor) 2002 Handbook of Optical Biomedical Diagnos-
tics. SPIE, Bellingham, WA.
[110] Vainberg, B. 1975 The short-wave asymptotic behavior of the solutions
of stationary problems, and the asymptotic behavior as t → ∞ of the
solutions of nonstationary problems. Russian Math. Surveys, 30, no. 2,
1–58.
[111] Vainberg, B. 1982. Asymptotics methods in the Equations of Mathe-
matical Physics. (Gordon & Breach)
[112] Vo-Dinh, T. (Editor) 2003 Biomedical Photonics Handbook. CRC, Boca
Raton, FL.
[113] Wang, L. V. &Wu, H. 2007 Biomedical Optics. Principles and Imaging.
Wiley-Interscience.
[114] Wang, X., Pang, Y., Ku, G., Xie, X., Stoica, G. & Wang, L. 2003
Noninvasive laser-induced photoacoustic tomography for structural and
functional in vivo imaging of the brain. Nature Biotechnology, 21, no.
7, 803–806.
[115] Xu, M. & Wang, L.-H. V. 2002 Time-domain reconstruction for ther-
moacoustic tomography in a spherical geometry. IEEE Trans. Med.
Imag. 21, 814–822.
[116] Xu, M. & Wang, L.-H. V. 2005 Universal back-projection algorithm for
photoacoustic computed tomography. Phys. Rev. E 71, 016706.
[117] Xu, M. & Wang, L.-H. V. 2006 Photoacoustic imaging in biomedicine.
Review of Scientific Instruments 77, 041101-01 – 041101-22.
[118] Xu, Y., Feng, D. & Wang, L.-H. V. 2002 Exact frequency-domain re-
construction for thermoacoustic tomography: I. Planar geometry. IEEE
Trans. Med. Imag. 21, 823–828.
[119] Xu, Y., Xu, M. & Wang, L.-H. V. 2002 Exact frequency-domain re-
construction for thermoacoustic tomography: II. Cylindrical geometry.
IEEE Trans. Med. Imag. 21, 829–833.
[120] Xu, Y., Wang, L., Ambartsoumian, G. & Kuchment, P. 2004 Recon-
structions in limited view thermoacoustic tomography. Medical Physics,
31(4), 724–733.
[121] Zobin, N. 1993. Unpublished.
	Introduction
	Thermoacoustic tomography
	Mathematical model of TAT: wave equation and the spherical mean transform
	The wave equation model
	Spherical mean model
	Main mathematical problems of TAT
	Uniqueness of reconstruction
	Constant speed case
	Non-uniqueness sets in R2.
	Higher dimensions
	Relations to other areas of analysis
	Uniqueness in the case of a variable sound speed
	Reconstruction: formulas and examples
	Constant sound speed
	Inversion formulas
	Fourier expansion methods
	Filtered backprojection methods
	Series solutions for arbitrary geometries
	Time reversal (backpropagation) methods
	Examples of reconstructions and additional remarks about the inversion formulas
	Reconstruction in the variable speed case
	``Analytic'' inversions
	Backpropagation
	Eigenfunction expansions
	Partial data. ``Visible'' and ``invisible'' singularities
	Range conditions
	Concluding remarks
	Variations of the TAT procedure
	Planar and linear transducers
	Direct imaging techniques
	Using contrast agents
	Passive thermoacoustic imaging
	Uniqueness
	Sketch of the proof of Theorem 10
	Some open problems concerning uniqueness
	Inversion
	Stability
	Range

  The paper presents a survey of mathematical problems, techniques, and
challenges arising in the Thermoacoustic and Photoacoustic Tomography.

Introduction
Computerized tomography has had a huge impact on medical diagnostics.
Numerous methods of tomographic medical imaging have been developed
and are being developed (e.g., the “standard” X-ray, single-photon emission,
positron emission, ultrasound, magnetic resonance, electrical impedance, op-
tical) [55, 59, 75, 76, 77]. The designers of these modalities strive to in-
crease the image resolution and contrast, and at the same time to reduce
the costs and negative health effects of these techniques. However, these
goals are usually rather contradictory. For instance, some cheap and safe
methods with good contrast (like optical or electrical impedance tomogra-
phy) suffer from low resolution, while some high resolution methods (such
as ultrasound imaging) often do not provide good contrast. Recently re-
searchers have been developing novel hybrid methods that combine different
physical types of signals, in hope to alleviate the deficiencies of each of the
types, while taking advantage of their strengths. The most successful exam-
ple of such a combination is the Thermoacoustic Tomography (TAT)
∗Mathematics Department, Texas A& M University, College Station, TX 77843-3368,
USA. kuchment@math.tamu.edu
†Mathematics Department, University of Arizona, Tucson, AZ 77843-3368, USA.
leonk@math.arizona.edu
http://arxiv.org/abs/0704.0286v2
1[62]. Albeit not being a common feature in clinics yet, TAT scanners are
actively researched, developed and already manufactured, for instance by
OptoSonics, Inc. (http://www.optosonics.com/), founded by the pioneer of
TAT R. Kruger.
After a substantial effort, major breakthroughs have been achieved in the
last couple of years in the mathematical modeling of TAT. The aim of this
article is to survey this recent progress and to describe the relevant models,
mathematical problems, and reconstruction procedures arising in TAT, and
to provide references to numerous research publications on this topic.
The main thrust of this text is toward mathematical methods; consid-
erations of the text length, as well as authors’ background do not let us
discuss in any detail industrial and physical set-ups and parameters of the
TAT technique, and limitations of the corresponding mathematical models.
Fortunately, the excellent recent surveys by M. Xu and L.-H. V. Wang [117]
and by A. A. Oraevsky and A. A. Karabutov [87, 88] accomplish all of these
tasks, and thus the reader is advised to consult with them for all such is-
sues (see also the recent textbook [113]). On the other hand, in spite of
the significant recent progress in mathematics of TAT, there is no compre-
hensive survey text addressing in details the relevant mathematical issues,
although the surveys [88, 117] do mention some mathematical reconstruction
techniques.
The structure of the paper is a follows: Section 2 contains a brief descrip-
tion of the TAT procedure. The next Section 3 provides the mathematical
formulation of the TAT problem. In general, it is formulated as an inverse
problem for the wave equation. However, in the case of the constant sound
speed, it can be also described in terms of a spherical mean operator (a
spherical analog of the Radon transform). The section also contains the list
of natural questions to be addressed concerning this model. These issues are
addressed then one by one in the following sections. In particular, Section
4 discusses uniqueness of reconstruction, i.e. the question of whether the
data collected in TAT is sufficient for recovery of the information of interest.
Albeit, for all practical purposes this issue is resolved in Corollary 2, we pro-
vide an additional discussion of unresolved uniqueness problems, which are
probably of more academic interest. Section 5 addresses inversion formulas
and algorithms. In Section 6 effects of having only partial data are discussed.
1TAT is also called Photoacoustic (PAT) or Optoacoustic (OAT) Tomography and is
sometimes abbreviated as TCT, which stands for Thermoacoustic Computed Tomography
http://www.optosonics.com/
Section 7 contains results concerning the so called range conditions, i.e. the
conditions that all ideal data must satisfy. Section 8 provides additional
remarks and discussions of the issues raised in the previous sections. The
paper ends with an Acknowledgments section and bibliography. Concerning
the latter, we need to mention that the engineering and biomedical literature
on TAT is rather vast and no attempt has been made in this text to create
a comprehensive bibliography of the topic from the engineering prospective.
The references in [87, 88, 109, 112, 117] to a large extent fill this gap. The
authors, however, have tried to present a sufficiently complete review of the
existing literature on mathematics of TAT.
2 Thermoacoustic tomography
In TAT, a short duration EM pulse is sent through a biological object (e.g.,
woman’s breast in mammography) with the aim of triggering a thermoa-
coustic response in the tissue. As it is explained in [117], the radiofrequency
(RF) and the visible light frequency ranges are currently considered to be the
most suitable for this purpose. Since mathematics works exactly the same
way in both of these frequency ranges, we will not make such distinction and
will be talking about just “an EM pulse”. E.g., in Figure 1 a microwave
pulse is assumed. In most cases the pulse is spatially wide, so that the
Figure 1: The TAT procedure.
whole object is more or less uniformly irradiated. Some part of EM energy
is absorbed throughout the object. The amount of energy absorbed at a
location strongly depends on local biological properties of the cells. Oxygen
saturation, concentration of hemoglobin, density of the microvascular net-
work (angiogenesis), ionic conductivity, and water content are among the
parameters that influence the absorption strongly [117]. Thus, if the energy
absorption distribution function f(x) were known, it would provide a great
diagnostic tool. For instance, it could be useful for detecting cancerous cells
that absorb several times more energy in the RF range than the healthy
ones [62, 88, 115, 117]. However, as an imaging tool neither RF waves, nor
visual light alone would provide acceptable resolution. In the RF case, this
is due to the long wave length. One can use shorter microwaves, but this
will be at the expense of the penetration depth. In the optical region, the
problem is with the multiple scattering of light. So, a different mechanism,
the so called Photoacoustic Effect [46, 107, 113, 117], is used to image f(x).
Namely, the EM energy absorption results in thermoelastic expansion and
thus in a pressure wave p(x, t) (an ultrasound signal) that can be measured
by transducers placed around the object. Now one can attempt to recover the
function f(x) (the image) from the measured data p(x, t). Such a measuring
scheme, utilizing two types of waves, brings about the high resolution of the
ultrasound diagnostics and the high contrast of EM waves. It overcomes the
adverse effect of the low contrast of ultrasound with respect to soft tissue.
In fact, such a low contrast is a good thing here, allowing one to assume in
the first approximation that the sound speed is constant. This often used
approximation is not always appropriate, but it is the most studied case at
the moment. Later on in this text we will describe some initial considerations
of the variable sound speed case, following [4, 57].
For this TAT method (and in particular, for the mathematical model
described below) to work, several conditions must be met. For instance, the
time duration of the EM pulse must be shorter than the time it takes the
sound wave to traverse the smallest feature that needs to be reconstructed.
The ultrasound detector must be able to resolve the time scale of the duration
of the EM pulse. On the other hand, the transducer must be also able to
detect much lower frequencies. Thus, one needs to have extra-wide-band
transducers, and these are currently available. One can find the technical
discussion of all these issues, for instance, in [88, 117]. In this text we will
assume that all these conditions are met and thus the mathematical models
described are applicable.
In the next section we present a mathematical description of the relation
between f(x) and p(x, t) (similar mathematical problems arise in sonar [73]
and radar [81] imaging, as well as in geophysics [27]).
3 Mathematical model of TAT:
wave equation and the spherical mean trans-
3.1 The wave equation model
We assume that the ultrasound speed at location x is equal to c(x). Then,
modulo some constant coefficients that we will assume all to be equal to 1,
the pressure wave p(x, t) satisfies the following problem for the standard wave
equation [28, 107, 115]:
ptt = c
2(x)∆xp, t ≥ 0, x ∈ R
p(x, 0) = f(x),
pt(x, 0) = 0
The goal is to find, using the data measured by transducers, the initial value
f(x) at t = 0 of the solution p(x, t).
In order to formalize what data is in fact measured, one needs to specify
what kind of transducers is used, as well as the geometry of the measurement.
By the geometry of the measurement we mean the distribution of locations
of transducers used to collect the data.
We briefly describe here the commonly considered measurement proce-
dure, which uses point detectors. Line and planar detectors have also been
suggested (see Section 8.1.1). It is too early to judge which one of them
will become most successful, but the one using point transducers has been
more thoroughly studied mathematically and experimentally, and thus will
be mostly addressed in this article. In this case, the transducers are assumed
to be point-like, i.e. of sufficiently small dimension. A transducer at time
t measures the average pressure over its surface at this time, which for the
small size of the transducer can be assumed to be just the value of p(y, t) at
the location y of the transducer. Dimension count shows immediately that in
order to have enough data for reconstruction of the function f(x), one needs
to collect data from the transducers’ locations y running over a surface S in
3. Thus, the data at the experimentalist’s disposal is the function g(y, t)
that coincides with the restriction of p(x, t) to the set of points y ∈ S.
Taking into account that the measurements produce the values g(y, t) of
the pressure p(x, t) of (1) on S × R+, the set of equations (1) extends to
become
ptt = c
2(x)∆xp, t ≥ 0, x ∈ R
p(x, 0) = f(x),
pt(x, 0) = 0
p(y, t) = g(y, t), y ∈ S × R+
Figure 2: An illustration to (2).
The problem now reduces to finding the initial value f(x) in (2) from the
knowledge of the lateral data g(x, t) (see Figure 3.1). A person familiar with
PDEs might suspect first that there is something wrong with this problem,
since we seem to have insufficient data for the recovery of the solution of the
wave equation in a cylinder from the lateral values alone. This, however,
is an illusion, since in fact there is a significant additional restriction: the
solution holds in the whole space, not just inside the cylinder S × R+. We
will see soon that in most cases, the data is sufficient for recovery of f(x).
3.2 Spherical mean model
We now introduce an alternative formulation of the problem that works in
the constant speed case only. We will assume that the units are chosen in
such a way that c(x) = 1. The known Poisson-Kirchhoff formula [25, Ch.
VI, Section 13.2, Formula (15)] for the solution of (1) gives
p(x, t) = c
(t(Rf)(x, t)) , (3)
where
(Rf)(x, r) =
|y|=1
f(x+ ry)dA(y) (4)
is the spherical mean operator applied to the function f(x), and dA is the
normalized area element on the unit sphere in R3. Hence, knowledge of the
function g(x, t) for x ∈ S and all t ≥ 0 essentially means knowledge of the
spherical mean Rf(x, t) at all points (x, t) ∈ S × R+. One thus is lead to
studying the spherical mean operator R : f → Rf and in particular its
restriction RS to the points x ∈ S only (these are the points where we place
transducers):
RSf(x, t) =
|y|=1
f(x+ ty)dA(y), x ∈ S, t ≥ 0. (5)
This explains why, in many works on TAT, the spherical mean operator has
been the model of choice. Albeit the (unrestricted) spherical mean operator
has been studied rather intensively and for a long time (e.g., [17, 25, 58]),
its version RS with the centers restricted to a subset S appears to have been
studied since early 1990s only [1]-[14], [16, 26, 30, 31, 33, 32, 34, 35, 36, 39,
41, 63, 64, 68, 69, 71, 72, 73, 77, 80, 82, 83, 89, 90, 91, 95, 96, 97, 104, 121]
and offers quite a few new and often hard questions.
In what follows, we will alternate between these two (PDE and integral
geometry) interpretations of the TAT model, since each of them has its own
advantages.
3.3 Main mathematical problems of TAT
We now formulate the typical list of problems one would like to address in
order to implement the TAT reconstruction.
1. For which sets S ∈ R3 is the data collected by transducers placed along
S sufficient for unique reconstruction of f? In terms of the spheri-
cal mean operator, the question is whether RS has zero kernel on an
appropriate class of functions, say continuous with compact supports.
2. If the data collected from S is sufficient, what are inversion formulas
and algorithms?
3. How stable is the inversion?
4. What happens if the data is “incomplete”?
5. What is the space of all possible “ideal” data g(t, y) collected on a sur-
face S? Mathematically (and in the constant sound speed case) it is
the question of describing the range of the operator RS in appropriate
function spaces. This question might seem to be unusual (for instance,
to people used to partial differential equations), but in tomography im-
portance of knowing the range of Radon type transforms is well known.
Such information is used to improve inversion algorithms, complete in-
complete data, discover and compensate for certain data errors, etc.
(e.g., [30, 38, 39, 40, 53, 54, 55, 76, 77, 90]).
4 Uniqueness of reconstruction
Many of the problems of interest to TAT can be formulated in any dimension
d, albeit the practical dimensions are only d = 3 and d = 2. We will consider
an arbitrary dimension d whenever we see this suitable.
Let S ⊂ Rd be the set of locations of transducers and f be a compactly
supported function (one can show that for purposes of uniqueness of recon-
struction problem, one can always assume that f is smooth [7]). Does the
absence of the signal on the transducers, i.e. g(t, y) = 0 for all t and y in
S, imply that f = 0? If the answer is a “yes,” we call S - a uniqueness
set, otherwise a non-uniqueness set. In other words, in terms of TAT, the
uniqueness sets are those that distributing transducers along them provides
enough data for unique reconstruction of the function f(x).
In terms of the wave equation, uniqueness sets are the sets of complete
observability, i.e. such that observing the motion on this set only, one gets
enough information to reconstruct the whole oscillation. In terms of the
spherical mean operator, the question is of whether the equality RSf = 0
implies that f = 0.
We will address this problem for the constant sound speed case first.
4.1 Constant speed case
As it has been discussed, the dimension count makes it clear that S must
be (d − 1)-dimensional, i.e. a surface in 3D or a curve in 2D. We will
also see that most of such surfaces are “good”, i.e. are uniqueness ones (or,
in other words, provide enough information for reconstruction). Thus, we
should rather discuss the problem of describing the “bad”, non-uniqueness
sets. The following simple statement is very important and not immediately
obvious.
Lemma 1 [7, 71, 72, 121] Any non-uniqueness set S is a set of zeros of a
(non-trivial) harmonic polynomial. In particular,
1. If there is no non-zero polynomial vanishing on S, then S is a unique-
ness set.
2. If there is no non-zero harmonic function vanishing on S, then S is a
uniqueness set.
The proof of this lemma is very simple. It works under the assumption
of exponential decay of the function f(x), not necessarily of compactness of
its support. It also introduces some polynomials that play significant role in
the whole analysis of the spherical mean operator RS.
Let k ≥ 0 be an integer. Consider the convolution
Qk(x) = |x|
2k ∗ f(x) =
|x− y|2kf(y)dy. (6)
This is clearly a polynomial of degree at most 2k. Rewriting the integral in
polar coordinates centered at x and using radiality of |x − y|, one sees that
Qk(x) is determined if we know the values Rf(x, t) of the spherical mean of
f centered at x:
Qk(x) = cd
t2k+d−1Rf(x, t)dt.
In particular, If RSf ≡ 0, then each polynomial Qk vanishes on S.
Another observation that is easy to justify is that if the function f is
exponentially decaying (e.g., is compactly supported), then if all polynomials
Qk vanish identically, the function itself must be equal to zero. (This is not
necessarily true anymore if f and its derivatives decay only faster than any
power, rather than exponentially.)
Thus, we conclude that if f is not identically equal to zero, then there is at
least one non-zero polynomial Qk. Since, as we discussed, equality RSf = 0
implies that Qk|S = 0, we conclude that S must be algebraic.
Now notice the following simple to verify equality (with a non-zero con-
stant ck):
∆Qk = ckQk−1, (7)
where ∆ is the Laplace operator. This implies that the lowest k non-zero
polynomial Qk is harmonic. Since Qk|S = 0, this proves the lemma.
Consider now the case when S is a closed (hyper-)surface (i.e., the bound-
ary of a bounded domain). Since, as it is well known, there is no non-zero
harmonic function in the domain that would vanish at the boundary (the
spectrum of the Dirichlet Laplace operator is strictly positive), we conclude
that such S is a uniqueness set for harmonic polynomials. Thus, we get the
following important
Corollary 2 [7, 63] Any closed surface is uniqueness set for the spherical
mean Radon transform.
An older alternative proof of this corollary provides an additional insight
into the problem. We thus sketch it here. Let us assume for simplicity
that the dimension d ≥ 3 is odd (even dimensions require a little bit more
work). Suppose that the closed surface S remains stationary (nodal) for the
oscillation described by (1). Since the oscillation is unconstrained and the
initial perturbation is compactly supported, after a finite time, the interior
of S will become stationary. On the other hand, we can think that S is
fixed (since it is not moving anyway). Then, the energy inside S must stay
constant. This is the contradiction that proves the statement of Corollary 2.
We will see in the next Section that the same method works in some cases
of variable sound speed, providing the needed uniqueness of reconstruction
result.
This corollary resolves the uniqueness problems for most practically used
geometries. It fails, however, if f does not decay sufficiently fast (see [3],
where it is shown in which Lp(Rd) classes of functions f(x) closed surfaces
remain uniqueness sets).
It also provides uniqueness for some “limited data” problems. For in-
stance, if S is an open (even tiny) piece of an analytic closed surface Σ, it
suffices. Indeed, if it did not, then it would be a part of an algebraic non-
uniqueness surface. Uniqueness of analytic continuation would show then
that the whole Σ is a non-uniqueness set, which we know to be incorrect.
This result, however, does not say that it would be practical to reconstruct
using observations from a tiny S. We will see later that this would not lead
to a satisfactory reconstructions, due to instabilities.
A geometry sometimes used is the planar one, i.e. detectors are placed
along a plane S (line in the 2D). In this case, there is no uniqueness of
reconstruction when the sound speed is constant. Indeed, if f(x) is odd with
respect to S, then clearly all measured data g(t, y) will vanish. However, it is
well known [25, 58] that functions even with respect to S can be recovered.
What saves the day in TAT is that the object to be imaged is located on one
side of S. Then, extending f(x) as an even function with respect to S, one
can still recover it from the data.
Although, for all practical purposes the uniqueness of reconstruction prob-
lem is essentially resolved by the Corollary 2, the complete understanding of
uniqueness problem has not been achieved yet. Thus, we include below some
known theoretical results and open problems.
4.1.1 Non-uniqueness sets in R2.
In this Section, we follow the results and exposition of [7, 71, 72] in discussing
uniqueness sets in 2D. What are simple examples of non-uniqueness sets? As
we have already mentioned, any line S (or a hyperplane in higher dimensions)
is a non-uniqueness set, since any function f odd with respect to S will clearly
produce no signal: RSf = 0. Analogously, consider a Coxeter system ΣN of
N lines passing through a point and forming equal angles (see Fig. 3).
Figure 3: Coxeter cross ΣN .
Choosing the intersection point as the pole and expanding functions into
Fourier series with respect to the polar angle, it is easy to discover existence
of an infinite dimensional space of functions that are odd with respect to
each of the N lines. Thus, such a cross ΣN is also a non-uniqueness set. Less
obviously, one can use the infinite dimensional freedom just mentioned to
add any finite set Φ of points still preserving non-uniqueness. The following
major and very non-trivial result was conjectured in [71, 72] and proven in
[7]. It shows that there are no other bad sets S besides the ones we have just
discovered:
Theorem 3 A set S ⊂ R2 is a non-uniqueness set for the spherical mean
transform in the space of compactly supported functions, if and only if
S ⊂ ωΣN ∪ Φ,
where ΣN is a Coxeter system of lines, ω is a rigid motion of the plane, and
Φ is a finite set.
A sketch of a rather intricate proof of this result is provided in Section
4.1.2 Higher dimensions
Here we present a believable conjecture of how the result should look like in
higher dimensions.
Conjecture 4 [7]A set S ⊂ Rd is a non-uniqueness set if and only if S ⊂
ωΣ ∪ Φ, where Σ is the surface of zeros of a homogeneous harmonic polyno-
mial, ω is a rigid motion of Rd, and Φ is an algebraic surface of dimension
at most d− 2.
Figure 4: A picture of a 3-dimensional non-uniqueness set.
The progress towards proving this conjecture has been slow, albeit some
partial cases have been treated ([1]-[12]). E.g., in some cases one can prove
that S is a ruled surface (i.e., consists of lines), but proving that these lines
(rules) pass through a common point remains a challenge. It is known,
though, that both the zero sets of homogeneous harmonic polynomials and
algebraic subsets of dimension at most d − 2 are non-uniqueness sets [2, 7],
and thus one should avoid using them as placements of transducers for TAT.
4.1.3 Relations to other areas of analysis
The problem of injectivity of RS has relations to a wide variety of areas of
analysis (see [1, 7] for many examples). In particular, the following interpre-
tation is important:
Theorem 5 [7, 63] The following statements are equivalent:
1. S ⊂ Rd is a non-uniqueness set for the spherical mean operator.
2. S is a nodal set for the wave equation, i.e. there exists a non-zero
compactly supported f such that the solution of the wave propagation
problem 
= ∆u,
u(x, 0) = 0,
ut(x, 0) = f(x)
vanishes on S for any moment of time.
3. S is a nodal set for the heat equation, i.e. there exists a non-zero
compactly supported f such that the solution of the problem
= ∆u,
u(x, 0) = f(x)
vanishes on S for any moment of time.
The interpretation in terms of the wave equation provides important PDE
tools and insights, which have lead to a recent progress [33, 12] (albeit it
has not lead yet to a complete alternative proof of Theorem 3). The rough
idea, originally introduced in [33], is that if S is a nodal set, then it might be
considered as the fixed boundary. In this case, the signals must go around S.
However, in fact, there is no obstacle, so signals can propagate along straight
lines. Thus, in order to avoid discrepancies in arrival times, S must be very
special. One can find details in [33] and in [12].
4.2 Uniqueness in the case of a variable sound speed
It is shown in [35, Theorem 4] that uniqueness of reconstruction also holds in
the case of a smoothly varying (strictly positive) sound speed, if the source
function f(x) is completely surrounded by the observation surface S (in other
words, if there is no US signal coming from outside of S). The proof uses
the celebrated unique continuation result by D. Tataru [108].
One can also establish uniqueness of reconstruction in the case of the
source not necessarily completely surrounded by S. However, here we need
to impose an additional non-trapping condition on the sound speed. We
assume that the sound speed is strictly positive c(x) > c > 0 and such that
c(x)− 1 has compact support, i.e. c(x) = 1 for large x.
Consider the Hamiltonian system in R2nx,ξ with the Hamiltonian H =
c2(x)
|ξ|2: 
x′t =
= c2(x)ξ
ξ′t = −
∇ (c2(x)) |ξ|2
x|t=0 = x0, ξ|t=0 = ξ0.
The solutions of this system are called bicharacteristics and their projections
into Rnx are rays.
We will assume that the following non-trapping condition holds:
all rays (with ξ0 6= 0) tend to infinity when t→ ∞.
Theorem 6 [4] Under the non-trapping conditions formulated above, com-
pactly supported function f(x) is uniquely determined by the data g measured
on S for all times. (No assumption of f being supported inside S is imposed.)
One should mention that ray trapping can occur for some sound speed
profiles. For instance, if c(x) = |x| for some range r1 < |x| < r2, then there
are rays trapped in this spherical shell. We are not sure what happens in
this case to the uniqueness of reconstruction statement of Theorem 6 and
inversion formula of Theorem 7.
5 Reconstruction: formulas and examples
Here we will address the procedures of actual reconstruction of the source
f(x) from the data g(t, y) measured by transducers.
5.1 Constant sound speed
We assume here that the sound speed is constant and normalized to be equal
to 1.
5.1.1 Inversion formulas
Before we move to our case of interest, which is spheres centered on a closed
surface S surrounding the object to be imaged, we briefly refer to related
but somewhat different work. Namely, the problem of recovering functions
from integrals over spheres centered on a (hyper)plane S has attracted a
lot of attention over the years. Albeit, as it has been mentioned before,
there is no uniqueness in this case (functions odd with respect to S are
annihilated), even functions can be recovered. Thus also functions supported
on one side of the plane can be as well, by means of their even extension.
Many explicit inversion formulas and procedures have been obtained for this
situation [16, 26, 31, 39, 41, 60, 77, 80, 89, 90, 101, 102, 103]. We will not
provide any details here, since this acquisition geometry is not very useful.
In particular, this is due to “invisibility” of some parts of the interfaces,
see Section 6, which arises from truncating the plane. The same problem
is encountered with some other unbounded acquisition surfaces, such as a
surface of an “infinitely” long cylinder.
Thus, it is more practical to place transducers along a closed surface
surrounding the object. The simplest surface of this type is a sphere.
5.1.2 Fourier expansion methods
Let us assume that S is the unit sphere in Rn. We would like to reconstruct
a function f(x) supported inside S from the known values of its spherical
integrals g(y, r) with the centers on S:
g(y, r) =
f(y + rω)rn−1dω, y ∈ S.
The first inversion procedures for the case of spherical acquisition were de-
scribed in [82] in 2D and in [83] in 3D. These solutions were obtained by
harmonic decomposition of the measured data and the sought function, and
by equating coefficients of the corresponding Fourier series.
In particular, the 2-D algorithm of [82] is based on the Fourier decompo-
sition of f and g in angular variables:
f(x) =
fk(ρ)e
ikϕ, x = (ρ cos(ϕ), ρ sin(ϕ)) (9)
g(y(θ), r) =
gm(r)e
ikθ, y = (R cos(θ), R sin(θ)).
Following [82] we consider the Hankel transform ĝm,J(λ) of the Fourier coef-
ficients gm(r) (divided by 2πr)
ĝm,J(λ) =
gm(r)J0(λr)dr = H0
gm(r)
. (10)
To simplify the presentation we introduce the convolution GJ(λ, y) of the
sought function with the Bessel function J0(λ|x− y|).
GJ(λ, y) =
f(x)J0(λ|x− y|)dx, (11)
One can notice that ĝm,J(λ) are the Fourier coefficients of GJ(λ, y) in θ:
ĝm,J(λ) =
GJ(λ, y)e
−imθdθ. (12)
Now coefficients fm(ρ) can be recovered from gm(r) by application of the
addition theorem for the Bessel function J0(λ|x− y|):
J0(λ|x− y|) =
Jm(λ|x|)Jm(λ|y|)e
−im(ϕ−θ). (13)
Indeed, by substituting equations (9) and (13) into (11), and (11) into (12)
one obtains
ĝm,J(λ) = 2πJm(λ|R|)
fm(ρ)Jm(λρ)ρdρ = Hm(fm(ρ)),
where Hm is the m-th order Hankel transform. Since the latter transform
is self-invertible, the coefficients fm(ρ) can be recovered by the following
formula
fm(ρ) = Hm
ĝm,J(λ)
Jm(λ|R|)
Jm(λ|R|)
gm(r)
, (14)
which is the main result of [82]. Function f(x) can now be reconstructed by
summing series (9).
Note that the above method requires a division of the Hankel transform
of the measured data by Bessel functions Jm that have infinitely many zeros.
Theoretically, there is no problem; the Hankel transform H0
gm(r)
has to
have zeros that would cancel those in the denominator. However, since the
measured data always contain some error, the exact cancelation is not likely
to happen, and one needs a sophisticated regularization scheme to keep the
total error bounded.
This problem can be avoided by replacing in (10) Bessel function J0 by
Hankel function H
ĝm,H(λ) =
gm(r)H
0 (λr)dr.
The addition theorem for H
0 (λ|x− y|) takes form
0 (λ|x− y|) =
Jm(λ|x|)H
m (λ|y|)e
−im(ϕ−θ),
and by proceeding as before one can obtain the following formula for fm(ρ):
fm(ρ) = Hm
ĝm,H(λ)
m (λ|R|)
m (λ|R|)
gm(r)H
0 (λr)dr
Unlike Jm, Hankel functions H
m (t) do not have zeros for all real values of t
and therefore problems with division by zeros do not arise in this amended
version of the method [82].
This derivation can be repeated in 3-D, with the exponentials eikθ replaced
by the spherical harmonics, and with cylindrical Bessel functions replaced by
their spherical counterparts. By doing this one will arrive at the Fourier
series method of [83]. Our use of Hankel function H
0 above is similar to
the way the authors of [83] utilized spherical Hankel function h
0 to avoid
the divisions by zero.
5.1.3 Filtered backprojection methods
The favorite way of inverting Radon transform for tomography purposes is
by using filtered backprojection type formulas, which involve filtration in
Fourier domain followed (or preceded) by a backprojection. In the case of
the set of spheres centered on a closed surface (e.g., sphere) S, one expects
such a formula to involve a filtration with respect to the radius variable and
then some integration over the set of spheres passing through the point of
interest. For quite a while, no such type formula had been discovered. This
did not prevent practitioners from reconstructions, since good approximate
inversion formulas (parametrices) could be developed, followed by an iterative
improvement of the reconstruction, see e.g. reconstruction procedures in
[114, 115, 118, 119, 120], and especially [96, 97].
The first set of exact inversion formulas of the filtered backprojection type
was discovered in [33]. These formulas were obtained only in odd dimensions.
Several different variations of such formulas (different in terms of exact order
of the filtration and backprojection steps) were developed. Let us denote by
g(p, r) = r2RSf the spherical integral, rather than the average, of f . Then
various versions of the 3D inversion formulas that reconstruct a function f(x)
supported inside S from its the spherical mean data RSf , read:
f(x) = − 1
g(y, |y − x|)dA(y),
f(x) = − 1
g(y, t)
) ∣∣∣∣∣
t=|y−x|
dA(y),
f(x) = − 1
g(y,t)
)) ∣∣∣∣∣
t=|y−x|
dA(y).
Recently, analogous formulas were obtained for even dimensions in [32]. De-
noting by g, as before the spherical integrals (rather than averages) of f , the
formulas of [32] in 2D look as follows:
f(x) =
g(y, t) log(t2 − |x− y|2) dt dl(y), (16)
f(x) =
g(y, t)
log(t2 − |x− y|2) dt dl(y), (17)
A different set of explicit inversion formulas that work in arbitrary dimensions
was presented in [69].
f(x) =
4(2π)n−1
n(y)h(y, |x− y|)dA(y). (18)
h(y, t) =
Y (λt)
J(λt′)g(y, t′)dt′
−J(λt)
Y (λt′)g(y, t′)dt′
λ2n−3dλ, (19)
J(t) =
Jn/2−1(t)
tn/2−1
, Y (t) =
Yn/2−1(t)
tn/2−1
Jn/2−1(t) and Yn/2−1(t) are respectively the Bessel and Neumann functions
of order n/2− 1, and n(y) is the vector of exterior normal to ∂B.
In 2-D equations (18), (19) can be simplified to yield the following recon-
struction formula
f(x) = −
g(y, t′)
|x− y|2 − t′2
 dl(y).
A similar simplification is also possible in 3D resulting in the formula
f(x) =
g(y, t)
) ∣∣∣∣∣
t=|y−x|
dA(y). (20)
Equation (20) is equivalent to one of the formulas derived in [116] for the 3D
case. It is interesting to notice that the “universal” formula of [116] holds for
all geometries when the backprojection type formulas are known: spherical,
cylindrical, and planar. It is not very likely that such explicit formulas would
be available for any closed surfaces S different from spheres (see a related
discussion in [15, 27]).
5.1.4 Series solutions for arbitrary geometries
Although, as we have just mentioned, we do not expect such explicit formulas
to be derived for non-spherical closed surfaces S, there is, however, a different
approach [70] that theoretically works for any closed S and that is practically
useful in some non-spherical geometries.
Let λ2m and um(x) be the eigenvalues and normalized eigenfunctions of
the Dirichlet Laplacian −∆ on the interior Ω of a closed surface S:
∆um(x) + λ
mum(x) = 0, x ∈ Ω, Ω ⊆ R
n, (21)
um(x) = 0, x ∈ S,
||um||
|um(x)|
2dx = 1.
As before, we would like to reconstruct a compactly supported function f(x)
from the known values of its spherical integrals g(y, r) with the centers on S:
g(y, r) =
f(y + rω)rn−1dω, y ∈ S.
We notice that um(x) is the solution of the Dirichlet problem for the Helmholtz
equation with zero boundary conditions and the wave number λm, and thus
it admits the Helmholtz representation
um(x) =
Φλm(|x− y|)
um(y)ds(y) x ∈ Ω, (22)
where Φλm(|x− y|) is a free-space rotationally invariant Green’s function of
the Helmholtz equation (21).
The eigenfunctions {um(x)}
0 form an orthonormal basis in L2(Ω). There-
fore, f(x) can be represented by the series
f(x) =
αmum(x) (23)
um(x)f(x)dx.
Since f(x) is C10 , series (23) converges pointwise. A reconstruction formula of
αm, and thus of f(x), will result if we substitute representation (22) into (23)
and interchange the order of integrations. Indeed, after a brief calculation
we will get
um(x)f(x)dx =
I(y, λm)
um(y)dA(x), (24)
where
I(y, λ) =
Φλ(|x− y|)f(x)dx. (25)
Certainly, the need to know the spectrum and eigenfunctions of the
Dirichlet Laplacian imposes a severe constraint on the surface S. However,
there are simple cases when the eigenfunctions are well known, and fast sum-
mation formulas for the corresponding series are available. Such is the case
of a cubic measuring surface S (see [70]); the eigenfunctions um are products
of sine functions
um(x) =
πm1x1
πm2x2
πm3x3
, (26)
where m = (m1, m2, m3), m1, m2, m3 ∈ N, and the eigenvalues are easily
found as well
λm = π
2|m|2/R2. (27)
Sum (23) is just a regular 3-D Fourier sine series easily computable by ap-
plication of the Fast Sine Fourier transform algorithm. The algorithmic
trick that allows one to calculate fast the coefficients αm consists in com-
puting first integrals (25) on a uniform mesh in λ. This is easily done by
a one-dimensional Fast Cosine Fourier transform algorithm, with Φλ(t) =
cos(λt)/t. The normal derivatives of um(x) are also products of sine func-
tions, this time two-dimensional ones. This, in turn, permits rapid evaluation
of integrals
I(y, λ) ∂
um(y)dA(x) for each mesh value of λ, and for each
one of the six faces ∂Ωi, i = 1, ..., 6 of the cube. Finally, the computation of
αm using equation (24) reduces to the interpolation in the spectral param-
eter λ, since the integrals in the right hand side of this equation have been
computed for the mesh values of this parameter (not for λm). Due to oscil-
latory nature of the integrals (25) a low order interpolation here would lead
to inaccurate reconstructions. Luckily, however, these integrals are analytic
functions of parameter λ (due to the finite support of g). Hence, high order
polynomial interpolation is applicable, and numerics yields very good results.
The algorithm we just described requires O(m3 logm) floating point op-
erations if the reconstruction is to be performed on an m×m×m Cartesian
grid, from comparably discretized data measured on a cubic surface. In prac-
tical terms, it yields reconstructions in the matter of several seconds on grids
with total number of nods exceeding a million [70].
5.1.5 Time reversal (backpropagation) methods
In the constant speed case, the following approach is possible in 3D: due
to the validity of the Huygens’ principle (i.e., the signal escapes from any
bounded domain in finite time), the pressure p(t, x) inside S will become
equal to zero for any time T larger than the time required to cross the domain
(i.e., time that it takes the sound to move along the diameter of S, which
for c = 1 equals the diameter). Thus, one can impose the zero conditions
on p(t, x) for t = T and solve the wave equation (2) back in time, using the
measured data g as the boundary values. The solution of this well posed
problem at t = 0 gives the desired source function f(x). Such methods have
been successfully implemented [22].
Although in 2D or in presence of sound speed variations, Huygens’ princi-
ple does not hold anymore, and thus the signal theoretically will stay forever,
one can find good approximate solutions using a similar approach [4, 18], see
discussions of the variable speed case below.
5.1.6 Examples of reconstructions and additional remarks about
the inversion formulas
• It is well known that different analytic inversion formulas in tomogra-
phy can behave differently in numerical implementation (e.g., in terms
of their stability), However, numerical implementation seems to show
that the analytic (backprojection type) formulas (15)-(20), in spite of
some of them being not equivalent, work equally well. See, for example
the results of an analytic formula reconstruction in 3D shown in Fig.
• It is worth noting that although formulas (15)-(16) and (18)-(20) will
yield identical results when applied to functions that can be represented
as the spherical mean Radon transform of a function supported inside
S, they are in general not equivalent when applied to functions with
larger supports. Simple examples (e.g., of f being the characteristic
Figure 5: A mathematical phantom in 3D (left) and its reconstruction using
an analytic inversion formula.
set of a large ball containing S) show that these two types of formulas
provide different reconstructions.
• An interesting observation is that backprojection formulas (15)-(20) do
not reconstruct the function f correctly inside the surface S, if f has
support reaching outside S. For instance, applying the reconstruction
formulas to the function RS(χ|x|≤3) leads to an incorrect reconstruction
of the value of f = χ|x|≤3 inside S = {|x| ≤ 1}. (Here by χV we denote
the characteristic function of the set V , i.e. it takes the value 1 in V
and zero outside. So, χ|x|≤3 is the characteristic function of the ball of
radius 3 centered at the origin.)
An another example: if one adds to the phantom shown in Fig. 5
two balls to the right of the surrounding sphere S, this leads to strong
artifacts, as seen on Fig. 6.
What is the reason for such a distortion? If one does not know in ad-
vance that f has support inside S, the backprojection formulas shown
before use insufficient information to recover a function with a larger
support, and thus uniqueness of reconstruction is lost. Then the for-
mulas misinterpret the data, wrongly assuming that they came form a
function supported inside S and thus reconstructing the function in-
correctly.
Notice that the series reconstruction of the preceding Section is free of
such problem. E.g., the reconstruction shown in Fig. 7 confirms this.
Figure 6: A perturbed reconstruction, due to presence of two additional balls
outside S (not shown on the picture).
Figure 7: In the phantom shown on the left, most disks are located outside
the square acquisition surface S indicated by the dotted line. This does not
perturb the reconstruction inside S (right).
5.2 Reconstruction in the variable speed case
We will assume here that the sound speed c(x) is smooth, positive, constant
for large x, and non-trapping. Although most analytic techniques we de-
scribed above do not work in the variable speed case, some formulas can be
derived and algorithms can be designed. This work is in a beginning stage
and the results described below most surely can and will be improved.
5.2.1 “Analytic” inversions
Let us denote by Ω the interior of the observation surface S, i.e. the area
where the object to be imaged is located. Consider in Ω the operator A =
−c2(x)∆ with zero Dirichlet conditions on the boundary S = ∂Ω. This
operator is self-adjoint, if considered in the weighted space L2(Ω; c−2(x)).
We also denote by E the operator of harmonic extension, which trans-
forms a function φ on S to a harmonic function on Ω which coincides with φ
on S.
The following result provides a formula for reconstructing f from the data
Theorem 7 [4] The function f(x) in (2) can be reconstructed in Ω as fol-
lows:
f(x) = (Eg|t=0)−
2 sin (τA
2 )E(gtt)(x, τ)dτ. (28)
The validity of this result hinges upon decay estimates for the solution (so
called local energy decay [29, 110, 111]), which hold under the non-trapping
condition. These estimates guarantee a qualified decay of the solution p(t, x)
inside any bounded region, e.g. in Ω, when time t increases. In odd dimen-
sions decay is exponential, but only polynomial in even dimensions. The
decay can be used instead of Huygens’ principle to solve the wave equation
backwards, starting at the infinite time. This leads to the formula (28).
Due to functions of the operator A being involved, it is not that clear how
explicit this formula can be made. For instance, it would be interesting to
see whether one can derive from (28) a backprojection inversion formula for
the case of a constant sound speed and S being a sphere (we have already
seen that such formulas are known).
5.2.2 Backpropagation
The exponential decay at large values of time can be used as follows: for a
sufficiently large T , one can assume that the solution is practically zero at
t = T . Thus, imposing zero initial conditions at t = T and solving in reverse
time direction, one arrives at t = 0 to an approximation of f(x) [18].
5.2.3 Eigenfunction expansions
One natural way to try to use the formula (28) is to use eigenfunction expan-
sion of the operator A in Ω (assuming that such expansion is known). This
immediately leads to the following result:
Theorem 8 Under the same conditions on the sound speed as before, func-
tion f(x) can be reconstructed inside Ω from the data g in (2), as the following
L2(B)-convergent series:
f(x) =
fkψk(x), (29)
where the Fourier coefficients fk can be recovered using one of the following
formulas:
fk = λ
k gk(0)− λ
sin (λkt)g
k(t)dt,
fk = λ
k gk(0) + λ
cos (λkt)g
k(t)dt, or
fk = −λ
sin (λkt)gk(t)dt = −λ
sin (λkt)g(x, t)
(x)dxdt,
gk(t) =
g(x, t)
(x)dx.
Here ν denotes the external normal to S.
One notices that this is a generalization to the variable sound speed case
of the expansion method of [70], discussed in Section 5.1.4. An interesting
feature is that, unlike in [70], we do not need to know the whole space Green’s
function for A (which is certainly not known).
It is not clear yet how feasible numerical implementation of this approach
could be.
6 Partial data. “Visible” and “invisible” sin-
gularities
Uniqueness of reconstruction does not imply practical recoverability, since the
reconstruction procedure might be severely unstable. This is well known to
be the case, for instance, in incomplete data situations in X-ray tomography,
and even for complete data problems in some imaging modalities, such as
the electrical impedance tomography [64, 68, 76, 77].
In order to describe the results below, we need to explain the notion
of the wave front set WF (f) of a function f(x). This set carries detailed
information on singularities of f(x). It consists of pairs (x, ξ) of a point x
in space and a wave vector (Fourier domain variable) ξ 6= 0. It is easier to
say what it means that a point (x0, ξ0) is not in the wave front set WF (f).
This means that one can smoothly cut-off f to zero at a small distance from
x0 in such a way that the Fourier transform φ̂f(ξ) of the resulting function
φ(x)f(x) decays faster than any power of ξ in directions that are close to
the direction of ξ0. We remind the reader that if this Fourier transform
decays that way in all directions, then f(x) is smooth near the point x0. So,
the wave front set contains pairs (x0, ξ0) such that f is not smooth near x0,
and ξo indicates why it is not: the Fourier transform does not decay well in
this direction. For instance, if f(x) consists of two smooth pieces joined non-
smoothly across a smooth interface Σ, then WF (f) contains pairs (x, ξ) such
that x is in Σ and ξ is normal to Σ at x. One can find simple introduction
to the notions of microlocal analysis, such as the wave front set, for instance
in [106].
Analysis done in [99] for the constant speed case (equivalently, for the
spherical mean transform RS), showed which parts of the wave front (and
thus singularities) of a function f can be recovered from its partial X-ray
data. An analog of this result also holds for the spherical mean transform
RS [73] (see also [120] for a practical discussion). We formulate it below in
an imprecise form (see [73] for precise formulation).
Theorem 9 [73] A wavefront set point (x, ξ) of f is “stably recoverable”
from RSf if and only if there is a circle (sphere in higher dimensions) centered
on S, passing through x, and normal to ξ at this point.
As we have already mentioned, this result does not exactly hold the way it is
formulated and needs to include some precise conditions (see [73, Theorem
3]). The statement is, for instance, correct if S is a smooth hypersurface and
the support of f lies on one side of the tangent plane to S at the center of
the sphere mentioned in the theorem.
Talking about jump singularities only (i.e., interfaces between smooth
regions inside the object to be imaged), this result says that in order for a
piece of the interface to be stably recoverable (dubbed “visible”), one should
have for each point of this interface, a sphere centered at S and tangent to
the interface at this point. Otherwise, the interface will be blurred away
(even if there is a uniqueness of reconstruction theorem). The reason is that
if all spheres of integration are transversal to the interface, the integration
smoothes off the singularity, and therefore its recovery becomes highly unsta-
ble (numerically, one has to deal with inversion of a matrix with exponentially
fast decaying singular values). The Figure 8 below shows an example of an
incomplete data reconstruction from spherical mean data. One sees clearly
the effect of disappearance of the parts of the boundaries that are not touched
tangentially by circles centered at transducers’ locations.
Figure 8: Effect of incomplete data: the phantom (left) and its incomplete
data reconstruction. The transducers were located along a 180o circular arc
(the left half of a large circle surrounding the squares).
7 Range conditions
As it has already been mentioned, the space of functions g(t, y) that could
arise as exact data measured by transducers (i.e., the range of the data),
is very small (of infinite codimension in the spaces of all functions of t >
0, y ∈ S). Knowing this space (range) is useful for many theoretical and
practical purposes (reconstruction algorithms, error corrections, incomplete
data completion, etc.), and thus has attracted a lot of attention (e.g., [30,
38, 39, 40, 53, 54, 64, 66, 67, 68, 74, 76, 77, 78, 90, 100].
For instance, for the standard Radon transform
f(x) → g(s, ω) =
x·ω=s
f(x)dx, |ω| = 1,
the range conditions on g(s, ω) are:
1. evenness: g(−s,−ω) = g(s, ω)
2. moment conditions: for any integer k ≥ 0, the kth moment
Gk(ω) =
skg(ω, s)ds
extends from the unit circle of vectors ω to a homogeneous polynomial
of degree k in ω.
The evenness condition is obviously necessary and is kind of “trivial”. It
seems that the only non-trivial conditions are the moment ones. However,
here the standard Radon transform misleads us, as it often happens. In fact,
for more general transforms of Radon type it is often easy (or easier) to find
analogs of the moment conditions, while analogs of the evenness conditions
are often elusive (see [64, 66, 67, 76, 77, 84] devoted to the case of SPECT
(single photon emission tomography)). The same happens in TAT.
Let us deal first with the case of a constant sound speed, when one can
think of the spherical mean transform RS instead of the wave equation model.
An analog of the moment conditions was already present implicitly (without
saying that these were range conditions) in [71, 72, 7] and explicitly formu-
lated as such in [95]. Indeed, our discussion in Section 4 of the polynomials
Qk provides the following conditions of the moment type:
Moment conditions [7, 71, 72, 95] on data g(p, r) = RSf(p, r) look as
follows: for any integer k ≥ 0, the moment
Mk(ω) =
r2k+d−1g(p, r)dr
can be extended from S to a (non-homogeneous) polynomial Qk(x) of degree
at most 2k.
These conditions, however, are incomplete, and in fact infinitely many
others, which play the role of an analog of evenness, need to be added.
Complete range descriptions for RS when S is a circle in 2D were discov-
ered in [13] and then in odd dimensions in [34]. They were then extended
to any dimension and interpreted in several different ways in [6]. These
conditions happen to be intimately related to PDEs and spectral theory.
Figure 9:
In order to describe these conditions, we need to introduce some notations.
Let B be the unit ball in Rd, S - the unit sphere, and C - the cylinder B×[0, 2]
(see Fig. 9).
We introduce the spherical mean operator RS as before:
RSf(x, t) =
|y|=1
f(x+ ty)dA(y), x ∈ S.
Several different range descriptions for RS were provided in [6], out of
which we only show a few:
Theorem 10 [6] The following three statements are equivalent:
1. The function g ∈ C∞0 (S × [0, 2]) is representable as RSf for some
f ∈ C∞0 (B). (In other words, g represents an ideal (free of errors) set
of TAT data.)
2. (a) The moment conditions are satisfied.
(b) Let −λ2 be any eigenvalue of the Laplace operator in B with zero
Dirichlet conditions and ψλ be the corresponding eigenfunction.
Then the following orthogonality condition is satisfied:
S×[0,2]
g(x, t)∂νψλ(x)jn/2−1(λt)t
n−1dxdt = 0. (31)
Here jp(z) = cp
Jp(z)
is the so called spherical Bessel function.
3. (a) The moment conditions are satisfied.
(b) Let ĝ(x, λ) =
g(x, t)jn/2−1(λt)t
n−1dt. Then, for any m ∈ Z, the
mth spherical harmonic term ĝm(x, λ) of ĝ(x, λ) vanishes at all
zeros λ 6= 0 of Bessel function Jm+n/2−1(λ).
Remark 11 [6]
1. In odd dimensions, moment conditions are not necessary, and thus con-
ditions 2(b) or 3(b) suffice. (A similar earlier result was established for
a related transform in [34].)
2. The range conditions (2) of the previous Theorem are also necessary
when S is the boundary of any bounded domain, not necessarily a
sphere.
3. An analog of these conditions can be derived for a variable sound speed
(without non-trapping conditions imposed).
8 Concluding remarks
8.1 Variations of the TAT procedure
8.1.1 Planar and linear transducers
Assuming that transducers are point-like, is clearly an approximation, and
in fact a transducer measures the average pressure over its area. It has
been rightfully claimed that the point approximation for transducers should
lead to some blurring in the reconstructions. This, as well as intricacies of
reconstructions from the data obtained by point transducers, triggered recent
proposals for different types of transducers (see [20, 21], [47]-[52], [92, 93]).
In these papers, it was suggested to use either planar, or line detectors.
In the first case [47], the detectors are assumed to be large and planar,
ideally assumed to be approximations of infinite planes that are placed tan-
gentially to a sphere containing the object. Thus, the data one collects is
the integrals of the pressure over these planes, for all values of t > 0. If one
takes the standard 3D Radon transform of the pressure p(x, t) with respect
to x:
P (x, t) 7→ q(s, t, ω) =
x·ω=s
p(x, t)dA(x),
where dA is the surface measure and ω is a unit vector in R3, this is well
known to reduce the 3D Laplace operator ∆x to the second derivative ∂
2/∂s2
[30, 38, 39, 40, 53, 54], and thus the 3D wave equation to the string vibra-
tion problem. The measured data provide the boundary conditions for this
problem. The initial conditions in (1) mean evenness with respect to time,
and thus the standard d’Alambert formula leads to the immediate realiza-
tion that the measured data is just the 3D Radon transform of f(x). Thus,
the reconstruction boils down to the well known inversion formulas for the
Radon transform.
Another proposal ([20, 21], [49]-[52], [92, 93]) is to use line detectors
that provide line integrals of the pressure p(x, t). Such detectors can be
implemented optically, using either Fabry-Perot [20], or Mach-Zehnder [93]
interferometers.
Suppose that the object is surrounded by a surface that is rotation in-
variant with respect to the z-axis. It is suggested to place the line detectors
perpendicular to the z-axis and tangential to the surface. The same consid-
eration as above then shows that after the 2D Radon (or X-ray, which in
2D is the same) transform in each plane orthogonal to z-axis, the 3D wave
equation converts into the 2D one for the Radon data. The measurements
provide the boundary data. Thus, the reconstruction boils down to solving
a 2D problem similar to the one in the case of point detectors, and then
inverting the 2D Radon transform.
Due to the recent nature of these two projects, it appears to be too early
to judge which one will be superior in the end. For instance, it is not clear
beforehand, whether the approximation of infinite size (length, area) of the
linear or planar detectors works better than the zero dimension approxima-
tion for point detectors. Further developments will resolve these questions.
8.1.2 Direct imaging techniques
Some direct imaging techniques have been suggested, which might not require
mathematical reconstructions. See, for instance, [79] about an acoustic lens
system.
8.1.3 Using contrast agents
Contrast agents to improve TAT imaging have been developed (e.g., [24]).
8.1.4 Passive thermoacoustic imaging
The TAT model we have considered can be called “active thermoacoustic to-
mography,” due to the set-up when the practitioner creates the signal. There
has been some recent development of the “passive thermoacoustic tomogra-
phy,” where the thermoacoustic signal is used to image the temperature
sources present inside the body. One can find a survey of this area in [94].
8.2 Uniqueness
8.2.1 Sketch of the proof of Theorem 10
We provide here a brief outline of the rather technical proof of Theorem 10.
Suppose that f is compactly supported, not identically zero, and such
that RSf = 0. Our previous considerations show that one can assume that
S is an algebraic curve (not a straight line) that is contained in the set of
zeros of a non-trivial harmonic polynomial. Now one touches the boundary
of the support of f from outside by a circle centered on S. Then microlocal
analysis of the operator RS (which happens to be an analytic Fourier Integral
Operator, FIO [19, 42, 43, 44, 45, 65, 98]) shows that, due to the equality
RSf = 0, at the tangency point the vector co-normal to the sphere should not
belong to the analytic wave front of f (microlocal regularity of solutions of
RSf = 0). This, for instance, can be also extracted from the results of [105].
On the other hand, a theorem by Hörmander and Kashiwara [56, Theorem
8.5.6] shows that this vector must be in the analytic wave front set, since
f = 0 on one side of the sphere (a microlocal version of uniqueness of analytic
continuation). This way, one gets a contradiction. Unfortunately, the life is
not so easy, and the proof sketched above does not go through smoothly, due
to possible cancelation of wavefronts at different tangency points. Then one
has to involve the geometry of zeros of harmonic polynomials [37] to exclude
the possibility of such a cancelation.
Thus, the proof uses microlocal analysis and geometry of zeros of har-
monic polynomials. Both these tools have their limitations. For instance,
the microlocal approach (at least, in the form it is used in [7]) does not al-
low considerations of non-compactly supported functions. Thus, the validity
of the Theorem for arbitrarily fast decaying, but not compactly supported,
functions is still not established, albeit it most certainly holds. On the other
hand, the geometric part does not work that well in dimensions higher than
two. Development of new approaches is apparently needed in order to over-
come these hurdles. A much simpler PDE approach has emerged recently
[33] (see also [12] and the next Section), albeit its achievements have been
limited so far.
8.2.2 Some open problems concerning uniqueness
As it has already been mentioned, one can consider the practical problems
about uniqueness resolved. However, the mathematical understanding of the
uniqueness problem for the restricted spherical mean operators RS is still
unsatisfactory. Here are some questions that still await their resolution:
1. Describe uniqueness sets in dimensions larger than 2 (prove the Con-
jecture 4). Recent limited progress, as well as variations on this theme
can be found in [1]-[12].
2. Prove Theorem 3 without using microlocal and harmonic polynomial
tools.
3. Prove Theorem 3 on uniqueness sets S under the condition of suffi-
ciently fast decay (rather than compactness of support) of the function.
Very little is known for the case of functions without compact support.
The main known result is of [3], which describes for which values of
1 ≤ p ≤ ∞ the result of Corollary 2 still holds:
Theorem 12 [3] Let S be the boundary of a bounded domain in Rd and
f ∈ Lp(Rd) such that RSf ≡ 0. If p ≤ 2d/(d−1), then f ≡ 0 (and thus
S is injectivity set for this space). This fails for any p > 2d/(d− 1).
8.3 Inversion
Albeit closed form (backprojection type) inversion formulas are available now
for the cases of S being a plane (and object on one side from it), cylinder,
and a sphere, there is still some mystery surrounding this issue.
1. Can one write a backprojection type inversion formula in the case of
the constant sound speed for a closed surface S which is not a sphere?
We suspect that the answer to this question is negative (see also related
discussion in [15, 27]).
2. The inversion formulas for S being a sphere assume that the object to
be imaged is inside S. One can check on simplest examples that if the
support of function f(x) reaches outside S, the inversion formulas do
not reconstruct the function correctly even inside of S. See [5] for a
discussion.
3. The I. Gelfand’s school of integral geometry has developed a marvelous
machinery of the so called κ operator, which provides a general ap-
proach to inversion and range descriptions for transforms of Radon type
[38, 39]. In particular, it has been applied to the case of integration of
various collections (“complexes”) of spheres in [39, 41]. This consider-
ation seems to suggest that one should not expect explicit closed form
inversion formulas for RS when S is a sphere. We, however, know that
such formulas have been discovered recently [33, 69]. This apparent
controversy has not been resolved.
4. Can one derive any more explicit analytic formulas from (28)?
5. Can the series expansion formulas of Theorem 8 be efficiently imple-
mented?
One can also mention that in some works [15, 23] it is suggested to use in
the TAT problem not only the values of the pressure measured by transducers
on the observation surface S, but its normal derivative to S as well. If one
knows both, then taking Fourier transform in the time variable and using the
whole space Green’s function for the Helmholtz equation leads immediately
to a reconstruction formula for the solution (which seems to be much simpler
than what is proposed in [23]). The problem is that this normal derivative is
not measured by TAT devices. Under some circumstances (e.g., when there
are no sources of ultrasound outside S), one can prove the theoretical pos-
sibility of recovering the missing normal derivative. This, however, does not
seem to us to be a plausible procedure. In rare cases (planar, cylindrical,
or spherical surface S), when involvement of the normal derivative can be
eliminated (e.g., [15, 27]), this might lead to feasible inversion algorithms,
but in these cases, as explained before in this text, explicit and nicely im-
plementable analytic inversion formulas are available. So, jury is still out on
this issue as well.
8.4 Stability
Stability of inversion when S is a sphere surrounding the support of f(x) is
the same as for the standard Radon transform, as the results of [91] and sec-
ond statement of Theorem 11 show. However, if the support reaches outside,
albeit Corollary 2 still guarantees uniqueness of reconstruction, stability (at
least for the parts outside S) is gone. Indeed, Theorem 9 shows that some
parts of singularities of f outside S will not be stably “visible.”
8.5 Range
As Theorem 9 states, the range conditions 2 and 3 of Theorem 10 are neces-
sary also for non-spherical closed surfaces S and for functions with support
outside S. They, however, are not expected to be sufficient, since Theorem 9
indicates that one might expect non-closed ranges in some cases. The same
applies for non-constant sound speed case.
Acknowledgments
The work of the first author was partially supported by the NSF DMS grants
0604778 and 0648786. The second author was partially supported by the
DOE grant DE-FG02-03ER25577 and NSF DMS grant 0312292. Part of
this work was completed when the first author was at the Isaac Newton
Institute for Mathematical Sciences. The authors express their gratitude to
the NSF, DOE and INI for this support. The authors thank M. Agranovsky
for extremely useful discussions, M. Anastasio, G. Beylkin and M. Klibanov
for providing preprints and references, and the reviewers for very helpful
remarks.
References
[1] Agranovsky, M. 1997 Radon transform on polynomial level sets and
related problems. Israel Math. Conf. Proc., 11, 1-21.
[2] Agranovsky, M. 2000 On a problem of injectivity for the Radon trans-
form on a paraboloid. Analysis, geometry, number theory: the mathe-
matics of Leon Ehrenpreis (Philadelphia, PA, 1998). Contemp. Math.
251, AMS, Provodence, RI, 1–14.
[3] Agranovsky, M., Berenstein, C., & Kuchment, P. 1996 Approximation
by spherical waves in Lp-spaces, J. Geom. Anal. 6, no. 3, 365–383.
[4] Agranovsky, M. & Kuchment, P. 2007 Uniqueness of reconstruction and
an inversion procedure for thermoacoustic and photoacoustic tomogra-
phy with variable sound speed. Inverse Problems 23, 2089–2102.
[5] Agranovsky, M., Kuchment, P. & Kunyansky, L. 2007 On reconstruc-
tion formulas and algorithms for the thermoacoustic and photoacoustic
tomography. Submitted.
[6] Agranovsky, M., Kuchment, P., & Quinto, E. T. 2007 Range descriptions
for the spherical mean Radon transform. J. Funct. Anal. 248, 344–386.
[7] Agranovsky, M. & Quinto, E. T. 1996 Injectivity sets for the Radon
transform over circles and complete systems of radial functions. Journal
of Functional Analysis, 139, 383–414.
[8] Agranovsky, M. & Quinto, E. T. 2001 Geometry of stationary sets for
the wave equation in Rn:the case of finitely suported initial data. Duke
Math. J., 107, no. 1, 57–84.
[9] Agranovsky, M. & Quinto, E. T. 2003 Stationary sets for the wave equa-
tion in crystallographic domains. Trans. AMS, 355, no. 6, 2439–2451.
[10] Agranovsky, M. & Quinto, E. T. 2006 Remarks on stationary sets for the
wave equation. Integral Geometry and Tomography, Contemp. Math.
405, 1–11.
[11] Agranovsky, M., Volchkov, V. V. & Zalcman, L. 1999 Conical uniqueness
sets for the spherical Radon transform. Bull. London Math. Soc., 31, no.
4, 363–372.
[12] Ambartsoumian, G. & Kuchment, P. 2005 On the injectivity of the
circular Radon transform. Inverse Problems 21, 473–485.
[13] Ambartsoumian, G. & Kuchment, P. 2006 A range description for the
planar circular Radon transform. SIAM J. Math. Anal. 38, no. 2, 681–
[14] Ambartsoumian, G. & Patch, S. 2007 Thermoacoustic tomography: nu-
merical results. Proceedings of SPIE 6437, Photons Plus Ultrasound:
Imaging and Sensing 2007: The Eighth Conference on Biomedical Ther-
moacoustics, Optoacoustics, and Acousto-optics, Alexander A. Oraevsky,
Lihong V. Wang, Editors, 64371B.
[15] Anastasio, M. A., Zhang, J., Modgil, D., and La Rivière, P. J. 2007
Application of inverse source concepts to photoacoustic tomography,
preprint.
[16] Andersson, L.-E.. 1988 On the determination of a function from spherical
averages. SIAM J. Math. Anal. 19 no. 1, 214–232.
[17] Asgeirsson, L. 1937 Über eine Mittelwerteigenschaft von Lösungen ho-
mogener linearer partieller Differentialgleichungen zweiter Ordnung mit
konstanten Koeffizienten. Ann. Math., 113, 321–346.
[18] Bangerth, W., Georgieva-Hristova, Y. & Kuchment, P. 2007 On recon-
struction in thermoacoustic tomography with variable speed, in prepa-
ration.
[19] Beylkin, G. 1984 The inversion problem and applications of the gener-
alized Radon transform. Comm. Pure Appl. Math. 37, 579–599.
[20] Burgholzer, P., Hofer, C., Paltauf, G., Haltmeier, M., & Scherzer, O.
2005 Thermoacoustic tomography with integrating area and line detec-
tors. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency
Control 52(9), 1577–1583.
[21] Burgholzer, P.,, Hofer, C., Matt, G. J., Paltauf, G., Haltmeier, M.
& Scherzer, O. 2006 Thermoacoustic tomography using a fiber-based
Fabry-Perot interferometer as an integrating line detector. Proc. SPIE
6086, 434–442.
[22] Burgholzer, P., Matt, G., Haltmeier, M. & Patlauf, G. 2007 Exact and
approximate imaging methods for photoacoustic tomography using an
arbitrary detection surface. Phys. Rev. E 75, 046706.
[23] Clason, C. and Klibanov, M. 2007 Quasireversibility method in ther-
moacoustic tomography in heterogeneous medium, preprint.
[24] Copland, J. A. et al. 2004 Bioconjugated gold nanoparticles as a molec-
ular based contrast agent: implications for imaging of deep tumors using
optoacoustic tomography. Molecular Imaging and Biology 6, no. 5, 341–
[25] Courant, R. & Hilbert, D. 1962 Methods of Mathematical Physics, Vol-
ume II Partial Differential Equations, Interscience, New York.
[26] Denisjuk, A. 1999 Integral geometry on the family of semi-spheres. Fract.
Calc. Appl. Anal. 2, no. 1, 31–46.
[27] Devaney, A. J. and Beylkin, G. 1984 Diffraction tomography using arbi-
trary transmitter and receiver surfaces. Ultrasonic Imaging 6, 181-193.
1984.
[28] Diebold, G. J., Sun, T. & Khan, M. I. 1991 Photoacoustic monopole
radiation in one, two, and three dimensions. Phys. Rev. Lett. 67, no.
24, 3384–3387.
[29] Egorov, Yu. V. & Shubin, M. A. 1992 Partial Differential Equations I.
Encyclopaedia of Mathematical Sciences, (Springer Verlag), 30, 1–259
[30] Ehrenpreis, L. 2003 The Universality of the Radon Transform, Oxford
Univ. Press.
[31] Fawcett, J. A. 1985 Inversion of n-dimensional spherical averages. SIAM
J. Appl. Math. 45, no. 2, 336–341.
[32] Finch, D., Haltmeier, M. & Rakesh 2007 Inversion of spherical
means and the wave equation in even dimensions. Preprint arXiv
math.AP/0701426.
[33] Finch, D., Patch, S. & Rakesh 2004 Determining a function from its
mean values over a family of spheres. SIAM J. Math. Anal. 35, no. 5,
1213–1240.
[34] Finch, D. & Rakesh 2006 The range of the spherical mean value operator
for functions supported in a ball. Inverse Problems 22, 923-938.
[35] Finch, D. & Rakesh 2007 Recovering a function from its spherical mean
values in two and three dimensions. Preprint.
http://arxiv.org/abs/math/0701426
[36] Finch, D. & Rakesh 2007 The spherical mean value operator with centers
on a sphere. Preprint. To appear in Inverse Problems.
[37] Flatto, L., Newman, D. J. & Shapiro, H. S. 1966 The level curves of
harmonic functions. Trans. Amer. Math. Soc. 123, 425–436.
[38] Gelfand, I., Gindikin, S. & Graev M. 1980 Integral geometry in affine
and projective spaces. J. Sov. Math. 18, 39–167.
[39] Gelfand, I., Gindikin, S. & Graev M. 2003 Selected Topics in Integral
Geometry. Transl. Math. Monogr. v. 220, Amer. Math. Soc., Providence
[40] Gelfand, I., Graev M. & Vilenknin, N. 1965 Generalized Functions, v.
5: Integral Geometry and Representation Theory, Acad. Press.
[41] Gindikin, S. 1995 Integral geometry on real quadrics, in Lie groups and
Lie algebras: E. B. Dynkin’s Seminar, 23–31, Amer. Math. Soc. Transl.
Ser. 2, 169, Amer. Math. Soc., Providence, RI.
[42] Greenleaf, A. & Uhlmann, G. 1990 Microlocal techniques in integral
geometry. Contemporary Math. 113, 149–155.
[43] Guillemin, V. 1975 Fourier integral operators from the Radon transform
point of view. Proc. Symposia in Pure Math., 27, 297–300.
[44] Guillemin, V. 1985 On some results of Gelfand in integral geometry.
Proc. Symposia in Pure Math., 43, 149–155.
[45] Guillemin, V. & Sternberg S. 1977 Geometric Asymptotics. Amer. Math.
Soc., Providence, RI.
[46] Gusev, V. E. & Karabutov, A. A. 1993 Laser Optoacoustics. American
Inst. of Physics, NY.
[47] Haltmeier, M., Burgholzer, P., Paltauf, G. & Scherzer, O. 2004 Ther-
moacoustic computed tomography with large planar receivers. Inverse
Problems 20, 1663–1673.
[48] Haltmeier, M., Schuster, T. & O. Scherzer. 2005 Filtered backprojec-
tion for thermoacoustic computed tomography in spherical geometry.
Mathematical Methods in the Applied Sciences, 28, 1919–1937.
[49] Haltmeier, M., Paltauf, G., Burgholzer, P. & Scherzer, O. 2005 Ther-
moacoustic Tomography with integrating line detectors. Proc. SPIE
5864:586402-8.
[50] Haltmeier, M., Burgholzer, P., Hofer, C., Paltauf, G., Nuster, R. &
Scherzer, O. 2005 Thermoacoustic tomography using integrating line
detectors. Ultrasonics Symposium 1, 166–169.
[51] Haltmeier, M., Scherzer, O., Burgholzer, P. & Paltauf, G. 2005 Ther-
moacoustic Computed Tomography with large planar receivers. ECMI
Newsletter 37, pp. 31-34. http://www.it.lut.fi/mat/EcmiNL/ecmi37/
[52] Haltmeier, M., & Fidler, T. Mathematical Challenges Aris-
ing in Thermoacoustic Tomography with Line Detectors, preprint
arXiv:math.AP/0610155.
[53] Helgason, S. 1980 The Radon Transform, Birkhäuser, Basel.
[54] Helgason, S. 2000 Groups and Geometric Analysis. Amer. Math. Soc.,
Providence, R.I.
[55] Herman, G.(Ed.) 1979 Image Reconstruction from Projections . Topics
in Applied Physics, v. 32, Springer Verlag, Berlin, New York.
[56] Hörmander, L. 1983 The Analysis of Linear Partial Differential Opera-
tors, vol. 1, Springer-Verlag, New York.
[57] Jin, X. & Wang, L. V. 2006 Thermoacoustic tomography with correction
for acoustic speed variations. Physics in Medicine and Biology 51, 6437–
6448.
[58] John, F. 1971 Plane Waves and Spherical Means Applied to Partial
Differential Equations, Dover.
[59] Kak, A. C. & Slaney, M. 2001 Principles of Computerized Tomographic
Imaging. SIAM, Philadelphia.
[60] Köstli K. P., Frenz, M., Bebie, H. & Weber H. P. 2001 Temporal back-
ward projection of optoacoustic pressure transients using Fourier trans-
form methods. Phys. Med. Biol. 46, 1863–1872
http://www.it.lut.fi/mat/EcmiNL/ecmi37/
http://arxiv.org/abs/math/0610155
[61] Kruger, R. A., Kiser, W. L., Reinecke, D. R. & Kruger, G. A. 2003 Ther-
moacoustic computed tomography using a conventional linear trans-
ducer array. Med. Phys. 30, no 5, 856–860.
[62] Kruger, R. A., Liu. P., Fang, Y. R. & Appledorn, C. R. 1995 Photoa-
coustic ultrasound (PAUS)reconstruction tomography. Med. Phys. 22,
1605–1609.
[63] Kuchment, P. 1993, unpublished.
[64] Kuchment, P. 2006 Generalized Transforms of Radon Type and Their
Applications. in [85], pp. 67–91.
[65] Kuchment, P., Lancaster, K. & Mogilevskaya, L. 1995 On local tomog-
raphy. Inverse Problems, 11, 571–589.
[66] Kuchment, P. & Lvin S. 1990 Paley-Wiener theorem for the exponential
Radon transform. Acta Applicandae Mathematicae, no.18, 251–260.
[67] Kuchment, P. & Lvin S. 1991 The Range of the Exponential Radon
Transform. Soviet Math Dokl, 42, no.1, 183–184.
[68] Kuchment, P. & Quinto, E. T. 2003. Some problems of integral geometry
arising in tomography. Chapter XI in [30].
[69] Kunyansky, L. 2007 Explicit inversion formulae for the spherical mean
Radon transform. Inverse problems 23 (2007), 737-783.
[70] Kunyansky, L. 2007 A series solution and a fast algorithm for the
inversion of the spherical mean Radon transform. Preprint arXiv
math.AP/0701236.
[71] Lin, V. & Pinkus, A. 1993 Fundamentality of ridge functions. J. Approx.
Theory, 75, 295–311.
[72] Lin, V. & Pinkus, A 1994 Approximation of multivariate functions. In
Advances in Computational Mathematics, H. P. Dikshit & C. A. Mic-
chelli, Eds., World Sci. Publ., 1–9.
[73] Louis, A. K. & Quinto, E. T. 2000 Local tomographic methods in Sonar.
In Surveys on solution methods for inverse problems, Springer, Vienna,
147–154.
http://arxiv.org/abs/math/0701236
[74] Lvin S. 1994 Data correction and restoration in emission tomography,
in E.T. Quinto, M. Cheney, and P. Kuchment (Editors), Tomography,
Impedance Imaging, and Integral Geometry, Lectures in Appl. Math.,
vol. 30, AMS, Providence, RI, 149–155.
[75] Mathematics and Physics of Emerging Biomedical Imaging,
The National Academies Press 1996. Available online at
http://www.nap.edu/catalog.php?record id=5066#toc.
[76] Natterer, F. 1986 The mathematics of computerized tomography, Wiley,
New York.
[77] Natterer, F. & Wübbeling, F. 2001 Mathematical Methods in Image Re-
construction, Monographs on Mathematical Modeling and Computation
5, SIAM, Philadelphia, PA.
[78] Nessibi, M. M., Rachdi, L. T. & Trimeche, K. 1995 Ranges and inversion
formulas for spherical mean operator and its dual. J. Math. Anal. Appl.
196, no. 3, 861–884.
[79] Niederhauser, J. J., Jaeger, M., Lemor, R., Weber, P. & Frenz, M.
2005 Combined ultrasound and optoacoustic system for real-time high-
contrast vascular imaging in vivo. IEEE Transactions on medical Imag-
ing 24, 436–440.
[80] Nilsson, S. 1997 Application of fast backprojection techniques for some
inverse problems of integral geometry. Linkoeping studies in science and
technology, Dissertation 499, Dept. of Mathematics, Linkoeping univer-
sity, Linkoeping, Sweden.
[81] Nolan, C. J. & Cheney, M. 2002 Synthetic aperture inversion. Inverse
Problems 18, 221–235.
[82] Norton, S. J. 1980 Reconstruction of a two-dimensional reflecting
medium over a circular domain: exact solution. J. Acoust. Soc. Am.
67, 1266–1273.
[83] Norton, S. J. & Linzer, M. 1981 Ultrasonic reflectivity imaging in three
dimensions: exact inverse scattering solutions for plane, cylindrical, and
spherical apertures. IEEE Transactions on Biomedical Engineering, 28,
200–202.
http://www.nap.edu/catalog.php?record_id=5066#toc
[84] Novikov, R. 2002 On the range characterization for the two-dimensional
attenuated X-ray transform. Inverse Problems 18, 677–700.
[85] Olafsson, G. & Quinto, E. T. (Editors), 2006 The Radon Transform, In-
verse Problems, and Tomography. American Mathematical Society Short
Course January 3–4, 2005, Atlanta, Georgia, Proc. Symp. Appl. Math.,
v. 63, AMS, RI.
[86] Oraevsky, A. A., Esenaliev, R. O., Jacques, S. L. Tittel, F. K. 1996
Laser optoacoustic tomography for medical diagnostics principles. Proc.
SPIE 2676, 22.
[87] Oraevsky, A. A. & Karabutov, A. A. 2002 In Handbook of Optical
Biomedical Diagonstics, edited by V. V. Tuchin, SPIE, Bellingham, WA,
Chap. 10.
[88] Oraevsky A. A.& A. A. Karabutov, A. A., 2003 Optoacoustic Tomogra-
phy, Ch. 34 In Biomedical Photonics Handbook, edited by T. Vo-Dinh,
CRC, Boca Raton, FL, Chap. 34, 34-1 – 34-34.
[89] Palamodov, V. P. 2000 Reconstruction from limited data of arc means.
J. Fourier Anal. Appl. 6, no. 1, 25–42.
[90] Palamodov, V. P. 2004 Reconstructive Integral Geometry, Birkhäuser,
Basel.
[91] Palamodov, V. 2006 Remarks on the general Funk transform. Preprint,
Tel Aviv University, August.
[92] Paltauf, G., Burgholzer, P., Haltmeier, M. & O. Scherzer 2005 Ther-
moacoustic Tomography using optical Line detection. Proc. SPIE 5864,
7–14.
[93] Paltauf, G., R. Nuster, Haltmeier, M. & Burgholzer, P. 2007(?) Ther-
moacoustic Computed Tomography using a Mach-Zehnder interferome-
ter as acoustic line detector. Submitted
[94] Passechnik, V. I., Anosov, A. A. & Bograchev, K. M. 2000 Fundamentals
and prospects of passive thermoacoustic tomography. Critical reviews in
Biomed. Eng. 28, no. 3&4, 603–640.
[95] Patch, S. K. 2004 Thermoacoustic tomography - consistency conditions
and the partial scan problem. Phys. Med. Biol. 49, 1–11.
[96] Popov D. A. & Sushko, D. V. 2002 A parametrix for the problem of
optical-acoustic tomography. Dokl. Math. 65, no. 1, 19–21.
[97] Popov D. A. & Sushko, D. V. 2004 Image restoration in optical-acoustic
tomography. Problems of Information Transmission 40, no. 3, 254–278.
[98] Quinto, E. T. 1980 The dependence of the generalized Radon transform
on defining measures. Trans. Amer. Math. Soc. 257, 331–346.
[99] Quinto, E. T. 1993 Singularities of the X-ray transform and limited data
tomography in R2 and R3. SIAM J. Math. Anal. 24, 1215–1225.
[100] Quinto, E. T. 2006 An introduction to X-ray tomography and Radon
transforms. In [85], 1–23.
[101] Ramm, A. G. 1985 Inversion of the backscattering data and a problem
of integral geometry. Phys. Lett. A 113, no. 4, 172–176.
[102] Ramm, A. G. 2002 Injectivity of the spherical means operator. C. R.
Math. Acad. Sci. Paris 335, no. 12, 1033–1038.
[103] Romanov, V. G. 1967 Reconstructing functions from integrals over a
family of curves. Sib. Mat. Zh. 7, 1206–1208.
[104] Schuster, T. & Quinto, E. T. 2005 On a regularization scheme for linear
operators in distribution spaces with an application to the spherical
Radon transform. SIAM J. Appl. Math. 65 (4), 1369–1387.
[105] Stefanov, P. & Uhlmann, G. Integral geometry of tensor fields on a class
of non-simple Riemannian manifolds. Preprint arXiv:math/0601178.
[106] Strichartz, Robert S. 2003 A Guide to Distribution Theory and Fourier
Transforms. World. Sci.
[107] Tam, A. C. 1986 Applications of photoacoustic sensing techniques. Rev.
Mod. Phys. 58, no. 2, 381–431.
[108] Tataru, D. 1995 Unique continuation for solutions to PDEs; between
Hörmander’s theorem and Holmgren’s theorem. Comm. PDE 20, 814–
http://arxiv.org/abs/math/0601178
[109] Tuchin, V. V. (Editor) 2002 Handbook of Optical Biomedical Diagnos-
tics. SPIE, Bellingham, WA.
[110] Vainberg, B. 1975 The short-wave asymptotic behavior of the solutions
of stationary problems, and the asymptotic behavior as t → ∞ of the
solutions of nonstationary problems. Russian Math. Surveys, 30, no. 2,
1–58.
[111] Vainberg, B. 1982. Asymptotics methods in the Equations of Mathe-
matical Physics. (Gordon & Breach)
[112] Vo-Dinh, T. (Editor) 2003 Biomedical Photonics Handbook. CRC, Boca
Raton, FL.
[113] Wang, L. V. &Wu, H. 2007 Biomedical Optics. Principles and Imaging.
Wiley-Interscience.
[114] Wang, X., Pang, Y., Ku, G., Xie, X., Stoica, G. & Wang, L. 2003
Noninvasive laser-induced photoacoustic tomography for structural and
functional in vivo imaging of the brain. Nature Biotechnology, 21, no.
7, 803–806.
[115] Xu, M. & Wang, L.-H. V. 2002 Time-domain reconstruction for ther-
moacoustic tomography in a spherical geometry. IEEE Trans. Med.
Imag. 21, 814–822.
[116] Xu, M. & Wang, L.-H. V. 2005 Universal back-projection algorithm for
photoacoustic computed tomography. Phys. Rev. E 71, 016706.
[117] Xu, M. & Wang, L.-H. V. 2006 Photoacoustic imaging in biomedicine.
Review of Scientific Instruments 77, 041101-01 – 041101-22.
[118] Xu, Y., Feng, D. & Wang, L.-H. V. 2002 Exact frequency-domain re-
construction for thermoacoustic tomography: I. Planar geometry. IEEE
Trans. Med. Imag. 21, 823–828.
[119] Xu, Y., Xu, M. & Wang, L.-H. V. 2002 Exact frequency-domain re-
construction for thermoacoustic tomography: II. Cylindrical geometry.
IEEE Trans. Med. Imag. 21, 829–833.
[120] Xu, Y., Wang, L., Ambartsoumian, G. & Kuchment, P. 2004 Recon-
structions in limited view thermoacoustic tomography. Medical Physics,
31(4), 724–733.
[121] Zobin, N. 1993. Unpublished.
	Introduction
	Thermoacoustic tomography
	Mathematical model of TAT: wave equation and the spherical mean transform
	The wave equation model
	Spherical mean model
	Main mathematical problems of TAT
	Uniqueness of reconstruction
	Constant speed case
	Non-uniqueness sets in R2.
	Higher dimensions
	Relations to other areas of analysis
	Uniqueness in the case of a variable sound speed
	Reconstruction: formulas and examples
	Constant sound speed
	Inversion formulas
	Fourier expansion methods
	Filtered backprojection methods
	Series solutions for arbitrary geometries
	Time reversal (backpropagation) methods
	Examples of reconstructions and additional remarks about the inversion formulas
	Reconstruction in the variable speed case
	``Analytic'' inversions
	Backpropagation
	Eigenfunction expansions
	Partial data. ``Visible'' and ``invisible'' singularities
	Range conditions
	Concluding remarks
	Variations of the TAT procedure
	Planar and linear transducers
	Direct imaging techniques
	Using contrast agents
	Passive thermoacoustic imaging
	Uniqueness
	Sketch of the proof of Theorem 10
	Some open problems concerning uniqueness
	Inversion
	Stability
	Range
