Full field inversion in photoacoustic tomography with variable sound speed

Recently, a novel measurement setup has been introduced to photoacoustic tomography, that collects data in the form of projections of the full 3D acoustic pressure distribution at a certain time instant. Existing imaging algorithms for this kind of data assume a constant speed of sound. This assumption is not always met in practice and thus leads to erroneous reconstructions. In this paper, we present a two-step reconstruction method for full field detection photoacoustic tomography that takes variable speed of sound into account. In the first step, by applying the inverse Radon transform, the pressure distribution at the measurement time is reconstructed point-wise from the projection data. In the second step, one solves a final time wave inversion problem where the initial pressure distribution is recovered from the known pressure distribution at the measurement time. For the latter problem, we derive an iterative solution approach, compute the required adjoint operator, and show its uniqueness and stability.


Introduction
Photoacoustic tomography (PAT) is a hybrid imaging modality that combines high spatial resolution of ultrasound and high contrast of optical tomography [3,24,47,46,49]. In PAT, a semitransparent sample is illuminated by a short laser pulse. As a result, parts of the optical energy are absorbed inside the sample. This causes an initial pressure distribution and a subsequent acoustic pressure wave. The pressure wave is detected outside the investigated object and used to recover an image of the interior.
In standard PAT, the induced waves are measured on a surface enclosing the investigated object. In the case of constant sound speed and when the observation surface exhibits a special geometry (planar, cylindrical, spherical), initial pressure distribution can be recovered by closed-form inversion formulas; see [1,8,10,11,17,13,14,16,23,26,28,29,27,48,49] and references therein. All these algorithms assume that the acoustic pressure is known point-wise on a detection surface. Due to finite width of the commonly used piezoelectric elements this assumption is only approximately satisfied. Therefore, the concept of integrating detectors has been invented as an alternative approach to PAT. Integrating detectors measure the integral of acoustic pressure over planes, lines or circles. Closed-form inversion formulas that incorporate integrated pressure data have been developed in [5,12,38,50].
Inspired by the concept of integrating line detectors, a full field detection method that is capable to image the whole acoustic field around an object has been invented in [36,37]. In full field detection PAT (FFD-PAT), a phase contrast method is used to obtain data in the form of 2D projections of the pressure field at a time instant T . If the measurement time T is sufficiently large, then the acoustic pressure has essentially left the object. As shown in [36,37], in the case of constant sound speed, projection data from different directions allow for a full 3D reconstruction of the initial pressure by Radon or Fourier transform techniques.
Existing image reconstruction methods for FFD-PAT assume a constant speed of sound. However, there are relevant cases when the assumption of constant speed of sound is inaccurate [22,25]. For example, it is known that acoustic properties vary within female human breasts. Consequently, for accurate image reconstruction, variable speed of sound has to be incorporated in the wave propagation model. Iterative methods are capable to deal with this assumption. In the case of standard PAT, such methods have been studied in [2,4,15,20,35,42]. Therein the spatially variable speed of sound is assumed to be a smooth function and bounded from below. Moreover, it is assumed to satisfy the so called nontrapping condition, which means that the supremum of the lengths of all geodesics connecting any two points inside the volume enclosed by the measurement surface S is finite. Under this assumption, it is known that the initial pressure can be stably reconstructed from pressure data restricted to S × [0, T ] provided that the measurement time T is sufficiently large.
In this paper, we study image reconstruction in FFD-PAT with a spatially variable speed of sound. We will give a precise mathematical formulation of FFD-PAT and describe the inverse problem we are dealing with (see Section 2). For its solution we propose a two-step process. In the first step, the acoustic pressure at time T is reconstructed pointwise from the full field data. In the second step, we recover the desired initial pressure from the pressure known for fixed measurement time T . The first step can be approximated by inverting the well-known Radon transform. The second step consists in a final time wave inversion problem with spatially varying speed of sound. To the best of our knowledge, the latter has not been addressed in the literature so far. For its solution, we develop iterative reconstruction methods based on an explicit computation of an adjoint problem. As main theoretical results, we establish uniqueness and stability of the final time wave inversion problem. In particular, this implies linear convergence for the proposed iterative reconstruction methods. An object is illuminated by a short pulse of electromagnetic radiation at t = 0. After a sufficiently large time T > 0, the acoustic pressure p( · , T ) has almost left the investigated object and linear projections of p( · , T ) along lines not intersecting B a and perpendicular to the CCD screen are recorded. After that, the object is rotated around the (0, 0, 1) axis and the measurement process is repeated.

Full field detection photoacoustic tomography
In this section, we describe a mathematical model for FFD-PAT including variable sound speed case, and state the inverse problem under consideration. Additionally, we outline the proposed two-step reconstruction procedure and formulate the final time inverse problem.

Mathematical model
In the case of variable sound speed, acoustic wave propagation in PAT is commonly described by the widely accepted model [15,42,20,22] Here c(x) > 0 is the sound speed at location x ∈ R 3 , and f ∈ C ∞ 0 (R 3 ) is the initial pressure distribution that encodes the inner structure of the object. Throughout this text it is assumed that the object is contained inside B a = x ∈ R 3 | x < a , the ball of radius a centered at the origin and that the sound speed is smooth, positive and has the constant value c 0 outside B a .
In FFD-PAT, linear projections (integrals along straight lines) of the 3D pressure field p( · , T ) for a fixed time T > 0 are recorded; compare Figure 1.1. This can be implemented using a special phase contrast method and a CCD-camera that recorders full field projections of the pressure field [36,37]. The projections are collected for rotation angles α ∈ [0, π] around the e 3 = (0, 0, 1) axis and are given by Here M a determines the set of admissible projections, where the defining condition s 2 + z 2 ≥ R 2 means that in practice only pressure integrals over those lines are recorded, which do not intersect the possible support of the imaged object.

Description of the inverse problem
In order to describe the inverse problem of FFD-PAT in a more compact way we introduce some further notation. First, we define the operator where p denotes the solution of (2.1)-(2.3) with initial data f . The operator A maps the initial data f to the solution (full field) of the wave equation (2.1) at the given measurement time T > 0. Second, we define the operator Note that for any fixed z ∈ R, the function (Rh)( · , z) is the Radon transform of h( · , z) in the horizontal plane R 2 × z . Finally, we define the restricted Radon transform where M a is defined in (2.4). For |z| ≤ R, (R a h)( · , z) is the exterior Radon transform of h( · , z) for lines not intersecting (x, y) ∈ R 2 | x 2 + y 2 < R 2 − z 2 ; compare Figure 2.1. Otherwise, (R a h)( · , z) coincides with the standard Radon transform of h( · , z).
Using the operator notation introduced above we can write the inverse problem of FFD-PAT in the form Evaluation of R a Af will be referred to as the forward problem in FFD-PAT. In this paper we study the solution of the inverse problem (2.8).

Two stage reconstruction
One possible approach to solve the inverse problem of FFD-PAT is to directly recover f from data in (2.8) via iterative methods. Typically, each iteration step will require the evaluation of R a A and (R a A) * = A * R * a . In this paper, we consider a two-step approach where we first invert R a via direct method and then use an iterative method to invert A. This avoids repeated and time consuming evaluation of R a and its adjoint.
The proposed two stage reconstruction method consists of the following: Inverse Radon transform: In this first reconstruction step, assume that projection data g a = R a Af are given. Assuming T > 0 to be sufficiently large, we consider the extension g : [0, π]×R 2 → R by g(α, x, z) = g a (α, x, z) for (α, x, z) ∈ M a and g(α, x, z) = 0 otherwise. We then define an approximation to Af by applying an inversion formula of the Radon transform in planes R 2 × {z}. Here we use the well-known filtered backprojection formula (see [31]) which yields where P.V. denotes the principal value integral.
Final time wave inversion: For the second step we assume that an approximation h Af to the 3D acoustic field at time T is given. This yields the final time wave inversion problem To the best of our knowledge, the problem has not been considered so far and its investigation will be the main theoretical focus of this work.
For solving the wave inversion problem (second step), we propose iterative solution methods that are described in detail in Section 3. Additionally, in Section A.1 we derive uniqueness and stability results for (2.9).
Another option for solving the first step would be to work with the exterior Radon transform [40,41,26]. However, we work with the standard Radon transform after replacing the missing values of RAf by zero, since they are approximately zero for large enough T . Theoretically, the smallness is supported by the following two facts. First, in the case of non-trapping sound speed the the known decay estimate for the solution of (2.1) states that the following.
Lemma 2.1 (Decay estimate [45]). Assume that the sound speed c is non-trapping and the initial data f is supported in B a . Then, for any (k, m) ∈ N 2 , the solution p of (2.1)-(2.3) satisfies Here δ > 0 is a constant only depending c and T , and C is a constant depending on the domain B a .
Second, in the case of constant sound speed, the Radon transform R reduces the initial value problem (2.1) to a two dimensional wave equation with initial data Rf which is supported in a disc of radius a. As the sound speed is assumed to be constant outside of B a in the constant sound speed case, RAf rapidly decays in the complement of M a . For non-trapping sound speed we numerically observed the same behaviour. Theoretically investigating this issue, however, is an open problem.

Final time wave inversion problem
In this section we study the final time wave inversion problem (2.9), where the forward operator A : f → p( · , T ) is defined in (2.5). According to standard results for the wave equation [44] the forward operator extends to a bounded linear operator A : L 2 (B a ) → L 2 (R 3 ). Below we establish uniqueness and stability results and derive an iterative reconstruction algorithm using the CG method.
For constant sound speed, recovering the function f from the solution at time T of (2.1) with initial data (0, f ) instead of (f, 0) is equivalent to the the inversion from spherical means at fixed radius. Uniqueness and in inversion method for this problem has been obtained in the classical book of Fritz John [21]. Neither for that case of initial data (f, 0) nor in the variable sound speed case we are not aware of related results.

Uniqueness and stability
The following theorem is the main theoretical result of this paper and states that the final time wave inversion problem (2.9) has a unique solution that stably depends on the right-hand side.
In particular, A : L 2 (B a ) → R(A) has a bounded inverse, where R(A) denotes the range of A. The latter result implies that standard iterative methods for (2.9) converge linearly, similar to the case of standard PAT [15].

Solution by the CG method
To find a solution of (2.9) we use the conjugate gradient (CG) method applied to the normal equation A * Af = A * h. The CG method has proven to be an accurate and fast reconstruction method for the PAT with variable sound speed [15]. Our numerical experiments confirm that the CG method is also efficient for FFD-PAT, where it reads as follows. (S1) Initialize: Using the injectivity and boundedness of A, Algorithm 3.2 generates a series of iterates f k that converge to the unique solution of the inverse source problem (2.9). The stability of (2.9) even implies that the CG method for FFD-PAT converges linearly. More precisely, the sequence (f k ) k∈N generated by Algorithm 3.2 satisfies the estimate (see [7]) where b is defined in (3.1).

The adjoint operator
The CG method requires knowledge of the adjoint operator A * : We show that the adjoint operator is again determined by the solution of a wave equation. More precisely, we have the following result: , consider the time reversed final state problem for the wave equation, and let χ Ba denote the indicator function of B a . Then, , with the final state given by (u( · , T ), u t ( · , T )) = (0, f ). Using the weak formulation (similar to [15]) for the wave equation shows that for every Two times integration by parts, rearranging terms and using the final state conditions (u( · , T ), u t ( · , T )) = (0, f ) yields By taking v as the solution of (2.1)-(2.3) this yields This implies A * g = χ Ba u t ( · , 0) = χ Ba q( · , 0) and completes the proof.
We can reformulate the adjoint operator as follows Then A * g = χ Ba ( · )q( · , T ).

Numerical experiments
In this section we present details on the implementation of CG method (Algorithm 3.2) for FFD-PAT, where the forward operator A and its adjoint A * given by the solution of (2.9) and (3.2), respectively. Numerical experiments are conducted for two variable and two trapping speed of sound models.  The ball B a is considered to be contained inside a discrete 3D cubic region. The side length l of the cube is chosen to be larger than 4R in order to contain the full field p( · , T ).

Discretization and data simulation
To implement the CG method, we have to discretize A and its adjoint A * . For that purpose we solve the forward and adjoint wave equation (2.9) and (3.2) on a cubical grid with nodes For the solution (2.9) and (3.2) we use the k-space method [6,32], which we briefly explain it in Appendix A. We denote by The discretized versions of A and its adjoint A * are defined by Here W N,M : R N ×N ×N → R N ×N ×N ×(M +1) denotes the discretized wave propagation defined by the k-space method using the discrete time steps jT /M for 0 ≤ j ≤ M .
The discrete and inverse Radon transforms R and R are computed by the standard Matlab implementation of the Radon transform and its inverse.

Sound speed models and phantom
In our numerical setup, we use four different variable sound speed models (A, B, C and D) which are shown in Figure 4.2. All variable sound speed models deviate within 30% from the background sound speed c 0 = 1. The two speed of sound models A and B (shown in the upper row) are acoustically non-trapping whereas the speed of sound models C and D (shown the bottom row) are trapping. The speed of sound models A and B have the form where a sum of Gaussian pulses centered at y j added to the background sound speed. The first non-trapping speed of sound model A consists of several pulses with small width, whereas the second model B is a single pulse with a very large width.
In the trapping case C we consider a cavity in the middle of region B a , which is the difference of the constant speed of sound with a Gaussian pulse. The sound speed D of sound is of the type Since in cases C and D, c is radially symmetric circles concentric to the origin are closed geodesics that never leave B a . Therefore, this sound speed cases serve as an test case for a trapping speed of sound.
We assume the sum of three solid spheres as initial pressure f , which is depicted in Figure 4     Once point-wise data p( · , T ) are approximated on R 3 we use the CG method outlined in section 3 to numerically compute the initial pressure f of system (2.1)-(2.3). Figure 4.7 shows slice images of the initial pressure corresponding to the different sound speed cases after four iterations of the CG method.

Application of the CG method
As expected reconstructions are better for the non-trapping speed of sound models. Figure 4.8 shows the output of the algorithm after one iteration. Finally, Figure  4.9 shows the reconstruction result, when a wrong speed of sound (namely constant value 1) is used. From this we can clearly see that not accounting for variable speed can introduce a significant error.

Conclusion
In this paper, we described a FFD-PAT method, where projection data of acoustic pressure are measured. For the first time, we consider this method for variable speed of sound. We developed a two-step reconstruction procedure that computes 3D acoustic pressure data point-wise in a first step and then uses them as input for an iterative algorithm in a second step. We prove uniqueness and stability estimates for the second step. Furthermore, in upcoming works we will also study   Proposition A.2. There is a constant C such that where K is a pseudo-differential operator of order at most −1.
Combining with a similar argument for A * − A − , we obtain that the principal symbol of A * A is α(x, ξ) = 1 2 . That is, A * A = 1 2 I + K, where K is a pseudodifferential operator of order at most −1. We conclude that f 2 L 2 ≤ 2( Af 2 L 2 + Kf 2 L 2 ) .
We are now ready to prove Theorem 3.1.
Proof of Theorem 3.1. Let us recall from Proposition A.2 where K is a pseudo-differential operator of order at most −1. Since K is compact and A is injective, applying [43, Theorem V.3.1], we obtain for some constant C ∈ (0, ∞). This finishes our proof.

A.2 k-space method
We briefly describe the k-space method for the 3D wave equation (2.1)-(2.3) as we use it for the numerical computation of A and A * . The k-space method is an attractive alternative to standard methods using finite differences, finite elements or pseudospectral methods, since it does not suffer from numerical dispersion [6,32]. w(x, t).
In the k-space method we use the time stepping formula where F x and F −1 ξ denote the Fourier transform and its inverse with respect to space and frequency variables x and ξ and h t is the time step size. This equivalent formulation motivates the following algorithm for numerically solving the wave equation.