Microlocally accurate solution of the inverse source problem of thermoacoustic tomography

We consider the inverse source problem for the free-space wave equation, with measurements made on an open bounded surface. Such problems arise in the thermo- and photoacoustic tomography (TAT/PAT) [1–3] and in several other coupled-physics modalities (see, for example, [4–6]). The forward problem can be modeled by the Cauchy problem for the standard wave equation in ${\mathbb{R}}^{d}$:

$$\begin{equation}{c}^{2}\left(x\right){\Delta}u\left(t,x\right)={u}_{tt}\left(t,x\right),\quad t{ >}0,\enspace x\in {\mathbb{R}}^{d},\end{equation} \tag{ 1 }$$

with the initial condition

$$\begin{equation}u\left(0,x\right)=f\left(x\right),\quad {u}_{t}\left(0,x\right)=0,\end{equation} \tag{ 2 }$$

where, in the case of TAT/PAT, u(t, x) represents the acoustic pressure in the tissues, c(x) is the speed of sound, and the initial pressure f(x) results from the thermoelastic effect. Function f is assumed to be compactly supported within a bounded region Ω₀ that itself is a subset of a larger bounded region Ω ⊃ Ω₀. Ideally, one would take the measurements g(t, z) defined as

$$\begin{equation}g\left(t,z\right)\equiv u\left(t,z\right),\quad z\in \partial {\Omega},\enspace t\in \left(0,\infty \right),\end{equation} \tag{ 3 }$$

where ∂Ω is the boundary Ω. However, here we consider the practically important case of measurements restricted to a subset Γ of ∂Ω; we will denote this data g^reduced(t, z):

$$\begin{equation}{g}^{\text{reduced}}\left(t,z\right)\equiv \begin{cases}g\left(t,z\right)=u\left(t,z\right),\hfill & \quad z\in {\Gamma}\subset \partial {\Omega},\hfill \\ 0,\hfill & \quad z\in \partial {\Omega}{\backslash}{\Gamma},\hfill \end{cases}\quad t\in \left(0,\infty \right).\end{equation} \tag{ 4 }$$

The two inverse source problems of interest to us consist of finding f(x) from known c(x) and either g(t, z) or g^reduced(t, z).

The inverse source problem with complete data g(t, z) is very well understood by now. In particular, the so-called time-reversal procedure (applicable under the non-trapping condition) is based on solving backwards in time the mixed initial/boundary value problem for the wave equation in the domain (0, ∞) × Ω, with uniform initial conditions at the limit t → ∞ and the boundary values g(t, z) [7–10]. Many other techniques have been developed for the inverse source problem with complete data, including a variety of explicit inversion formulas applicable if c(x) is constant (we will not attempt to overview these results here).

However, in practical situations one can measure data only on a subset Γ of ∂Ω, and has to work with g^reduced. For such data very few reconstruction technique are presently known, most of them iterative in nature. In the present paper under the assumption of the constant speed of sound, we develop an explicit non-iterative image reconstruction procedure that recovers f(x) up to an infinitely smooth error term. Our numerical experiments show that this error is also quite small from the practical standpoint, although we are able to prove only the smoothness of it.

We make the following assumptions about the data acquisition scheme. The region Ω containing Ω₀ is convex and has a smooth boundary ∂Ω. Without further loss of generality we assume that Ω is contained inside the unit sphere ${\mathbb{S}}^{d-1}$ centered at the origin (see figure 1), and that the constant speed of sound c is equal to 1. The support Ω₀ of f lies within Ω and above the plane x^d = 0. The measurement surface Γ is a subset of ∂Ω containing all points lying above the plain x^d = −2δ, where δ is a fixed positive number:

$$\begin{equation*}{\Gamma}=\left\{x\enspace :\enspace x\in \partial {\Omega}\;\text{and}\;{x}^{d}{\geqslant}-2\delta \right\}.\end{equation*}$$

Our goal is to solve

Zoom In Zoom Out Reset image size

Problem 1 (Inverse source problem with reduced data).  Find f(x) given data g^reduced(t, z) defined by equations (1), (2) and (4) .

However, it will be helpful for us to start the discussion with

Problem 2 (Inverse source problem with complete data).  Find f(x) given data g(t, z) defined by equations (1)–(3).

The solvability and stability of the inverse problem with reduced data have been thoroughly studied under various assumptions on the speed of sound c(x) (see [7, 8, 10]). Importantly, the stability of the reconstruction is guaranteed under the so-called visibility condition [11], which for the case of constant speed of sound can be formulated as follows

Visibility condition. For every point x ∈ Ω₀, every straight line passing through x intersects Γ at least once.

Our data acquisition scheme satisfies this condition.

The most straightforward constructive approach to the problem with the reduced data is to apply the standard time-reversal technique to g^reduced to obtain an approximation to f. It is shown in [10, 12] that, under the visibility condition, such approximation is equivalent to applying a pseudo-differential operator (ψDO) of order zero to f. A Neumann series approach was used in [13] to invert such an operator; however, convergence of this approach is not well understood theoretically. Instead of the time reversal one can apply the operator adjoint to the operator defining the direct map of to g^reduced [14]. This still leads to expensive iterations, although they are guaranteed to converge eventually.

Non-iterative approaches to this problem include [15–18]. The numerical algorithm [17] produces good results in 2D; however, it requires further theoretical investigation and it might become expensive for 3D computations. The method presented in [18] yields a theoretically exact reconstruction of f, but only under a certain sub-optimal geometrical condition: the support of f must be significantly smaller than the region that otherwise satisfies the visibility condition. Finally, in [15, 16] the source term f is reconstructed up to an error term given by a ψDO of order −1 applied to f.

In the preset paper we first develop a novel explicit procedure that yields theoretically exact reconstruction of the Radon projections of f if the complete data g are given. When applied to reduced data, this algorithm produces the Radon projections accurate up to an infinitely smooth term. This allows for an explicit reconstruction of f up to a smooth error. In sections 3 and 4 we formulate and prove the main theorems defining our method. An asymptotically fast algorithm for the particular case of the circular arc acquisition curve in 2D is presented in section 5. In the last section we demonstrate the performance of this approach in numerical simulations.

2.1. Theoretically exact reconstruction from complete data

Let us first introduce an explicit reconstruction procedure that is theoretically exact in the case of complete data. In other words, in this section we are solving problem 2.

2.1.1. Exterior problem

Unlike the classical approach that consists in solving backwards in time the wave equation in a domain (0, ∞) × Ω, we first solve an exterior Dirichlet problem for the wave equation in the domain ${Q}^{\text{ext}}\equiv \left(0,\infty \right){\times}\left({\mathbb{R}}^{d}{\backslash}\overline{{\Omega}}\right)$. In order words, we find solution u(t, x) of the initial/boundary value problem

$$\begin{equation}\begin{aligned}\hfill {u}_{tt}\left(t,x\right)-{\Delta}u\left(t,x\right)& =0,\quad \left(t,x\right)\in {Q}^{\text{ext}},\hfill \\ \hfill u\left(0,x\right)& =0,\quad {u}_{t}\left(0,x\right)=0,\enspace x\in {\mathbb{R}}^{d}{\backslash}{\Omega},\hfill \\ \hfill u\left(t,z\right)& =g\left(t,z\right),\quad z\in \partial {\Omega},\enspace t\in \left(0,\infty \right).\hfill \end{aligned}\end{equation} \tag{ 5 }$$

Notice that the solution of this problem coincides on Q^ext with the solution of the problem (1) and (2), which allows us to use the same notation u(t, x) for both functions.

2.1.2. Radon transform

The Radon transform is the main tool we use. To make this paper self-contained we collect in this section some basic facts about it. Consider a continuous function q(x) compactly supported within some bounded region ${{\Omega}}^{\text{large}}\subset {\mathbb{R}}^{d}$ (within the ball B^d(0, r) of radius r centered at the origin) and extended by zero to the rest of ${\mathbb{R}}^{d}$. For a direction $\omega \in {\mathbb{S}}^{d-1}$, and $p\in \mathbb{R}$, the values h(ω, p) of the Radon transform $\mathcal{R}q\left(\omega ,p\right)$ are given by the formula

$$\begin{equation}h\left(\omega ,p\right)=\left[\mathcal{R}q\right]\left(\omega ,p\right)\equiv \underset{{\mathbb{R}}^{d}}{\int }q\left(x\right)\delta \left(x\cdot \omega -p\right)\mathrm{d}x,\quad p\in \mathbb{R},\end{equation} \tag{ 6 }$$

where δ(⋅) is the Dirac's delta function ([19]). The Radon transform is thus defined on the cylinder $Z={\mathbb{S}}^{d-1}{\times}\mathbb{R}$. The obvious but important property of $\mathcal{R}$ is symmetry

$$\begin{equation}\mathcal{R}q\left(\omega ,p\right)=\mathcal{R}q{\left(-\omega ,-p\right)}^{2}.\end{equation} \tag{ 7 }$$

Let us introduce Sobolev norms on ${\mathbb{R}}^{d}$ and Z:

$$\begin{align*}\hfill {\Vert}q{{\Vert}}_{{H}^{s}\left({\mathbb{R}}^{d}\right)}^{2}& =\underset{{\mathbb{R}}^{d}}{\int }{\left[1+\vert \xi {\vert }^{2}\right]}^{2}\vert \hat{q}\left(\xi \right){\vert }^{2}\mathrm{d}\xi ,\hfill \\ \hfill {\Vert}h{{\Vert}}_{{H}^{s}\left(Z\right)}^{2}& =\underset{{\mathbb{S}}^{d-1}}{\int }\underset{\mathbb{R}}{\int }{\left[1+\vert \sigma {\vert }^{2}\right]}^{2}\vert \hat{h}\left(\omega ,\sigma \right){\vert }^{2}\mathrm{d}\sigma \enspace \mathrm{d}\omega ,\hfill \end{align*}$$

where $\hat{q}\left(\xi \right)$ is the Fourier transform of q(x) and (ω, σ) is the one-dimensional Fourier transform of h(ω, p) in the second variable. Then there exist constants c(s, d) and C(s, d) (see [19], chapter 2, sections 5) such that

$$\begin{equation*}c\left(s,d\right){\Vert}q{{\Vert}}_{{H}^{s}\left({\mathbb{R}}^{d}\right)}{\leqslant}{\Vert}\mathcal{R}q{{\Vert}}_{{H}^{s+\left(d-1\right)/2}\left(Z\right)}{\leqslant}C\left(s,d\right){\Vert}q{{\Vert}}_{{H}^{s}\left({\mathbb{R}}^{d}\right)}^{2}.\end{equation*}$$

Using these estimates, one extends the Radon transform to compactly supported distributions on ${\mathbb{R}}^{d}$.

2.1.3. Reconstruction procedure

We will denote by F(ω, p) the Radon transform of f

$$\begin{equation}F\left(\omega ,p\right)\equiv \left[\mathcal{R}f\right]\left(\omega ,p\right).\end{equation} \tag{ 8 }$$

Recall that u(t, x) is the solution of (1) and (2) with the initial condition f(x). Let us define the Radon projections $\left[\mathcal{R}u\right]\left(t,\omega ,p\right)$ of u(t, x) for each fixed t ⩾ 0 by (6). It is well known that, for a fixed ω, such projections solve the 1D wave equation:

$$\begin{equation}\frac{{\partial }^{2}}{\partial {t}^{2}}\left[\mathcal{R}u\right]\left(t,\omega ,p\right)=\frac{{\partial }^{2}}{\partial {p}^{2}}\left[\mathcal{R}u\right]\left(t,\omega ,p\right).\end{equation} \tag{ 9 }$$

The above equation is understood either in the classical sense or in the sense of distributions, depending on the smoothness of u. Further, by taking into account (2) one can confirm that

$$\begin{equation}\left[\mathcal{R}u\right]\left(0,\omega ,p\right)=F\left(\omega ,p\right),\quad \frac{\partial }{\partial t}\left[\mathcal{R}u\right]\left(0,\omega ,p\right)=0.\end{equation} \tag{ 10 }$$

Then, the solution of the Cauchy problem (9) and (10) is given by d'Alembert's formula,

$$\begin{equation*}\left[\mathcal{R}u\right]\left(t,\omega ,p\right)=\frac{1}{2}\left[F\left(\omega ,p+t\right)+F\left(\omega ,p-t\right)\right].\end{equation*}$$

Since f(x) is compactly supported in B(0, 1), its Radon projections F(ω, p) are supported on (−1, 1) in p. Therefore, for t = 2,

$$\begin{equation*}\left[\mathcal{R}u\right]\left(2,\omega ,p\right)=\frac{1}{2}F\left(\omega ,p-2\right),\quad p\in \left[1,3\right],\end{equation*}$$

or

$$\begin{equation*}F\left(\omega ,p\right)=2\left[\mathcal{R}u\right]\left(2,\omega ,p+2\right),\quad p\in \left[-1,1\right].\end{equation*}$$

Thus, once u(2, x) has been computed, one easily reconstructs the Radon projections F of f from $\mathcal{R}u$.

The methods for recovering a function from its Radon projections are well-known [19]. Consider the Hilbert transform $\mathcal{H}$ whose action a function of a 1D variable h(t) is given by the formula

$$\begin{equation*}\left[\mathcal{H}h\right]\left(s\right)=\underset{\mathbb{R}}{\int }\frac{h\left(t\right)}{s-t}\mathrm{d}t,\end{equation*}$$

where the integral is understood in the principal value sense. The adjoint ${\mathcal{R}}^{\#}$ of the Radon transform (6) is defined by its action on a function h(ω, p) supported on Z as follows:

$$\begin{equation*}\left[{\mathcal{R}}^{\#}h\right]\left(x\right)=\underset{{\mathbb{S}}^{d-1}}{\int }h\left(\omega ,\omega \cdot x\right)\mathrm{d}\omega .\end{equation*}$$

The inverse Radon transform ${\mathcal{R}}^{-1}$ is the left inverse of $\mathcal{R}$, and is in general not unique. For the purposes of the present paper we define ${\mathcal{R}}^{-1}$ as a composition of the differentiation in p, Hilbert transform $\mathcal{H}$ in the second variable (p), and the adjoint Radon transform ${\mathcal{R}}^{\#}$ as follows ([19], chapter 2, section 2):

$$\begin{equation}\left({\mathcal{R}}^{-1}F\right)\left(x\right)\equiv \frac{1}{2{\left(2\pi \right)}^{d-1}}\cdot \begin{cases}^{\frac{d-2}{2}}\left[{\mathcal{R}}^{\#}\mathcal{H}\frac{{\partial }^{d-1}}{\partial {p}^{d-1}}F\right]\left(x\right),\quad \hfill & d\;\text{is}\;\text{even,}\hfill \\ {\left(-1\right)}^{\frac{d-1}{2}}\left[{\mathcal{R}}^{\#}\frac{{\partial }^{d-1}}{\partial {p}^{d-1}}F\right]\left(x\right),\quad \hfill & d\;\text{is}\;\text{odd.}\hfill \end{cases}\end{equation} \tag{ 11 }$$

If $F\left(\omega ,p\right)=\left[\mathcal{R}f\right]\left(\omega ,p\right)$, then $f\left(x\right)=\left[{\mathcal{R}}^{-1}F\right]\left(x\right)$. From the numerical standpoint, the problem of reconstructing a function from its Radon projections has been studied extensively, and efficient and accurate algorithms for implementing (11) are well known by now [19].

It should be noted that the image reconstruction based on the Radon transform of the solution of the wave equation has been considered recently in [20]. In that work, a non-traditional data acquisition scheme yields the values of the whole acoustic field u(t, x), instead of boundary values g(t, z). The problem we are solving here is different, but our technique bears some resemblance of the methods of [20].

From the numerical standpoint it may be preferable to compute the derivative of the Radon projections:

$$\begin{equation*}\frac{\partial }{\partial p}F\left(\omega ,p\right)=\frac{\partial }{\partial p}2\left[\mathcal{R}u\right]\left(2,\omega ,p+2\right)=-2\left[\mathcal{R}\frac{\partial u}{\partial t}\right]\left(2,\omega ,p+2\right).\end{equation*}$$

Formulas (11) can be easily modified to process $\frac{\partial }{\partial p}F$ instead of F, since the Hilbert transform commutes with differentiation. However, for simplicity of the presentation we will concentrate on reconstruction of F.

2.2. Reconstruction from reduced data

In order to deal with problem 1, we first solve a modified initial/boundary problem in Q^ext as follows. For a fixed δ > 0 (see figure 1) we introduce a function $\psi \in {C}^{\infty }\left(\mathbb{R}\right)$ such that 0 ⩽ ψ ⩽ 1, ψ = 1 for τ ⩾ −δ, and ψ = 0 for τ ⩽ −2δ. We define the function h ∈ H^s((0, ∞) × ∂Ω) by setting

$$\begin{equation*}h\left(t,x\right)\equiv \psi \left({x}_{3}\right){g}^{\text{reduced}}\left(t,x\right)=\psi \left({x}_{3}\right)g\left(t,x\right).\end{equation*}$$

Consider now the exterior initial/boundary problem with homogeneous initial data and Dirichlet data h for the wave equation:

$$\begin{equation}\begin{aligned}\hfill {w}_{tt}\left(t,x\right)-{\Delta}w\left(t,x\right)& =0,\quad \left(t,x\right)\in {Q}^{\text{ext}},\hfill \\ \hfill w\left(0,x\right)& =0,\quad {w}_{t}\left(0,x\right)=0,\quad x\in {\mathbb{R}}^{d}{\backslash}\overline{{\Omega}},\hfill \\ \hfill w\left(t,z\right)& =h\left(t,z\right),\quad \left(t,z\right)\in \left(0,\infty \right){\times}\partial {\Omega}.\hfill \end{aligned}\end{equation} \tag{ 12 }$$

For h ∈ H^s((0, ∞) × ∂Ω), $s\in \mathbb{R}$, there exists a unique solution $w\in {H}_{\text{loc}}^{s}\left({Q}^{\text{ext}}\right)$. This result can be proved using the techniques developed in [21, section 24.1]. The local existence theorem (lemma 24.1.6) is restricted to non-negative s only because of a right-hand side term which vanishes in our case. Furthermore, if h ∈ L²((0, ∞) × ∂Ω), then the solution of our initial-boundary value problem is continuous in time, that is $w\in C\left(\left[0,\infty \right),{L}^{2}\left({\mathbb{R}}^{d}{\backslash}{\Omega}\right)\right)$ [22].

Solution w(t, x) serves as a crude approximation of the exact solution u(t, x) (compare this problem to problem (5)). Let us consider w(t, x) at t = 2. Due to the finite (unit) speed of propagation of waves, w(2, x) is finitely supported within the ball of radius 3 centered at the origin. Then, the (exterior) Radon transform $\left[\mathcal{R}w\right]\left(2,\omega ,p\right)$ of w(2, x) can be defined by (6) for 1 ⩽ |p| ⩽ 3. We only use a subset of these data and define $\tilde {F}\left(\omega ,p\right)$ as follows:

$$\begin{equation}\tilde {F}\left(\omega ,p\right)\equiv \left[\mathcal{R}w\right]\left(2,\omega ,p+2\right),\quad p\in \left[-1,1\right],\enspace \omega \in {\mathbb{S}}^{d-1}.\end{equation} \tag{ 13 }$$

When considered on the whole sphere ${\mathbb{S}}^{d-1}$ (in ω) the function $\tilde {F}\left(\omega ,p\right)$ is a rather poor approximation to F(ω, p). However, if one considers only the subset ${\mathbb{S}}_{+}^{d-1}\equiv \left\{\omega \in {\mathbb{S}}^{d-1}\enspace :\enspace \right.$$\left.{w}_{d}{ >}-\delta /2\right\}$, one finds that $\tilde {F}\left(\omega ,p\right)$ is equal to F(ω, p) modulo infinitely smooth function of (ω, p). This fact is proven in theorem 8 in section 4 of this paper.

Motivated by the symmetry relation (7), we define G(ω, p) as follows:

$$\begin{equation}G\left(\omega ,p\right)=\begin{cases}\tilde {F}\left(\omega ,p\right),\quad \hfill & {\omega }_{d}{\geqslant}0,p\in \left[-1,1\right],\hfill \\ \tilde {F}\left(-\omega ,-p\right),\quad \hfill & {\omega }_{d}{< }0,p\in \left[-1,1\right].\hfill \end{cases}\end{equation} \tag{ 14 }$$

We use G(ω, p) as an approximation for F(ω, p) and apply the inverse Radon transform ${\mathcal{R}}^{-1}$ (defined by (11)) to G(ω, p) to obtain an approximation $\tilde {f}\left(x\right)$ of f(x):

$$\begin{equation*}\tilde {f}\left(x\right)=\left[{\mathcal{R}}^{-1}G\right]\left(x\right).\end{equation*}$$

One can show that $\tilde {f}\left(x\right)=f\left(x\right)+\kappa \left(x\right)$, where κ(x) is a C^∞(B(0, 1)). This is the main theoretical result of the present paper, given by proposition 9 in section 4. Although our analysis is microlocal in nature, and does not yield a norm estimate on the error κ(x), our numerical simulations show that it is quite small in practical terms.

Our general reconstruction procedure consists of the following steps:

• (a)  
  Solve the initial/boundary value problem (12) on time interval [0, 2];
• (b)  
  Compute approximate projections $\tilde {F}$ using (13), and find G using (14);
• (c)  
  Compute ${\mathcal{R}}^{-1}G$ using (11), to obtain $\tilde {f}$.

Theoretical foundations of our technique are presented in sections 3 and 4 below.

Computational efficiency of our scheme depends on particular algorithms used on each of the steps. On the first step, known methods for solving the exterior initial/boundary value problems can be used, including, in particular, finite difference schemes [23] and boundary layer techniques [24]. The second step is straightforward. The last step, inversion of the Radon transform, is well-studied by now, and efficient implementations of the filtered backprojection algorithm (11) are known.

However, for the geometries that allow for separation of variables in the wave equation, it may be possible to have a faster implementations of steps 2 and 3. In section 5 we present a very fast algorithm applicable when the acquisition surface is a circular arc (in 2D). Finally, the last section of the paper illustrates our findings by the results of numerical simulations.

Let us analyze the propagation of singularities in the solution u(t, x) of the problem (1) and (2), with c(x) ≡ 1 and supp f = Ω₀. The solution for this Cauchy problem is given by the formula

$$\begin{align}\hfill u\left(t,x\right)& =\frac{1}{{\left(2\pi \right)}^{3}}{\int }_{{{\Omega}}_{0}}{\int }_{{\mathbb{R}}^{3}}{\mathrm{e}}^{\text{i}\left[\left(x-y\right)\cdot \xi +t\vert \xi \vert \right]}f\left(y\right)\enspace \mathrm{d}\xi \enspace \mathrm{d}y+\frac{1}{{\left(2\pi \right)}^{3}}{\int }_{{{\Omega}}_{0}}{\int }_{{\mathbb{R}}^{3}}{\mathrm{e}}^{\text{i}\left[\left(x-y\right)\cdot \xi -t\vert \xi \vert \right]}f\left(y\right)\enspace \mathrm{d}\xi \enspace \mathrm{d}y\hfill \\ \hfill & =\frac{1}{{\left(2\pi \right)}^{3}}{\int }_{{{\Omega}}_{0}}{\int }_{{\mathbb{R}}^{3}}\left[{\mathrm{e}}^{\text{i}t\vert \xi \vert }+{\mathrm{e}}^{-\mathrm{i}t\vert \xi \vert }\right]{\mathrm{e}}^{\text{i}\left(x-y\right)\cdot \xi }f\left(y\right)\enspace \mathrm{d}\xi \enspace \mathrm{d}y.\hfill \end{align} \tag{ 15 }$$

Consider the mapping $\mathcal{A}:f\to g$ where $g=u{\left\vert \right.}_{{\mathbb{R}}_{+}{\times}\partial {\Omega}}$. This operator is a sum of two Fourier integral operators. Since Ω₀ is a compact subset of B(0, 1), the operator $\mathcal{A}$ has a non-degenerate phase function and the mapping f → g is continuous from ${H}_{0}^{s}\left({{\Omega}}_{0}\right)$ into ${H}^{s}\left({\mathbb{R}}_{+}{\times}\partial {\Omega}\right)$ for $s\in \mathbb{R}$, see for example [10]. Hence, for $f\in {H}_{0}^{s}\left({{\Omega}}_{0}\right)$ we obtain h ∈ H^s((0, ∞) × ∂Ω), where h is the function introduced at the beginning of the previous subsection.

At first we will study the wave front set of the Dirichlet trace g. In what follows we abbreviate ${\mathbb{R}}_{+}\equiv \left(0,\infty \right)$. Since g is defined on ${\mathbb{R}}_{+}{\times}\partial {\Omega}$ its wave front set is a subset of ${T}^{{\ast}}\left({\mathbb{R}}_{+}{\times}\partial {\Omega}\right)$, the cotangent bundle of the surface ${\mathbb{R}}_{+}{\times}\partial {\Omega}$. The first result is of a more general nature.

Lemma 3. Let $p\in {H}^{s}\left({\mathbb{R}}^{d}\right)$ and let S be a smooth hypersurface in ${\mathbb{R}}^{d}$. Suppose that WF(p) ∩ N(S) = ∅ where N(S) is the conormal bundle of S. Then, the trace operator $p\to q=p{\left\vert \right.}_{S}$ extends to a continuous operator from ${H}^{s}\left({\mathbb{R}}^{d}\right)$ to H^s−1/2(S). Furthermore,

$$\begin{equation*}\mathrm{W}\mathrm{F}\left(q\right)\subset \left\{\left(x,\xi -\left(\xi \cdot {\nu }_{x}\right){\nu }_{x}\right)\in {T}^{{\ast}}S{\backslash}0\enspace :\enspace \left(x,\xi \right)\in \mathrm{W}\mathrm{F}\left(p\right)\right\}.\end{equation*}$$

Here ν_x is a unit normal vector at x ∈ S and ξ − (ξ ⋅ ν_x)ν_x is the orthogonal projection of ξ onto the cotangent space of S at x.

Proof. The first statement of the lemma is just the trace theorem in Sobolev spaces provided s > 1/2. The condition on the wave front set allows us to extend the trace theorem to all $s\in \mathbb{R}$. The concept of the wave front set is independent of the choice of coordinates. Locally, we may always assume that S = {x_d = 0}. In this case we will write x' = (x₁, ..., x_d−1) and correspondingly, we have ξ = (ξ', ξ_d). In local coordinates, the second statement of the lemma is

$$\begin{equation*}\mathrm{W}\mathrm{F}\left(q\right)=\left\{\left({x}^{\prime },{\xi }^{\prime }\right)\in {T}^{{\ast}}{\mathbb{R}}^{d-1}{\backslash}0\enspace :\enspace \left({x}^{\prime },0;{\xi }^{\prime },{\xi }_{d}\right)\in \mathrm{W}\mathrm{F}\left(p\right)\right\},\end{equation*}$$

and the conormal bundle is

$$\begin{equation*}N\left(S\right)=\left\{\left({x}^{\prime },0;0,{\eta }_{d}\right)\in {T}^{{\ast}}{\mathbb{R}}^{d}{\backslash}0\enspace :\enspace \left({x}^{\prime },0\right)\in S\right\}.\end{equation*}$$

Fix x' ∈ S. Since WF(p) ∩ N(S) = ∅ we know that there exists a $\varphi \in {C}_{0}^{\infty }\left({\mathbb{R}}^{d}\right)$ satisfying φ(x', 0) = 1, such that for all $N\in \mathbb{N}$ and ξ in some conic neighborhoods of the vectors $\left(0,\dots ,0,{\pm}1\right)\in {\mathbb{R}}^{d}$ we have

$$\begin{equation*}\vert \hat{\varphi p}\left(\xi \right)\vert {\leqslant}{C}_{N}{\left(1+\vert \xi \vert \right)}^{-N}.\end{equation*}$$

Furthermore, away from these conic neighborhoods, we have |ξ_d| ≲ |ξ'|. For s ⩽ 1/2, we have $\varphi p\in {H}_{\left(1,s-1\right)}\left({\mathbb{R}}^{d}\right)$. For a definition of these Sobolev spaces we refer to definition B.1.10 [21]. This local result can be made global and yields $p\in {H}_{\left(1,s-1\right)}\left({\mathbb{R}}^{d}\right)$. Using theorem B.1.11 [21], we infer that q ∈ H^s−1/2(S). Recall that for $p\in \mathcal{S}\left({\mathbb{R}}^{d}\right)$ we have

$$\begin{equation*}q\left({x}^{\prime }\right)=p\left({x}^{\prime },0\right)=\frac{1}{{\left(2\pi \right)}^{d}}{\int }_{{\mathbb{R}}^{d}}{\mathrm{e}}^{\text{i}{x}^{\prime }\cdot {\xi }^{\prime }}\hat{p}\left(\xi \right)\enspace \mathrm{d}\xi ,\end{equation*}$$

and hence,

$$\begin{equation}\hat{q}\left({\xi }^{\prime }\right)=\frac{1}{2\pi }{\int }_{\mathbb{R}}\hat{p}\left({\xi }^{\prime },{\xi }_{d}\right)\enspace \mathrm{d}{\xi }_{d}.\end{equation} \tag{ 16 }$$

By continuity this last formula extends to Fourier transforms which are integrable in ξ_d with values in ${\mathcal{S}}^{\prime }\left({\mathbb{R}}^{d}\right)$. Now suppose that the points (x', 0; ξ', ξ_d) ∉ WF(p) for all ${\xi }_{d}\in \mathbb{R}$. By the definition of the wave front set, there exists a function $\varphi \in {C}_{0}^{\infty }\left({\mathbb{R}}^{d}\right)$ such that φ(x', 0) = 1 and a conic neighborhood of U(ξ') in ${\mathbb{R}}^{d-1}$ such that for all $N\in \mathbb{N}$

$$\begin{equation*}\left(\hat{\varphi p}\right)\left({\eta }^{\prime },{\xi }_{d}\right){\leqslant}{C}_{N}{\left(1+\vert {\eta }^{\prime }{\vert }^{2}+\vert {\xi }_{d}{\vert }^{2}\right)}^{-N}\quad \;\text{for}\;\text{all}\enspace {\eta }^{\prime }\in U\left({\xi }^{\prime }\right)\enspace \text{and}\enspace {\xi }_{d}\in \mathbb{R}.\end{equation*}$$

Now define $\psi \in {C}_{0}^{\infty }\left(S\right)$ by ψ(y') = φ(y', 0). Then we know that for all η' ∈ U(ξ')

$$\begin{equation*}\vert \left(\hat{\psi q}\right)\left({\eta }^{\prime }\right)\vert {\leqslant}\frac{1}{2\pi }{\int }_{\mathbb{R}}\vert \hat{\varphi p}\left({\eta }^{\prime },{\xi }_{d}\right)\vert \enspace \mathrm{d}{\xi }_{d}{\leqslant}\frac{{C}_{N}}{2\pi }{\int }_{\mathbb{R}}{\left(1+\vert {\eta }^{\prime }{\vert }^{2}+{\xi }_{d}^{2}\right)}^{-N}\enspace \mathrm{d}{\xi }_{d}{\leqslant}{C}_{N}^{\prime }{\left(1+\vert {\xi }^{\prime }{\vert }^{2}\right)}^{-N+1/2},\end{equation*}$$

where ${C}_{N}^{\prime }={C}_{N}{\int }_{\mathbb{R}}{\left(1+{s}^{2}\right)}^{-N}\enspace \mathrm{d}s/\left(2\pi \right)$. Hence, we conclude that (x', ξ') ∉ WF(q) and conclude that

$$\begin{equation*}\mathrm{W}\mathrm{F}\left(q\right)\subset \left\{\left({x}^{\prime },{\xi }^{\prime }\right)\in {T}^{{\ast}}{\mathbb{R}}^{d-1}{\backslash}0\enspace :\enspace \left({x}^{\prime },0;{\xi }^{\prime },{\xi }_{d}\right)\in \mathrm{W}\mathrm{F}\left(p\right)\right\}.\end{equation*}$$

□

The wave front set of the solution u to (15) can be described very precisely. Over every point (x, η) ∈ WF(f) there are two null bicharacteristics, that is lines in $\gamma \left(s\right)=\left(t\left(s\right),x\left(s\right);\tau \left(s\right),\eta \left(s\right)\right)\subset {T}^{{\ast}}{\mathbb{R}}^{d+1}$, $s\in \mathbb{R}$, which satisfy the following initial value problem

$$\begin{equation*}\dot {t}=-2\tau ,\quad \dot {x}=2\eta ,\quad \dot {\tau }=0,\quad \dot {\eta }=0,\quad t\left(0\right)=0,\quad x\left(0\right)=x,\quad \tau \left(0\right)={\pm}\vert \eta \vert ,\quad \eta \left(0\right)=\eta .\end{equation*}$$

Note that the initial condition for τ is inferred from the fact that null bicharacteristics are in the characteristic variety. One can use t as a parameter and obtain the following solutions

$$\begin{equation}{\gamma }_{+}\left(t\right)=\left(t,x+t\eta /\vert \eta \vert ;-\vert \eta \vert ,\eta \right)\quad \text{and}\quad {\gamma }_{-}\left(s\right)=\left(t,x-t\eta /\vert \eta \vert ;\vert \eta \vert ,\eta \right)\enspace ,\enspace t\in \mathbb{R}.\end{equation} \tag{ 17 }$$

Note that γ₊ has the direction vector (|η|, η) for the first d + 1 components where as γ₋ has the direction (|η|, −η).

Since we are interested in the forward problem, we will restrict ourselves to t ⩾ 0. The wave front set of u consists exactly of all bicharacteristics whose initial data are in the wave front set of f, see for example [21, chapter XXIII]. The next result describes the wave front set of the trace of the wave u along ∂Ω.

Proposition 4. Let f ∈ H^s(Ω₀) for some $s\in \mathbb{R}$ and let $u\in C\left(\left[0,\infty \right),{H}^{s}\left({\mathbb{R}}^{d}\right)\right)$ be the unique solution to (1) and (2) and $g=u{\left\vert \right.}_{{\mathbb{R}}_{+}{\times}\partial {\Omega}}$. Then,

$$\begin{equation*}\begin{aligned}\hfill \mathrm{W}\mathrm{F}\left(g\right)\subset & \left\{\left(\vert y-x\vert ,y;\mp \alpha \vert y-x\vert ,{\pm}\alpha \left[y-x-\left(\left(y-x\right)\cdot {\nu }_{y}\right){\nu }_{y}\right]\right)\in {T}^{{\ast}}\left({\mathbb{R}}_{+}{\times}\partial {\Omega}\right){\backslash}0:\right.\hfill \\ \hfill & \left.y\in \partial {\Omega},\alpha { >}0\quad \;\text{and}\;\quad \left(x,{\pm}\left[y-x\right]\right)\in \mathrm{W}\mathrm{F}\left(f\right)\right\}.\hfill \end{aligned}\end{equation*}$$

Furthermore, every point in WF(g) is a hyperbolic point. This is to say that (t, y; τ, ξ) ∈ WF(g) implies |τ| > |ξ|.

Proof. We know that

$$\begin{equation}\mathrm{W}\mathrm{F}\left(u\right)=\left\{\left(t,x{\pm}t\eta /\vert \eta \vert ,\mp \vert \eta \vert ,\eta \right)\enspace :\enspace t{\geqslant}0\quad \;\text{and}\;\quad \left(x,\eta \right)\in \mathrm{W}\mathrm{F}\left(f\right)\right\}\subset {T}^{{\ast}}{\mathbb{R}}_{+}^{d+1}{\backslash}0\end{equation} \tag{ 18 }$$

and observe that $\mathrm{W}\mathrm{F}\left(u\right)\cap N\left({\mathbb{R}}_{+}{\times}\partial {\Omega}\right)=\varnothing$. Hence, the previous lemma is applicable.

If y ∈ ∂Ω, then y = x ± tη/|η| implies that ±α(y − x) = η for some α > 0 and t = |y − x|. Hence,

$$\begin{align*}\hfill & \mathrm{W}\mathrm{F}\left(u\right)\cap {T}_{{\mathbb{R}}_{+}{\times}\partial {\Omega}}^{{\ast}}{\mathbb{R}}^{d+1}\hfill \\ \hfill & =\left\{\left(\vert y-x\vert ,y;\mp \alpha \vert y-x\vert ,{\pm}\alpha \left(y-x\right)\right)\enspace :\enspace y\in \partial {\Omega},\alpha { >}0\quad \;\text{and}\;\quad \left(x,{\pm}\left[y-x\right]\right)\in \mathrm{W}\mathrm{F}\left(f\right)\right\}\hfill \end{align*}$$

Now we apply lemma 3 and the proof of the first statement is finished.

For the second statement, let (t, y; τ, ξ) ∈ WF(g). Then there exists a point (t, y; τ, η) ∈ WF(u) where

$$\begin{equation*}\xi =\eta -\left(\eta \cdot {\nu }_{y}\right){\nu }_{y}\quad \;\text{and}\;\quad \vert \tau \vert =\vert \eta \vert .\end{equation*}$$

Since the support Ω₀ is a compact subset of the convex set Ω, no bicharacteristic tangential to ${\mathbb{R}}_{+}{\times}\partial {\Omega}$ can meet Ω₀. Hence, $\eta \notin {T}_{y}^{{\ast}}\partial {\Omega}$ and thus |ξ| < |η| = |τ|. □

Remark 5. The inclusion in proposition 4 can be replaced by an equal sign. This follows from the fact that the operator $\mathcal{A}$ is a sum of two elliptic Fourier integral operators.

Now comes the second step. The unique solution w(t, x) to the exterior problem (12) has the following interesting property. The wave fronts of u and w are the same in the exterior upper half space

$$\begin{equation*}{D}_{+}=\left\{x\in {\mathbb{R}}^{d}\enspace :\enspace {x}_{d}{\geqslant}-\delta \right\}.\end{equation*}$$

More precisely, we have

Theorem 6. For $f\in {H}_{0}^{s}\left({{\Omega}}_{0}\right)$, suppose that u is the solution to the Cauchy problem (1) and (2) and that w is the solution to the exterior initial-boundary value problem (12). Then

$$\begin{equation*}\mathrm{W}\mathrm{F}\left(u-w\right)\cap {T}^{{\ast}}\left\{\left(t,x\right)\in {Q}^{\text{ext}}\enspace :\enspace {x}_{d}{\geqslant}-\delta \right\}=\varnothing \end{equation*}$$

and

$$\begin{equation*}\mathrm{W}\mathrm{F}\left(u-w\right)\left(2,\cdot \right)\cap {T}^{{\ast}}{D}_{+}=\varnothing .\end{equation*}$$

Proof. Let v be the solution of the exterior Cauchy–Dirichlet problem with zero initial data and Dirichlet data (1 − ψ)g, i.e. v = u − w. We will show that

$$\begin{equation*}\mathrm{W}\mathrm{F}\left(v\right)\cap {T}^{{\ast}}\left\{\left(t,x\right)\in {Q}^{\text{ext}}\enspace :\enspace {x}_{d}{\geqslant}-\delta \right\}=\varnothing ,\end{equation*}$$

i.e. the function v is infinitely often differentiable in the upper half space. According to theorem 2.1 [25], the singularities of v depend on the wave front set of the Dirichlet data (1 − ψ)g. Since ψ ∈ C^∞ we know that $\mathrm{W}\mathrm{F}\left(\left(1-\psi \right)g\right)\subset \mathrm{W}\mathrm{F}\left(g\right)\cap {T}^{{\ast}}\left({\mathbb{R}}_{+}{\times}\partial {{\Omega}}_{-}\right)$, where $\partial {{\Omega}}_{-}=\left\{x\in \partial {\Omega}\enspace :\enspace {x}_{d}{< }-\delta \right\}$.

According to proposition 4 the wave front set of g depends only on the wave front set of f. In particular all points in WF(g) are hyperbolic points. Over each hyperbolic points there are exactly two bicharacteristics, an incoming one and an outgoing one. Since the initial data are zero, the solution is identically zero along all incoming bicharacteristics.

Let now (t, y; τ, ξ) ∈ WF(g) and denote the exterior unit normal vector to Ω at y by ν_y. Then, according to proposition 4, we have either τ > 0 (bicharacteristic γ₋) or τ < 0 (bicharacteristic γ₊). The spatial frequency of the outgoing bicharacteristic is $\eta =\xi +\sqrt{{\tau }^{2}-\vert \xi {\vert }^{2}}{\nu }_{y}$ in the case γ₊ and $\eta =\xi -\sqrt{{\tau }^{2}-\vert \xi {\vert }^{2}}{\nu }_{y}$ in the case of γ₋, compare with formula (17).

We will consider only the case τ > 0, the case τ < 0 is similar. Because of proposition 4 there exists a point x ∈ Ω₀ and α > 0 such that

$$\begin{equation*}t=\vert y-x\vert ,\quad \tau =\alpha \vert y-x\vert ,\quad \xi =-\alpha \left[y-x-\left(\left(y-x\right)\cdot {\nu }_{y}\right){\nu }_{y}\right].\end{equation*}$$

Compute now

$$\begin{equation*}{\tau }^{2}-\vert \xi {\vert }^{2}={\alpha }^{2}{\left[\left(y-x\right)\cdot {\nu }_{y}\right]}^{2}\end{equation*}$$

and observe that the convexity of Ω implies that

$$\begin{equation*}\sqrt{{\tau }^{2}-\vert \xi {\vert }^{2}}=\alpha \left[\left(y-x\right)\cdot {\nu }_{y}\right].\end{equation*}$$

The outgoing bicharacteristic γ₋ starting at the point (t, y) must have the direction

$$\begin{equation*}\left(\tau ,\sqrt{{\tau }^{2}-\vert \xi {\vert }^{2}}{\nu }_{y}-\xi \right)=\alpha \left(\vert y-x\vert ,y-x\right),\end{equation*}$$

compare with (17). Since x ∈ Ω₀ and $y\in \partial {{\Omega}}_{-}$, the projection of the bicharacteristic γ₋ into space-time is a straight line which cannot meet {(t, x) ∈ Q^ext: x_d ⩾ −δ} for t > 0.

The second statement follows from the first one and lemma 3. □

We have defined the Radon transform earlier, by equation (6).

Proposition 7. The Radon transform can be extended to a bounded linear operator from ${\mathcal{E}}^{\prime }\left({\mathbb{R}}^{d}\right)$ to ${\mathcal{D}}^{\prime }\left({\mathbb{S}}^{d-1}{\times}\mathbb{R}\right)$ and

$$\begin{align}\hfill \mathrm{W}\mathrm{F}\left(\mathcal{R}f\right)\subset & \left\{\left({\pm}\omega ,{\pm}\omega \cdot x;\mp \alpha \left[x-\left(x\cdot \omega \right)\omega \right],{\pm}\alpha \right)\in {T}^{{\ast}}\left({\mathbb{S}}^{d-1}{\times}\mathbb{R}\right){\backslash}0\right.\hfill \\ \hfill & \left.\enspace :\omega \in {\mathbb{S}}^{d-1},\alpha { >}0\quad \;\text{and}\;\quad \left(x,\omega \right)\in \mathrm{W}\mathrm{F}\left(f\right)\right\}.\hfill \end{align} \tag{ 19 }$$

For d = 2 this is theorem 3.1 in [26].

Proof. A large part of the proof is an application of theorem 8.2.13 in [27]. We use local coordinates ω = φ(y'), ${y}^{\prime }\in {Y}^{\prime }\subset {\mathbb{R}}^{d-1}$ where we demand that φ(Y') is at most a hemisphere. The kernel of the Radon transform is

$$\begin{equation*}K\left(y,x\right)=\delta \left(\varphi \left(z\right)\cdot x-p\right),\end{equation*}$$

where we abbreviate $y=\left({y}^{\prime },p\right)\in {\mathbb{R}}^{d}$ and see that the kernel is a composition of the delta function and the function m(y, x) = φ(y') ⋅ x − p. Note that $K\in {\mathcal{D}}^{\prime }\left(Y{\times}X\right)$ where $Y={Y}^{\prime }{\times}\mathbb{R}$ and $X={\mathbb{R}}^{d}$. Recall that $\mathrm{W}\mathrm{F}\left(\delta \right)=\left\{\left(0,\tau \right)\in {T}^{{\ast}}\mathbb{R}{\backslash}0\enspace :\enspace \tau \in \mathbb{R}\right\}$. The wave front set of this kernel is analyzed using theorem 8.2.4 in [27]. We have

$$\begin{equation*}{m}^{\prime }\left(y,x\right)=\left[\begin{bmatrix}{c}\hfill {\varphi }^{\prime }{\left({y}^{\prime }\right)}^{\top }x\hfill \\ -1\hfill \\ \varphi \left({y}^{\prime }\right)\hfill \end{bmatrix}\right]\end{equation*}$$

and hence,

$$\begin{equation}\mathrm{W}\mathrm{F}\left(K\right)\subset \left\{\left(\left[\begin{bmatrix}{c}{y}^{\prime }\hfill \\ \hfill \varphi \left({y}^{\prime }\right)\cdot x\hfill \end{bmatrix}\right],x;\left[\begin{bmatrix}{c}\hfill \tau {\varphi }^{\prime }{\left({y}^{\prime }\right)}^{\top }x\hfill \\ -\tau \hfill \end{bmatrix}\right],\tau \varphi \left({y}^{\prime }\right)\right)\in {T}^{{\ast}}\left(Y{\times}X\right){\backslash}0\enspace :\enspace \tau \in \mathbb{R},{y}^{\prime }\in {Y}^{\prime },x\in {\mathbb{R}}^{d}\right\}.\end{equation} \tag{ 20 }$$

Here φ' is the Jacobian matrix of φ. Now theorem 8.2.13 in [27] is used to finish the proof. At first note that

$$\begin{equation*}{\mathrm{W}\mathrm{F}}^{\prime }{\left(K\right)}_{X}=\left\{\left(x,\xi \right)\in {T}^{{\ast}}X{\backslash}0\enspace :\enspace \left(y\cdot x,x;0,-\xi \right)\in \mathrm{W}\mathrm{F}\left(K\right)\quad \;\text{for}\enspace \text{some}\enspace y\in Y\right\}=\varnothing ,\end{equation*}$$

since the first fiber component in (20) cannot be zero. Hence, the Radon transform can be continuously extended to a linear operator from ${\mathcal{E}}^{\prime }\left(X\right)$ into ${\mathcal{D}}^{\prime }\left(Y\right)$. Similarly, we have

$$\begin{equation*}\mathrm{W}\mathrm{F}{\left(K\right)}_{Y}=\left\{\left(y,\eta \right)\in {T}^{{\ast}}X{\backslash}0\enspace :\enspace \left(y,x;\eta ,0\right)\in \mathrm{W}\mathrm{F}\left(K\right)\quad \;\text{for}\enspace \text{some}\enspace y\in Y\right\}=\varnothing ,\end{equation*}$$

and hence,

$$\begin{align}\hfill & \mathrm{W}\mathrm{F}\left(\mathcal{R}\left(f\right)\right)\subset {\mathrm{W}\mathrm{F}}^{\prime }\left(K\right){\circ}\mathrm{W}\mathrm{F}\left(f\right)=\left\{\left(y,\eta \right)\enspace :\enspace \left(y,x;\eta ,-\xi \right)\in \mathrm{W}\mathrm{F}\left(K\right)\enspace \text{for}\enspace \text{some}\enspace \left(x,\xi \right)\in \mathrm{W}\mathrm{F}\left(f\right)\right\}\hfill \\ \hfill & \quad =\left\{\left(\left[\begin{bmatrix}{c}{y}^{\prime }\hfill \\ \hfill \varphi \left({y}^{\prime }\right)\cdot x\hfill \end{bmatrix}\right],\left[\begin{bmatrix}{c}\hfill -\tau {\varphi }^{\prime }{\left({y}^{\prime }\right)}^{\top }x\hfill \\ \tau \hfill \end{bmatrix}\right]\right)\in {T}^{{\ast}}\left(Y\right){\backslash}0\enspace :\enspace \tau \in \mathbb{R},{y}^{\prime }\in {Y}^{\prime },\left(x,\varphi \left({y}^{\prime }\right)\right)\in \mathrm{W}\mathrm{F}\left(f\right)\right\}.\hfill \end{align} \tag{ 21 }$$

Finally we return from the local coordinate system in Y to $\left(\omega ,p\right)\in {\mathbb{S}}^{d-1}{\times}\mathbb{R}$. Here it is important to note that the co-vector ${\varphi }^{\prime }{\left({y}^{\prime }\right)}^{\top }x\in {T}_{{y}^{\prime }}^{{\ast}}{Y}^{\prime }$ gets mapped to $x-\left(\omega \cdot x\right)\omega \in {T}_{\omega }^{{\ast}}{\mathbb{S}}^{d-1}$. Furthermore, if (x, ω) ∈ WF(f), then (x, τω) ∈ WF(f) for τ > 0 and (x, −τω) ∈ WF(f) for all τ < 0. □

Actually, a more precise result holds true, there is equality in (19). This can be inferred from the theory of Fourier integral operators and is explained in [28, p 59].

The Radon transforms of the solutions u(t, x) to the Cauchy problem (1) and (2) and w(t, x) to the Cauchy–Dirichlet problem (12) carry t as a parameter. Recall ${\mathbb{S}}_{+}^{d-1}=\left\{\omega \in {\mathbb{S}}^{d-1}\enspace :\enspace {\omega }_{d}{ >}-\delta /2\right\}$.

Theorem 8. Let $f\in {H}_{0}^{s}\left({{\Omega}}_{0}\right)$ and suppose that u is the solution to the Cauchy problem (1) and (2). Let w be the solution to the exterior initial-boundary value problem (12). Then,

$$\begin{equation*}\mathrm{W}\mathrm{F}\left(\left[\mathcal{R}\left(u-w\right)\right]\left(2,\cdot \right)\right)\cap {T}^{{\ast}}\left({\mathbb{S}}_{+}^{d-1}{\times}\left(1,3\right)\right)=\varnothing .\end{equation*}$$

Proof. As before, let v = u − w. We have to show that $\mathrm{W}\mathrm{F}\left(\left[\mathcal{R}v\right]\left(2,\cdot \right)\right)\cap {T}^{{\ast}}\left({\mathbb{S}}_{+}^{d-1}{\times}\left(1,3\right)\right)=\varnothing$. Let $\omega \in {\mathbb{S}}_{+}^{d-1}$. proposition 7 can be rephrased as follows: the wave front set of the Radon transform of v(2, ⋅) can have points with base point (ω, p) only if the wave front set of v(2, ⋅) has points of the form (x, ±ω) for some $x\in {\mathbb{R}}^{d}$ satisfying ω ⋅ x = p. According to theorem 6, the wave front sets of v(2, ⋅) is empty over the set D₊ whereas over ${D}_{-}{:=}\left\{x\in {\mathbb{R}}^{d}{\backslash}\overline{{\Omega}}\enspace :\enspace {x}_{d}{< }-\delta \right\}$ the wave front set of v(2, ⋅) is a subset of the wave front set of u(2, ⋅).

However, we can show that the wave front set of u over D₋ has no points of the form (x, ±ω) with x satisfying ω ⋅ x = p for some p > − 1. The wave front set of u is given in (18). If $\left(x,{\pm}\omega \right)\in \mathrm{W}\mathrm{F}\left(u\left(2,\cdot \right)\right)$, then (y, ±ω) ∈ WF(f) where x = y ± 2ω. Note that x = y + 2ω is not possible since supp f = Ω₀ implies y_d ⩾ 0 and thus x_d > −δ. Hence, x ⋅ ω = y ⋅ ω − 2 < −1. □

Recall the definitions of $\tilde {F}$ and G given in (13) and (14). We will be interested in the error

$$\begin{equation*}{E}_{G}\left(\omega ,p\right)=F\left(\omega ,p\right)-G\left(\omega ,p\right).\end{equation*}$$

For $\omega \in {\mathbb{S}}_{+}^{d-1}$ and p ∈ (−1, 1) we have just proved that $F\left(\omega ,p\right)-\tilde {F}\left(\omega ,p\right)$ is an infinitely smooth function. Furthermore, by symmetry, for ω_d < δ/2 and p ∈ (−1, 1), we infer that $F\left(\omega ,p\right)-\tilde {F}\left(-\omega ,-p\right)$ is smooth.

Let H denote the Heaviside function of a one-dimensional variable

$$\begin{equation*}H\left(s\right)\equiv \begin{cases}1,\quad \hfill & s{\geqslant}0,\hfill \\ 0,\quad \hfill & s{< }0.\hfill \end{cases}\end{equation*}$$

Then,

$$\begin{equation}{E}_{G}\left(\omega ,p\right)=F\left(\omega \right)-H\left({\omega }_{d}\right)\tilde {F}\left(\omega ,p\right)-\left[1-H\left({\omega }_{d}\right)\right]\tilde {F}\left(-\omega ,-p\right)\end{equation} \tag{ 22 }$$

and this function defined on ${\mathbb{S}}^{d-1}{\times}\left(-1,1\right)$ is an infinitely smooth function except for a discontinuity on the set $\mathbb{E}{\times}\left(-1,1\right)$ where $\mathbb{E}$ is the equator $\mathbb{E}\equiv \left\{\omega \enspace :\enspace \omega \in {\mathbb{S}}^{d-1},{\omega }_{d}=0\right\}$. The next result shows that these singularities disappear, once the inverse Radon transform is applied. In order to apply the inverse Radon transform we need to extend E_G to all $p\in \mathbb{R}$. This is done by setting G ≡ 0 for all |p| > 1. Abusing notation, we will denote the resulting extension by E_G as well.

Proposition 9. The inverse Radon transform ${\mathcal{R}}^{-1}{E}_{G}$ restricted to the open unit ball is a smooth function κ ∈ C^∞(B(0, 1)).

Proof. The function E_G has two locations of singularities: one is on $\mathbb{E}{\times}\left(-1,1\right)$, due to (22) and the other is at |p| = 1 due to the extension by zero outside of p ∈ [−1, 1].

We discuss at first $\mathrm{W}\mathrm{F}\left({E}_{G}\right)\cap {T}^{{\ast}}\left({\mathbb{S}}^{d-1}{\times}\left(-1,1\right)\right)$. According to formula (22), we have

$$\begin{equation*}\begin{aligned}\hfill & \mathrm{W}\mathrm{F}\left({E}_{G}\right)\cap {T}^{{\ast}}\left({\mathbb{S}}^{d-1}{\times}\left(-1,1\right)\right)\subset \left\{\left(\left[\begin{bmatrix}{c}\hfill \omega \hfill \\ \hfill p\hfill \end{bmatrix}\right],\left[\begin{bmatrix}{c}\hfill \xi \hfill \\ \hfill 0\hfill \end{bmatrix}\right]\right)\in {T}^{{\ast}}\left({\mathbb{S}}^{d-1}{\times}\left(-1,1\right)\right){\backslash}0\enspace :\right.\hfill \\ \hfill & \left.{\omega }_{d}=0,\xi \perp \dot{c}\left(0\right),c\left(0\right)=\omega ,\quad \;\text{for}\enspace \text{all}\enspace \enspace c\left(t\right)\in \mathbb{E}\right\}.\hfill \end{aligned}\end{equation*}$$

Here $c:\left(-\alpha ,\alpha \right)\to \mathbb{E}$ denotes regular curves in $\mathbb{E}$. In view of (19) we see that this set has an empty intersection with the wave front set of the Radon transform of any distribution.

From the theory of elliptic Fourier integral operators [29, theorem 4.2.5] one knows that the inverse transform ${\mathcal{R}}^{-1}:{\mathcal{E}}^{\prime }\left({\mathbb{S}}^{d-1}{\times}\mathbb{R}\right)\to {\mathcal{D}}^{\prime }\left({\mathbb{R}}^{d}\right)$ has the same canonical relation as does $\mathcal{R}$, only with the variables switched around. The application of the inverse Radon transform will ignore those points in the wave front set of E_G.

In view of (19), the singularities of E_G at |p| = 1 will correspond to singularities of its inverse Radon transform outside of the unit ball.

But this means that $\mathrm{W}\mathrm{F}\left({\mathcal{R}}^{-1}\left[{E}_{G}\right]\right)\cap {T}^{{\ast}}B\left(0,1\right)=\varnothing$. □

The three-step image reconstruction procedure proposed in this paper is detailed at the end of section 2.2. Practicality of this technique depends, in particular, on the computational complexity of steps 1 and 2. (Efficient algorithms for implementing step 3 are well known). Let us assume that the domain Ω is discretized using m^d nodes, and the number of time steps is $\mathcal{O}\left(m\right)$. The optimal number of projections is then $\mathcal{O}\left({m}^{d-1}\right)$ with the number of samples in a projection equal to $\mathcal{O}\left(m\right)$. If the first step is implemented using finite differences, it will have computational complexity $\mathcal{O}\left({m}^{d+1}\right)$ floating point operations (flops). A straightforward computation of the Radon projections requires $\mathcal{O}\left({m}^{d+2}\right)$ flops.

However, for special geometries it might be possible to design significantly faster algorithms. In the present section we develop an asymptotically fast ($\mathcal{O}\left({m}^{2}\mathrm{log}m\right)$ flops) algorithm applicable in 2D when the acquisition curve is a circular arc.

To be precise, we consider the following geometry. The region Ω is the unit disk B(0, 1) centered at the origin, Ω₀ is the upper half of B(0, 1), i.e.

$$\begin{equation*}{{\Omega}}_{0}\equiv \left\{x:\vert x\vert {< }1\quad \text{and}\quad {x}_{d}{ >}0\right\}.\end{equation*}$$

The data acquisition surface Γ is an arc of a unit circle

$$\begin{equation*}{\Gamma}\equiv \left\{z:\vert z\vert =1\quad \text{and}\quad {z}_{d}{ >}-2\delta .\right\}\end{equation*}$$

The geometry is shown in figure 2, together with the phantom we use in the numerical simulations described further below (for which we chose 2δ to be equal to arctan(10°)).

Zoom In Zoom Out Reset image size

Due to the circular symmetry of Ω, solving the exterior boundary value problem (12) commutes with rotation. In addition, it commutes with the shifts in time variable. It is well known that the Radon transform commutes with rotation, too. Thus, one can expect that the operator that maps boundary data h(ω, t) into the Radon projections is a convolution, and that the Fourier transform in t combined with Fourier series in the angular variable would diagonalize it. This is indeed the case as shown below.

We start by noticing that the wave equation separates in the polar coordinate system. In particular, the following combinations W_k(t, λ, r, θ) of functions represent radiating solutions of the wave equation (propagating away from the origin):

$$\begin{equation*}{W}_{k}\left(t,\lambda ,r,\theta \right)\equiv \frac{{H}_{\vert k\vert }^{\left(1\right)}\left(\lambda r\right)\mathrm{exp}\left(\mathrm{i}k\theta \right)\mathrm{exp}\left(-\mathrm{i}\lambda t\right)}{{H}_{\vert k\vert }^{\left(1\right)}\left(\lambda \right)},\quad \lambda \in \mathbb{R}{\backslash}0,\enspace k\in \mathbb{Z}.\end{equation*}$$

where ${H}_{k}^{\left(1\right)}\left(s\right)$ is the Hankel function of the first kind of order k, λ is any non-zero frequency, and r and t are the polar coordinates corresponding to the variable x, i.e. x = r(cos θ, sin θ). Let us extend W_k(t, λ, r, θ) to λ = 0 by continuity (see [30], section 10.7):

$$\begin{equation*}{W}_{k}\left(t,0,r,\theta \right)\equiv \frac{1}{{r}^{\vert k\vert }}\mathrm{exp}\left(\mathrm{i}k\theta \right),\quad k\in \mathbb{Z}.\end{equation*}$$

At λ = 0 functions W_k(t, 0, r, θ) are stationary solutions of the wave equations (i.e., harmonic functions), bounded at infinity. The values of W_k's on the cylinder r = 1, θ ∈ [0, 2π], $t\in \mathbb{R}$ are functions exp(ikθ) exp(−iλt); they represent a Fourier basis.

We will use functions W_k(t, λ, r, θ) to solve the forward problem (12). Let us extend the boundary condition h(t, x(1, θ)) in the time variable by 0 to all t < 0. Since solutions of the wave equation are causal, and we are only interested in w(t, x) at t = T, we can smoothly cut-off h(t, x(1, θ)) to 0 on the interval [T, T + b], and define $h\left(\right.t,x\left(1,\theta \right)$ to be 0 for all t > T + b. There is enough freedom in defining this cut-off to enforce the condition

$$\begin{equation*}\underset{0}{\overset{2\pi }{\int }}\underset{\mathbb{R}}{\int }h\left(t,x\left(1,\theta \right)\right)\mathrm{d}t\enspace \mathrm{d}\theta =0.\end{equation*}$$

Now, solution of problem (12) can be formally represented as

$$\begin{equation}w\left(t,x\right)\equiv {w}^{{\ast}}\left(t,r,\theta \right)=\underset{\mathbb{R}}{\int }\left[\sum _{k=-\infty }^{\infty }{b}_{k}\left(\lambda \right){W}_{k}\left(t,\lambda ,r,\theta \right)\right]\mathrm{d}\lambda .\end{equation} \tag{ 23 }$$

where coefficients a_k(λ) are to be determined from the boundary condition

$$\begin{equation}{w}^{{\ast}}\left(t,1,\theta \right)=h\left(t,x\left(1,\theta \right)\right),\end{equation} \tag{ 24 }$$

on $\mathbb{R}{\times}{\mathbb{S}}^{1}$, which yields

$$\begin{equation}{b}_{k}\left(\lambda \right)=\frac{1}{{\left(2\pi \right)}^{2}}\underset{0}{\overset{2\pi }{\int }}\left[\underset{\mathbb{R}}{\int }h\left(t,x\left(1,\theta \right)\right)\mathrm{exp}\left(\mathrm{i}\lambda t\right)\mathrm{d}t\right]\mathrm{exp}\left(-\mathrm{i}k\theta \right)\mathrm{d}\theta .\end{equation} \tag{ 25 }$$

We notice, for future use, that each b_k(λ) is continuous (in fact, real-analytic) due to the bounded support of h in time variable.

From the numerical standpoint it is somewhat easier to find instead of the Radon projections their derivative in p. To that end we will use $\frac{\partial }{\partial t}{w}^{{\ast}}\left(t,r,\theta \right)$; we only need its value at t = T:

$$\begin{align}\hfill \frac{\partial w}{\partial t}\left(T,x\right)& =\frac{\partial {w}^{{\ast}}}{\partial t}\left(T,r,\theta \right)=-i\underset{\mathbb{R}}{\int }\left[\sum _{k=-\infty }^{\infty }{b}_{k}\left(\lambda \right)\lambda {W}_{k}\left(T,\lambda ,r,\theta \right)\right]\mathrm{d}\lambda .\hfill \\ \hfill & =-i\underset{\mathbb{R}{\backslash}0}{\int }\left[\sum _{k=-\infty }^{\infty }{b}_{k}\left(\lambda \right)\lambda \frac{{H}_{\vert k\vert }^{\left(1\right)}\left(\lambda r\right)}{{H}_{\vert k\vert }^{\left(1\right)}\left(\lambda \right)}\mathrm{exp}\left(\mathrm{i}k\theta \right)\mathrm{exp}\left(-\mathrm{i}\lambda T\right)\right]\mathrm{d}\lambda .\hfill \end{align} \tag{ 26 }$$

In the above equation we are able to exclude λ = 0 from the integration domain since the integrand is zero at this point.

Our goal is to approximate $\frac{\partial }{\partial p}F\left(\omega ,p\right)$ by the Radon transform of $\frac{\partial w}{\partial t}$:

$$\begin{equation*}\frac{\partial }{\partial p}F\left(\omega ,p\right)\approx -2\left[\mathcal{R}\frac{\partial w}{\partial t}\right]\left(2,\omega ,p+2\right)\end{equation*}$$

However, computing $\frac{\partial w}{\partial t}\left(T,x\right)$ by using equation (26) and by numerically evaluating its Radon transform would be relatively slow. Instead, we take advantage of the fact that the line integrals of the circular wave functions ${H}_{\vert k\vert }^{\left(1\right)}\left({\lambda }_{l}r\right)\mathrm{exp}\left(\mathrm{i}k\theta \right)$ can be found analytically as closed form expressions.

Define I_k(ω(α), p, λ), with ω(α) = (cos α, sin α) as the value of the Radon transform of the function $\frac{{H}_{\vert k\vert }^{\left(1\right)}\left(\lambda r\left(x\right)\right)}{{H}_{\vert k\vert }^{\left(1\right)}\left(\lambda \right)}\mathrm{exp}\left(\right.\mathrm{i}k\theta \left(x\right)$ at (ω(α), p):

$$\begin{align*}\hfill {I}_{k}\left(\omega \left(\alpha \right),p,\lambda \right)& =\mathcal{R}\left[\frac{{H}_{\vert k\vert }^{\left(1\right)}\left(\lambda r\left(x\right)\right)}{{H}_{\vert k\vert }^{\left(1\right)}\left(\lambda \right)}\mathrm{exp}\left(\mathrm{i}k\theta \left(x\right)\right)\right]\left(\omega \left(\alpha \right),p\right)\hfill \\ \hfill & =\frac{1}{{H}_{\vert k\vert }^{\left(1\right)}\left(\lambda \right)}\underset{\mathbb{R}}{\int }{H}_{\vert k\vert }^{\left(1\right)}\left(\lambda r\left(x\right)\right)\mathrm{exp}\left(\mathrm{i}k\theta \left(x\right)\right)\delta \left(x\cdot \omega \left(\alpha \right)-p\right)\mathrm{d}x.\hfill \end{align*}$$

Compute first

$$\begin{equation*}{H}_{\vert k\vert }^{\left(1\right)}\left(\lambda \right){I}_{k}\left(\omega \left(0\right),p,\lambda \right)=2\underset{0}{\overset{\infty }{\int }}{H}_{\vert k\vert }^{\left(1\right)}\left(\lambda \sqrt{{p}^{2}+{x}_{2}^{2}}\right)\mathrm{cos}\left(k\enspace \mathrm{arcsin}\left(\frac{{x}_{2}}{r}\right)\right)\enspace \mathrm{d}{x}_{2}.\end{equation*}$$

Using formulas 6.941.3 and 6.941.4 in [31], one obtains, for λ > 0,

$$\begin{equation*}{I}_{k}\left(\omega \left(0\right),p,\lambda \right)=2\frac{{\left(-i\right)}^{\vert k\vert }}{\lambda {H}_{\vert k\vert }^{\left(1\right)}\left(\lambda \right)}\mathrm{exp}\left(\mathrm{i}\lambda p\right).\end{equation*}$$

Then, for general α,

$$\begin{equation}{I}_{k}\left(\omega \left(\alpha \right),p,\lambda \right)=2\frac{{\left(-i\right)}^{\vert k\vert }}{\lambda {H}_{\vert k\vert }^{\left(1\right)}\left(\lambda \right)}\mathrm{exp}\left(\mathrm{i}\lambda p\right)\mathrm{exp}\left(\mathrm{i}k\alpha \right).\end{equation} \tag{ 27 }$$

By utilizing the well-known formulas for Hankel functions of a negative argument (see [30], section10.11.7)

$$\begin{equation*}{H}_{\vert k\vert }^{\left(1\right)}\left(-s\right)={\left(-1\right)}^{\vert k\vert -1}{H}_{\vert k\vert }^{\left(2\right)}\left(s\right)={\left(-1\right)}^{\vert k\vert -1}\overline{{H}_{\vert k\vert }^{\left(1\right)}\left(s\right)},\end{equation*}$$

one verifies that (27) is valid for all λ ≠ 0.

In order to turn the above relations into an algorithm, let us assume that M point-like transducers are placed at equispaced point θ_j = jΔθ, j = 0, 1, 2, ‥M, and the measurements are done on an equispaced time grid t_i = iΔt, i = 1, 2, ..., N. We approximate b_k(λ) on a uniform grid of frequencies λ_l = lΔλ, l = −l_Nyq, ..., l_Nyq, where l_NyqΔλ equals to the Nyquist frequency of the discretization in time. This is can be done by discretizing (25) using the trapezoid rule, and applying the fast Fourier transform (FFT) algorithm both in θ and in t. Such an approximation of the outer integral of a periodic function in (25) is quite standard. However, the use of the FFT in t is equivalent to making the boundary condition periodic in time. In order to reduce the error introduced by this effect, the computation of the inner integral should be done with enough zero-padding. On the other hand, computing b_k's using the 2D FFT makes this step asymptotically (and practically) very fast.

Once coefficients b_k(λ_l) are obtained, one can approximate the time derivative of w*(t, r, θ) at t = T:

$$\begin{equation}\frac{\partial {w}^{{\ast}}}{\partial t}\left(T,r,\theta \right)\approx -\sum _{\begin{subarray}{c}l=-{l}_{\mathrm{N}\mathrm{y}\mathrm{q}},\\ l\ne 0\end{subarray}}^{{l}_{\mathrm{N}\mathrm{y}\mathrm{q}}}\sum _{k=-M/2}^{M/2-1}i{\lambda }_{l}{b}_{k}\left({\lambda }_{l}\right){W}_{k}\left(T,{\lambda }_{l},r,\theta \right).\end{equation} \tag{ 28 }$$

Combining this result with (28) and substituting T = 2 we obtain

$$\begin{align}\hfill \frac{\partial }{\partial p}\mathcal{R}\left[w\left(2,x\right)\right]\left(\omega \left(\alpha \right),p+2\right)& \approx \sum _{k=-M/2}^{M/2-1}\left[\sum _{\begin{subarray}{c}l=-{l}_{Nyq},\\ l\ne 0\end{subarray}}^{{l}_{\mathrm{N}\mathrm{y}\mathrm{q}}}\frac{2i{\left(-i\right)}^{\vert k\vert }{b}_{k}\left({\lambda }_{l}\right)}{{H}_{\vert k\vert }^{\left(1\right)}\left({\lambda }_{l}\right)}\mathrm{exp}\left(\mathrm{i}{\lambda }_{l}p\right)\right]\mathrm{exp}\left(\mathrm{i}k\alpha \right)\hfill \\ \hfill & =\sum _{k=-M/2}^{M/2-1}\left[\sum _{l=-{l}_{\mathrm{N}\mathrm{y}\mathrm{q}},}^{{l}_{\mathrm{N}\mathrm{y}\mathrm{q}}}{a}_{k,l}\mathrm{exp}\left(\mathrm{i}{\lambda }_{l}p\right)\right]\mathrm{exp}\left(\mathrm{i}k\alpha \right),\hfill \end{align} \tag{ 29 }$$

with

$$\begin{equation}{a}_{k,l}=\begin{cases}\frac{2i{\left(-i\right)}^{\vert k\vert }{b}_{k}\left({\lambda }_{l}\right)}{{H}_{\vert k\vert }^{\left(1\right)}\left({\lambda }_{l}\right)},\quad \hfill & l\ne 0\hfill \\ 0,\quad \hfill & l=0\hfill \end{cases}.\end{equation} \tag{ 30 }$$

One easily recognizes that formula (29) is a 2D discrete Fourier transform, and thus it can be efficiently implemented using the 2D FFT. Our algorithm then consists in

• (a)  
  Computing coefficients b_k(λ_l) (equation (25) using 2D FFT;
• (b)  
  Computing coefficients a_k,l using (30);
• (c)  
  Computing $\frac{\partial }{\partial p}\mathcal{R}\left[w\left(2,x\right)\right]\left(\omega \left(\alpha \right),p+2\right)$ given by (29) using 2D FFT;
• (d)  
  Computing $\frac{\partial }{\partial p}F\left(\omega ,p\right)\approx 2\frac{\partial }{\partial p}\mathcal{R}\left[w\left(2,x\right)\right]\left(\omega \left(\alpha \right),p+2\right)$ for α ∈ [0, π];
• (e)  
  Computing $\frac{\partial }{\partial p}F\left(\omega ,p\right)\approx -\frac{\partial }{\partial p}F\left(-\omega ,-p\right)$ for α ∈ (π, 2π).

An optional sixth step is to compute F(ω, p) from $\frac{\partial }{\partial p}F\left(\omega ,p\right)$ by numerical anti-differentiation, while enforcing conditions F(ω, 1) = F(ω, −1) = 0 in the least square sense. This step is optional since in the formula (11) (first raw) the Hilbert transform commutes with differentiation, and therefore (11) can be easily modified to reconstruct f(x) from $\frac{\partial }{\partial p}F\left(\omega ,p\right)$.

If one assumes (realistically) that M, N, and l_Nyq are of order of $\mathcal{O}\left(m\right)$, the algorithm requires $\mathcal{O}\left({m}^{2}\mathrm{log}\enspace n\right)$ flops, since the most time consuming operations here are the 2D FFTs.

The algorithm we presented shares important computational steps with the 2D algorithm presented in [32] and with the method applicable to the reduced circular geometry in [18]. We would like to point out that the technique of [32] was developed to solve the inverse source problem with full data; reduced data were not considered. The algorithm of [18] yields a theoretically exact reconstruction from reduced data, but at the expense of a significantly reduced support of source term f(x) as compared to the theoretically visible region. In the present method, the theoretically visible region {x: x ∈ B(0, 1)} and ${x}_{2}{ >}-2\delta \left.\right\}$ is only slightly larger than Ω₀, since δ can be chosen to be much smaller than 1. However, the present technique is not theoretically exact.

In this section performance of our algorithm is demonstrated in numerical simulations. The model acquisition scheme can be seen in figure 2. Measurements where modeled by using 256 samples over time interval [0, 2] (slightly more were used to produce a smooth cut-off after T = 2). There were 1024 detector locations uniformly distributed over the angle interval between 0 and 360°. These were used to model full measurements. When we modeled the reduced measurements, all the values corresponding to angle interval [190°, 350°] were set to zero. We implemented the algorithm described in the previous section in MATLAB. The time required to compute $\frac{\partial }{\partial p}F\left(\omega ,p\right)$ on a 256 × 1024 grid using one thread of the Intel Xeon processor E5620 running at 2.4 GHz was about 1 s. About 85% of that time was used to compute values of the Hankel functions utilized by the algorithm. This suggests that, if needed, computations can be significantly sped up by either pre-computing the Hankel functions, or by implementing the algorithm in C or Fortran that compute these functions significantly faster (cf. the computing times reported in [32]).

Figure 3(b) presents the solution u(t, x) at t = 2 of the problem (5) with full data. Recall that theoretically, in ${\mathbb{R}}^{d}{\backslash}B\left(0,1\right)$ this function coincides with the solution of the forward problem (1) and (2) with the initial condition f(x) shown in part (a) of the figure. Solution of the initial/boundary problem (12) with reduced data is demonstrated in figure 3(c). The difference w(2, x) − u(2, x) (i.e. the error induced by the loss of data) is plotted in the part (d) of the figure. One can see that the singularities in this difference are all located within the lower half space. Moreover, the error within the upper half-space is smooth and quite small.

Zoom In Zoom Out Reset image size

Figure 4(a) presents the Radon projections F(ω, p) computed using our algorithm from the full data. Our algorithm is theoretically exact in this case, and F(ω, p) is indistinguishable from the exact Radon projections. Figure 4(b) demonstrates approximate projections $\tilde {F}\left(\omega \left(\alpha \right),p\right)$ computed from the reduced data. One can see that on the interval α ∈ [0, π] the approximate projections are quite close to F(ω, p). Part (c) of this figure shows projections G(ω(α), p) computed from $\tilde {F}\left(\omega \left(\alpha \right),p\right)$ using symmetry relation (14). One can notice a small jump of values occurring at α = π. However, according to proposition 9 this singularity will not propagate into the reconstructed image.

Zoom In Zoom Out Reset image size

Figures 5 and 6 demonstrate the results of applying the inverse Radon transform to the correct projections F(ω, p) and to our approximation G(ω(α), p). In detail, theoretically exact reconstruction ${\mathcal{R}}^{-1}F$ is shown in figure 5(a). Part (b) shows the image obtained by a naïve reconstruction from reduced data; it was computed using the fast 2D algorithm from [32]. Part (c) of figure 5 presents microlocally accurate reconstruction ${\mathcal{R}}^{-1}G$. Visually it differs little from ${\mathcal{R}}^{-1}F$ exhibited in part (a). The error is about 6% in the relative L_∞ norm and about 3% in the relative L₂ norm. The difference ${\mathcal{R}}^{-1}F-{\mathcal{R}}^{-1}G$ can be seen in figure 5(d), plotted using a much finer color scale. It appears to be smooth, in accordance with the analysis presented in the previous sections.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have also plotted in figure 6 the graphs of the horizontal (along the line x₂ = 0.2) and central vertical cross-sections the phantom f(x) and images shown in figures 5 (b) and (c). One can see that while the naïve reconstruction is incorrect even qualitatively, our technique yields a close fit to the ground truth f(x).

The authors would like to thank L Friedlander and E T Quinto for helpful discussions and useful literature references. The last two authors gratefully acknowledge partial support from the NSF, award NSF/DMS 1814592. M Eller and L Kunyansky are grateful to the University of Lübeck for organizing the workshop 'Modeling, analysis, and approximation theory toward applications in tomography and inverse problems' in June 2016 supported by the Volkswagen Foundation, that initiated their collaboration. Finally, the authors thank the anonymous referees for suggestions on improving this paper.