**1. Introduction**

Photoacoustic tomography (PAT) is a hybrid imaging modality that combines high spatial resolution of ultrasound and high contrast of optical tomography [1–5]. 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 pressure is measured on a detection surface as a function of time. In the case of constant sound speed and when the observation surface exhibits a special geometry (planar, cylindrical, spherical), the initial pressure distribution can be recovered by closed-form inversion formulas; see [6–20] and references therein. While some algorithms take the size and other properties of the detectors into account [21–24], most algorithms assume that the acoustic pressure is known point-wise on a detection surface for all times in the measurement interval. Due to the finite width of the commonly used piezoelectric elements this assumption is only approximately satisfied. One approach to address this issue uses the concept of integrating detectors [25]. Integrating detectors

measure integrals of the acoustic pressure over planes, lines or circles and several inversion methods have been derived in the literature [26–29].

Inspired by the concept of integrating line detectors, a full field detection and reconstruction method that uses projections of the full acoustic field surrounding the sample has been developed in [30,31]. The pressure field within the depth of field is imaged by the use of an optical phase contrast imaging system with a CCD camera. In this way, one obtains 2D projections of the pressure field at a time instant *T*. Reconstructing 2D projections of the initial pressure has been investigated in [32,33]. In [30,31] we go one step further and show that projection data from different directions allow for a full 3D reconstruction of the initial pressure by Radon or Fourier transform techniques. We refer to this approach as full field detection photoacoustic tomography (FFD-PAT). A photograph of the FFD-PAT setup developed in Graz is shown in Figure 1. For practical aspects and a detailed description of the working principle of the FFD-PAT detection system we refer to the original works [30,31]. Existing image reconstruction methods for FFD-PAT data assume a constant speed of sound. However, there are relevant cases when the assumption of constant speed of sound is inaccurate [34,35]. 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 [36–41] for spatially variable speed of sound, that is smooth and bounded from below. Moreover, it is assumed to satisfy the so-called non-trapping 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 these assumptions, it is known that the initial pressure can be stably reconstructed from pressure data given on *S* × [0, *<sup>T</sup>*].

**Figure 1.** FFD-PAT SETUP DEVELOPED IN GRAZ. The sample to be imaged is mounted and rotated on stage A. Full field projections of *p*(· , *T*) are recorded with the CCD camera B using an optical phase contrast technique. For a detailed description of the working principle, see the original works [30,31].

In this paper, for the first time, 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 FFD data. In the second step, we recover the desired initial pressure distribution from the known valued of the pressure at 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. We develop iterative reconstruction methods based on an explicit computation of an adjoint problem. As main theoretical results of this paper, we establish uniqueness and stability of the final time wave inversion problem (see Section 3). In particular, this implies linear convergence for the proposed iterative reconstruction methods. In Section 4, we present numerical results demonstrating that the proposed algorithm is efficient and stable. Moreover, the presented numerical results clearly highlight the importance of taking sound speed variations into account in FFD-PAT image reconstruction.

#### **2. Full Field Detection Photoacoustic Tomography**

In this section, we describe a mathematical model for FFD-PAT including the 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.

## *2.1. Mathematical Model*

In the case of variable sound speed, acoustic wave propagation in PAT is commonly described by the initial value problem [34,38,39,41]

$$p\_{lt}(\mathbf{x},t) - c^2(\mathbf{x})\Delta p(\mathbf{x},t) = 0,\tag{1} \tag{1}$$

$$p(\mathbf{x},0) = f(\mathbf{x}), \tag{2}$$

$$p\_l(\mathbf{x},0) = 0, \quad \qquad \qquad \mathbf{x} \in \mathbb{R}^3. \tag{3}$$

here *c*(**x**) > 0 is the sound speed at location **x** ∈ R<sup>3</sup> and *f* ∈ *<sup>C</sup>*<sup>∞</sup>0 (R<sup>3</sup>) is the initial pressure distribution that encodes the inner structure of the object. It is assumed that the object is contained inside the ball *B*(0, *R*) = **x** ∈ R<sup>3</sup> | **x** < *<sup>R</sup>* , of radius *R* centered at the origin, and that the sound speed is smooth, positive and has the constant value *c*0 outside *B*(0, *<sup>R</sup>*). We denote by *<sup>C</sup>*<sup>∞</sup>0 (*B*(0, *R*)) the set of all smooth functions *f* : R<sup>3</sup> → R that vanish outside *B*(0, *<sup>R</sup>*).

In FFD-PAT, linear projections (integrals along straight lines) of the 3D pressure field *p*(· , *T*) for a fixed time *T* > 0 are recorded. As illustrated in Figure 1, this can be implemented using a special phase contrast method and a CCD-camera that records full field projections of the pressure field [30,31]. The linear projections are collected for rotation angles *θ* ∈ [0, *π*] around the *e*3 = (0, 0, 1) axis and are given by

$$\lg\_R(\theta, s, z) = \int\_{\mathbb{R}} p(s \cos(\theta) - t \sin(\theta), s \sin(\theta) + t \cos(\theta), z, T) \, dt \, \tag{4}$$

for (*<sup>θ</sup>*,*s*, *z*) ∈ *MR* := (*<sup>θ</sup>*,*s*, *z*) ∈ [0, *π*] × R<sup>2</sup> | *s*2 + *z*2 ≥ *<sup>R</sup>*<sup>2</sup> . Here *MR* determines the set of admissible projections and the defining condition *s*2 + *z*2 ≥ *R*<sup>2</sup> means that in practice only pressure integrals over those lines are recorded, which do not intersect the possible support *B*(0, *R*) of the imaged object; compare Figure 2.

**Figure 2.** HORIZONTAL CROSS-SECTION OF THE FFD-PAT SETUP.. Full field projections are measured over lines that do not intersect the ball *B*(0, *<sup>R</sup>*). Consequently, for any plane R<sup>2</sup> × {*z*}, integrals of *p*(· , *T*) over those lines are given that do not intersect the disc *Dz*(0, *R*) := *B*(0, *R*) ∩ (R<sup>2</sup> × {*z*}). For planes with |*z*| ≤ *R*, this yields the exterior problem for the Radon transform.

#### *2.2. Description of the Inverse Problem*

In order to describe the inverse problem of FFD-PAT in a compact way, we introduce some further notation. First, we define the operator

$$\mathbf{W}\_T \colon \mathbb{C}\_0^\infty(B(0, R)) \to \mathbb{C}\_0^\infty(\mathbb{R}^3) \colon f \mapsto p(\cdot, T) \,,\tag{5}$$

where *p* denotes the solution of (1)–(3). The operator **W***T* maps the initial data *f* to the solution (full field) of the wave Equations (1)–(3) at the given measurement time *T* > 0. Second, we define the X-ray transform 

$$\begin{aligned} \mathsf{X} & \mathsf{X} \colon \mathsf{C}\_{0}^{\infty} \left( \mathbb{R}^{3} \right) \to \mathsf{C}\_{0}^{\infty} \left( [0, \pi] \times \mathbb{R}^{2} \right) \\ & \left( \mathsf{X} h \right) (\theta, s, z) := \int\_{\mathbb{R}} h (s \cos(\theta) - t \sin(\theta), s \sin(\theta) + t \cos(\theta), z) \, dt \end{aligned} \tag{6}$$

Note that for any fixed *z* ∈ R, the function (**X***h*)(· , *z*) is the Radon transform of *h*(· , *z*) in the horizontal plane R<sup>2</sup> × *<sup>z</sup>* . Finally, we define the restricted X-ray transform

$$\mathcal{X}\_{\mathbb{R}} \colon \mathbb{C}\_0^{\infty} \left( \mathbb{R}^3 \right) \to \mathbb{C}^{\infty} \left( M\_R \right) : h \mapsto \left( \mathbf{X} h \right) \vert\_{M\_R} \,. \tag{7}$$

For |*z*| ≤ *R*, (**<sup>X</sup>***Rh*)(· , *z*) is the exterior Radon transform of *h*(· , *z*) and consist of integrals over lines not intersecting the disc (*<sup>x</sup>*, *y*) ∈ R<sup>2</sup> | *x*2 + *y*2 < *R*<sup>2</sup> − *z*<sup>2</sup> ; compare Figure 2. Otherwise, (**<sup>X</sup>***Rh*)(· , *z*) coincides with the standard Radon transform of *h*(· , *<sup>z</sup>*).

Using the operator notation introduced above we can write the inverse problem of FFD-PAT in the form

$$\text{Recover } f \text{ from data } \quad g\_R = \mathbf{X}\_R \mathbf{W}\_T f + \text{noise}.\tag{8}$$

Evaluation of **X***R***W***T f* is referred to as the forward problem in FFD-PAT. In this paper, we study the solution of the inverse problem (8).
