A Direct Method for Photoacoustic Tomography with Inhomogeneous Sound Speed

The standard approach for photoacoustic imaging with variable speed of sound is time reversal, which consists in solving a well-posed final-boundary value problem for the wave equation backwards in time. This paper investigates the iterative Landweber regularization algorithm, where convergence is guaranteed by standard regularization theory, notably also in cases of trapping sound speed or for short measurement times. We formulate and solve the direct and inverse problem on the whole Euclidean space, what is common in standard photoacoustic imaging, but not for time-reversal algorithms, where the problems are considered on a domain enclosed by the measurement devices. We formulate both the direct and adjoint photoacoustic operator as the solution of an interior and an exterior differential equation which are coupled by transmission conditions. The prior is solved numerically using a Galerkin scheme in space and finite difference discretization in time, while the latter consists in solving a boundary integral equation. We therefore use a BEM-FEM approach for numerical solution of the forward operators. We analyze this method, prove convergence, and provide numerical tests. Moreover, we compare the approach to time reversal.


Introduction
Photoacoustic Imaging (PAI) is a novel imaging technique that allows for three dimensional imaging of small biological or medical specimens with high spatial resolution. It utilizes that an object expands and emits ultrasound waves when it is exposed to a short pulse of electromagnetic radiation (see e.g. [43,40]). The emitted ultrasound is assumed to be proportional to the electromagnetic absorption density which provides detailed anatomical and functional information. PAI aims for visualization of the absorption density by using measurements of the emitted wave outside of the object.
Opposed to the standard PAI problem [24,39,40] we assume a spatially varying speed of sound. The underlying mathematical model of the wave propagation with spatially varying sound speed c(x) is the acoustic wave equation y (x, t) − ∆y(x, t) = 0 in IR n × (0, ∞) , y (x, 0) = 0 in IR n (1) where f denotes the absorption density, which is proportional to the absorbed electromagnetic energy.
The inverse problem of photoacoustics with variable wave speed consists in determining the function f from measurement data m of y on a surface ∂Ω over time (0, T ). That is, it consists in solving the equation With constant sound speed there exists a variety of analytical reconstruction formulas and numerical inversion techniques. We mention [45,44,42,41,46], see also the surveys [24,40,25,38]. In the case of inhomogeneous sound speed, exact reconstruction formulas have been derived in [3]. A numerical method, providing approximations, is time reversal [13,17,16,21], . Thereby the measurement data serve as Dirichlet boundary data, and the initial conditions at final observation time T are assumed to be identically zero. On top of the original algorithm of Fink [13], the approach in [35] suggests a harmonic extension of the boundary data at time T as initialization data. The reconstruction obtained by time-reversal, lets say f T , approximates the true solution f as T increases to infinity. The exact f can be reconstructed if y(·, T ) ≡ 0, which happens, aside from trivial cases, only in odd dimensions and for homogeneous sound speed. In all other cases the results are error-prone. Moreover, this method produces approximations of f only under non-trapping conditions on c (see [21,20]). Stefanov and Uhlmann [35] showed that if diam(Ω) = 2T 0 (where diam(Ω) denotes the diameter of Ω with respect to the Riemannian metric c −2 dx), then an observation time T > T 0 is sufficient for a unique reconstruction of f . However, for stability of the inverse problem one needs longer measurement times and a non-trapping speed of sound condition. In fact, if the measurement time T is larger than a certain threshold, depending on the longest geodesic in the metric induced by c, the algorithm presented in [35] provides a theoretically exact reconstruction in terms of a Neumann series that contains multiple, subsequent time reversal and forward propagation of the data term. A computational realization of this approach has been investigated in [32]. This algorithm serves as a benchmark for our proposed algorithm.
The presented approach for PAI inversion with variable sound speed relies on linear regularization theory [15]. Specifically, we obtain regularized convergence to the minimum norm solution even for short measurement times. Moreover, we obtain a new reconstruction method that does not require an artificial cut-off of the measurement data, nor harmonic extension of the data at the final observation time T .
For the numerical computations, we decouple the wave equation into an interior part (solved by finite element methods), and an exterior part (with homogeneous sound speed) that is rewritten in terms of a boundary integral formulation. We then solve the coupled BEM-FEM system numerically. Note that by this approach we use exact, non-reflecting boundary conditions [1,11] and therefore avoid the necessity of a perfectly matched layer to deal with the cut-off outside the domain of interest. The results are compared to conventional time reversal and the Neumann-series approach.
The paper is organized as follows: In Section 1 we formulate the direct photoacoustic operator L in a suitable choice of function spaces, and derive the adjoint. In Section 2 we give a short overview about Landweber iteration and review some regularization results regarding convergence and convergence rates in view of PAI reconstruction. Moreover, we discuss the relation to the Neumann series approach in [35]. In Section 3, we formulate the transmission problem for the wave equation used for numerical computation of both the direct and the adjoint problem. We state the boundary relations used for taking into account the exterior domain. In addition, we briefly describe the used discretization. Finally, Section 4 provides a comparison of our reconstruction algorithm with the state-of-the-art reconstruction by time reversal and the enhanced time reversal Neumann series method from [35].
All along this paper we use the following notation and abbreviations: Let Ω be a non-empty, open, bounded and connected domain in IR n with C 1 -boundary ∂Ω. The vector n(x), with x ∈ ∂Ω, denotes the outward pointing unit normal vector. We use the following sets Ω + := IR n \Ω, Ω − := Ω and Σ := ∂Ω × (0, T ) .
We use the following Hilbert spaces: • ForΩ = Ω or IR n : -Let H 1 0 (Ω) be the closure of differentiable functions on IR n with compact support inΩ, associated with the inner product -H 1 (Ω) denotes the standard Sobolev space with inner product • L 2 (∂Ω) denotes the standard Hilbert space of square integrable functions on ∂Ω with inner product L 2 (Σ) denotes the standard Hilbert space of square integrable functions on Σ with inner product • The induced norms are denoted by · L 2 (Ω) , · L 2 (Σ) , · H 1 (Ω) and · H 1 0 (Ω) , respectively.
For Ω bounded, on H 1 0 (Ω), the norms · H 1 0 (Ω) and · H 1 (Ω) are equivalent (see [2, Theorem 6.28]): • The trace operator γ Ω : H 1 (IR n ) → L 2 (∂Ω) restricts functions defined on IR n onto ∂Ω, respectively. This operator is the decomposition of the standard trace operator γ : and the restriction operator and thus as a composition of two bounded operators [2, Theorem 5.22] bounded. We abbreviate the norm with Notation 2 The absorption density f and the sound speed c 2 are supposed to satisfy: • Without loss of generality we assume that c ≡ 1 in Ω + .
• The absorption density function f ∈ H 1 0 (Ω) is compactly supported in Ω: For the sake of simplicity of notation we omit space and time arguments of functions whenever this is convenient.

Direct Problem of Wave-Propagation
We analyze the wave operator L mapping the absorption density f onto the solution y of the wave equation (1) restricted to Σ. That is In the following we show that L is bounded. Let us write Computing the derivative of E with respect to t and taking into account (1) gives Consequently which implies that and for every t ∈ (0, T ).
Lemma 1.1 Let y be the solution of (1), then with Proof. First, we note that for arbitrary t ∈ (0, T ), it follows from (8) that: This, together with (9) shows that for all t ∈ (0, T ): In the following we prove boundedness of L: Proof. For given f let y be the solution of (1)). From (6) it follows that the solution y of (10) is in H 1 (IR n ) for every t > 0. Thus from (4) and (10) it follows that which gives the assertion.

Remark 1.3 (Injectivity of L)
In order to obtain injectivity of L, we need T to be sufficiently large. To specify this, we define where dist(x, ∂Ω) is the distance of x to the closest point x ∈ ∂Ω with respect to the Riemannian metric c −2 dx (see also [32]). From [35,Thm. 2] it follows that if In the following we characterize the adjoint of L : H 1 0 (Ω) → L 2 (Σ) on a dense subset of L 2 (Σ). Because we know from elementary functional analysis that L * : That is where ∆ is the Laplace-operator with homogeneous Dirichlet boundary conditions.
In the following we derive the adjoint L * : of the operator L, which is required for the implementation of the Landweber iteration below. where and z := z(h) is the weak solution of Here Proof. For h ∈ C ∞ ((0, T ) × ∂Ω) the existence of a weak solution of (42) is proven in the Appendix. Taking v = y where y denotes the solution of (1) it follows that 2 Landweber Iteration for Solving the Inverse Problem of Photoacoustics The photoacoustic imaging problem rewrites as the solution of the operator equation: If the null-space of L is non-trivial, then iterative regularization algorithms, in general, when m is an element of the range of L reconstruct the minimum norm solution where L † denotes the Moore-Penrose inverse of L (see [31] for a survey on Moore-Penrose inverses). We propose to use the Landweber's iteration for solving (18), because it can be compared with time reversal methods, which are the standard references in this field. More efficient regularization algorithms are at hand [19], but these are less intuitive to be compared with time reversal.
In the following we review properties of the Landweber iteration in an abstract setting (see [15,9]). We use the same notation for the abstract operator and the photoacoustic operator and measurement data m, m δ , respectively, in order to have direct connection.

Abstract Landweber Regularization
Everything that is formulated below is based on the following assumption: Assumption 2.1 Let L : H 1 → H 2 be an operator between Hilbert spaces H 1 and H 2 satisfying ω L 2 ≤ 1 for some ω > 0. Moreover, we assume that data m δ of m is available (typically considered as noisy data), which satisfy Then the Landweber iteration reads as follows: In case δ = 0, that is, if m δ = m, then we write f k instead of f δ k . Let τ > 1 be some fixed constant. The Landweber iteration is only iterated for k = 1, 2, . . . as long as The index, where (22) is violated for the first time is denoted by k * (δ, m δ ).
The following theorem shows that the Landweber iteration converges to the bestapproximating solution: Theorem 2.2 Let m ∈ R(L) (note that the range of L equals the domain of L † ).
In the following we prove properties of the wave-operator L, such that we can apply the general regularization results.
• If T > T 0 , the reconstruction is unique, and f δ k * (δ) converges to the unique solution.
Proof. First, we note that from (11) it follows that ω L 2 ≤ 1. Then, the first two items follow directly from Theorems 2.2.
From the injectivity of L for T > T 0 (see Remark 1.3) it follows that f † = f , which implies unique reconstruction.

Comparison with Time Reversal
We compare our approach with different variants of time reversal. We formally define the time reversal operator:L where z is a solution of The fundamental differences betweenL and L * are thatL is defined for functions with support in Ω and that therefore L * requires a transmission condition in its definition.
Stefanov and Uhlmann [35] modified the time reversal approach in the following sense: Rather than assuming (in most cases unjustified) the initial data z(T ) ≡ 0, they used the harmonic extension of the data term h(s, T ), for s ∈ ∂Ω, as initial datum at T . That is, for −∆φ = 0 in Ω , with φ(·) = m(·, T ) on ∂Ω the modified time-reversal operatorL is defined by the solution of equation They were able to show that under non-trapping conditions and for sufficiently large measurement time T , there exists a compact operator K : H 1 0 (Ω) → H 1 0 (Ω) satisfying K < 1, andL Therefore, the initial condition f can be expanded into the Neumann series By induction, it is easy to see that the m-th iterate can be written as where At this point, we emphasize on the structural similarity between (29) and (21).

Remark 2.4
We emphasize that for time-reversal there is no theory on stopping in case of error-prone data, such as we have available for the Landweber iteration.

Numerical realization of L and L *
We solve (1) and (16) with the same numerical framework. By changing the variable t → T − t in (16), both equations can be rewritten as the transmission problem: where for (1) ρ ≡ 0 and v 0 = f ∈ H 1 0 (Ω) and for (16) we have ρ ≡ h ∈ L 2 (Σ) and v 0 ≡ 0 .
Let v − 0 = v 0 in Ω − , v + 0 = 0 in Ω + , then v ± = v| Ω ± satisfy, respectively: Let G denote the fundamental solution of the standard wave equation with c 2 ≡ 1 in IR n . It is defined in IR n × IR, and its explicit expression with H denoting the Heaviside step function, and δ(·) being the 3D Dirac delta distribution (see, e.g., [14]).
The (retarded) single-and double-layer potentials for (x, t) ∈ Σ are defined by Because we assume that c = 1 outside of Ω and because we assume that ∂Ω is a C 1 boundary, it follows that v + satisfies (see [18]): The transmission conditions imply that Therefore the equation for v = v − on Ω = Ω − can be rewritten as follows: The numerical solution is based on a weak formulation of (35). Integrating over Ω, multiplying by w ∈ H 1 (Ω) and integration by parts of the first line of (35) gives Additionally we introduce the unknown function λ := ∂v ∂n defined on Σ. We are therefore searching for a solution (v, λ) ∈ H 1 (Ω)×H −1/2 (∂Ω) that satisfies for almost all t ∈ (0, T ) the system For detailed analysis of the discretization of equations of that type (36) we refer to [11], since we closely follow their approach. Next, we discuss the time discretization of the integrals V[λ] and K[v]. We consider a uniform time discretization of the interval (0, T ) into N steps of length ∆ t = T /N , and defining the discrete time levels t n = n∆ t . Following [12], by using Lubich's convolution quadrature formula [29], we approximate V[λ] and K[v] at the time steps t n , n = 0, . . . , N by where the coefficients w V n and w K n satisfy where K V (r, s) = K 0 (rs) , K K (r, s) = −sK 1 (rs) ∂r ∂n , and K 0 (·), K 1 (·) are the second kind modified Bessel functions of , respectively. The function γ is given by γ(z) = 3/2 − 2z + 1/2z 2 and is the associated characteristic quotient of the used backward differentiation formula method of order two. For the involved constants we choose L = 2N and β = 1/2N , where is the machine precision (see [29,12] for more details).
In the first equation in (36), the second time derivative is approximated using the second order central difference expression where v n := v(., t n ). The first time derivative occurring in the initial condition is also discretized by central differences, namely From that we can restate the initial conditions as Using this, the explicit Euler discretization of (36) is stated as for all 1 ≤ n ≤ N . In space Ω is triangulated, and we use piecewise quadratic basis functions, supplemented by one cubic function for v. The functions Λ and v| ∂Ω are discretized by the use of piecewise linear functions. With this ansatz a mass-lumped integration scheme can be robustly implemented, which is not the case for purely piecewise quadratic functions. The details and numerical analysis to this scheme can be found in [6].

Numerical experiments and results
Here, we present numerical results for the Landweber reconstructions and compare it with standard time-reversal reconstructions and the Neumann-series-approach (28) -always computed with the same discretization. We concentrate in particular on the cases where classical time reversal techniques have their theoretical and practical drawbacks, that is, so-called trapping speed geometries and short measurement times. In these cases, neither time reversal nor the Neumann-series (28) provide theoretical convergence. As we have shown in Corollary 2.3, the Landweber reconstruction converges to the least squares solution, regardless of the chosen measurement time. In what follows, we also show the practical applicability of our method in such cases.
To create the data sets, we use the discretization of the coupled method introduced in Section 3. The forward computations are performed in a circle of radius 1. For both simulation and reconstructions, the chosen finite element space is that of continuous, piecewise quadratic functions in space. For the spatial discretization of the boundary terms, we take piecewise linear basis functions. Note that for convenience and to optimize the computational effort, we assume ∂Ω to be a circle of Radius R. For the time discretization, wo choose the step size ∆t = h/(15c max ) for the simulation and ∆t r = h r /(14c max ). Here c max is the maximum speed of sound. The time reversal reconstructions are obtained by solving the initial boundary value problem (24) with homogeneous initial values. In the images, when we used the harmonic extension time-reversal (26), the heading is Neumann. The Landweber reconstruction is performed by the scheme described in Corollary 2.3. The choices of c in both the trapping and non-trapping case have been taken from [32].

Non-trapping sound speed
Note that here and below, the speed of sound is smoothed out near the boundaries to satisfy our smoothness assumptions. The time reversal method already gives reliable reconstructions if the measurement time T is sufficiently large. For the test example, a measurement time of T = 4T 0 (see Remark 1.3) is enough to provide a quantitatively reasonable reconstruction by time reversal. This is illustrated in Figure 3. We  therefore are interested in the case where the measurement time is shorter, e.g. T as near as possible at T 0 . In Figure 4 we compare the time reversal reconstruction and the Neumann series approach with our method, using T = 1.2T 0 . The Landweber reconstruction is stopped at a suitable stage of iteration. In practice, the improvement is very large in the first steps, while it needs a lot of iterations to satisfy the discrepancy principle from Theorem 2.2. The main differences are clearly visible to the naked eye: The time reversal reconstruction fails to compute the central point in the image and produces artifacts. These artifacts can be avoided by the use of the harmonic extension as in (24). We here present iterate j = 5 of the Neumann series (28). Also the Landweber reconstruction avoids to amplify the artifacts in the image center. Moreover, it delivers the correct quantitative values of the initial pressure, whereas time reversal underestimates these values. However, the smoothing step naturally included in every Landweber iteration seems to make this approach more stable than the time reversal Neumann series. In fact, at least with our method of numerical wave propagation, we have to stop the Neumann series after 5-10 iterates before numerical errors are amplified too much.

Remark 4.1 (Convergence rates in practice)
The main difference between the convergence rates of the Neumann series and the Landweber iteration lies in the division by c 2 after every backpropagation step, indicated by (15).
In Figure 5 we show reconstructions using the non-trapping speed and measurement time T = 1.2T 0 . The measured error y δ − Lf k L 2 (Σ) keeps decreasing till k = 50. However, the major visible improvements seem to occur within the first 20 iterations. We therefore use the picture of iteration number 20 for our practical comparisons with the other methods.

Trapping sound speed
In the second example, we want to deal with the trapping speed c(x) = 1 + 0.5 sin(−3πx 1 ) cos(3πx), x ∈ B R−ε , In this case, there are geodesics present that do not propagate singularities to the measurement surface within finite time. The Landweber approach is the only one that gives a theoretical convergence result in this case. In practice, we see that conventional time reversal, at least for T = 1.2T 0 , fails to give a detailed reconstruction. The Neumann series approach and the Landweber iteration behave similarly, again with   the advantage of the Neumann series giving faster convergence, whereas the Landweber gives a regularized solution that seems to be more robust against numerical errors and noise (figures 6, 7 and 8). and z m c 2 , w k L 2 (I R n ) + z m , w k H 1 0 (I R n ) = h, w k L 2 (Σ) , for all 0 < t < T , for all k = 1, . . . , m .
2. The following estimates are different to [10, Theorem 1, Section 7.2.2] and thus included -it is essential to consider an additional time integration. As in [10] we use z m (t) =