*4.1. Discretization*

To implement the steepest descent method (or other gradient type schemes), one has to discretize the final time wave operator **W***T* and its adjoint **<sup>W</sup>**<sup>∗</sup>*T*. For that purpose we solve the forward and adjoint wave Equations (9) and (16) on a cubical grid with side length 2*L* and nodes **x**[**i**] := <sup>−</sup>(*<sup>L</sup>*, *L*, *L*) + **i** 2*L*/*Nx* for **i** = (*<sup>i</sup>*1, *i*2, *i*3) ∈ 0, ... , *Nx* − 1 with the k-space method [52,53], which we briefly recall in Appendix A. Note that the implementation of the k-space method yields a <sup>2</sup>*L*-periodic solution. The parameter *L* is chosen such that *a* + *T* < *L* which implies that inside *B*(0, *<sup>R</sup>*), for *t* ∈ [0, *<sup>T</sup>*], the solution of (16) coincides with its <sup>2</sup>*L*-periodic extension.

We denote by F ⊆ R*Nx*×*Nx*×*Nx* the set of all **f** ∈ R*Nx*×*Nx*×*Nx* with **f**[**i**] = 0 for **x**[**i**] ∈ *B*(0, *<sup>R</sup>*). The discretized versions of **W***T* and its adjoint **<sup>W</sup>**<sup>∗</sup>*T* are defined by W : F → R*Nx*×*Nx*×*Nx* : **f** → (K**f**)(· , *Nt*) and W<sup>T</sup> : R*Nx*×*Nx*×*Nx* → F : **h** → *<sup>χ</sup>R*(K**h**)(· , *Nt*), respectively, where K: R*N*×*N*×*<sup>N</sup>* → R*Nx*×*Nx*×*Nx*×(*Nt*+<sup>1</sup>) denotes the discretized wave propagation defined by the k-space method using the discrete time steps *jT*/*Nt* for 0 ≤ *j* ≤ *Nt*, and *χR* is the discretized indicator function of the ball *B*(0, *<sup>R</sup>*). The linear projection **X** and its left inverse **X** are implemented using the MATLAB build in functions for the Radon and inverse Radon transforms, respectively, with *Nθ* equally spaced projection angles covering [0, *<sup>π</sup>*]. The resulting discretized transforms are denoted by X and X , respectively.

## *4.2. Data Simulation*

For the presented numerical results, we use *R* = 0.4, *T* = 1, *L* = 1.5, *N* = 300, *M* = 300 and collect FFD-data for *Nθ* = 300 projection directions. As initial pressure **f**- we take the sum of three solid spheres, as show in the left picture in Figure 3. The used sound speed **c**- is shown in the right picture in Figure 3, and consist of two Gaussian peaks with opposite signs, added to the constant sound speed **<sup>c</sup>**0[**i**] = 1. The top row in Figure 4 shows the simulated FFD-PAT data **g**- = X W**f**- (left) and the noisy FFD-PAT data **g** := **g**- + noise (right) at projection angle *θ* = 0 for both cases. To generate the noisy FFD data we have added Gaussian white noise with a standard deviation of 10% of the maximal value of **g**-, resulting in a relative *L*2-data error noise/**g**- 85.3%. For the first reconstruction step we apply the inverse Radon transform to the FFD-PAT data, resulting in the approximate 3D final time pressure X **g** W**f**-. The reconstruction of final time pressure from noisy data is shown in the bottom right image in Figure 4. For comparison purpose, the simulated final time pressure **h**- := W**f**- is shown in the bottom left image in Figure 4.

**Figure 3.** PHANTOM AND SOUND SPEED USED FOR THE NUMERICAL SIMULATIONS. (**Left**): Initial pressure **f**- (the phantom to be recovered). (**Right**): Spatially varying sound speed **c**-.

**Figure 4.** FFD-PAT MEASUREMENT DATA AND FINAL TIME PRESSURE. (**Top left**): Simulated FFD-PAT data **g**- for projection direction *θ* = 0. (**Top right**): Noisy FFD-PAT data **g** for projection direction *θ* = 0. (**Bottom left**): Simulated final time pressure **h**- = W**f**-. (**Bottom right**): Recovered final time pressure X **g h**- from noisy FFD-PAT data.

## *4.3. Reconstruction Results*

To recover the initial pressure **f**-, we apply the steepest descent method (Algorithm 1) with *a*- = 0.8 to the point-wise approximation **h** W**f**- that is recovered in step one. Reconstruction results after 5 iterations are shown in Figure 5. The top row shows reconstruction results from simulated data and the middle row shows reconstruction results from noisy data, both using the correct sound speed **c**- for the reconstruction process. One observes, that the reconstruction results from simulated as well as from noisy data are quite accurate. This demonstrates that the proposed algorithm is efficient, accurate and stable with respect to noise. In both cases, the whole reconstruction procedure takes about 90 min on a standard desktop PC.

**Figure 5.** RECONSTRUCTIONS **f**rec USING FIVE ITERATIONS OF THE STEEPEST DESCENT ALGORITHM (**LEFT**) AND THE DIFFERENCES |**<sup>f</sup>**rec − **f**-| TO THE TRUE PHANTOM (**RIGHT**). (**Top**): Using simulated data **g**- and correct sound speed **c**-. (**Middle**): Using noisy data **g** and correct sound speed **c**-. (**Bottom**): Using simulated data **g**- data and wrong sound speed **c**0.

The bottom row in Figure 5 shows reconstruction results using the wrong constant sound speed **c**0 in Algorithm 1 (while data generation still uses the inhomogeneous sound speed **c**-). In this case, the reconstruction error is much larger, showing the relevance of integrating sound speed variations in image reconstruction. Finally, in Figure 6 we show reconstruction results with the previous Fourier method that has been derived in [31] under a constant sound speed assumption. The results with the Fourier method show a large reconstruction error and are very similar to the reconstruction result using the steepest descent method assuming constant sound speed. This demonstrates that the artifacts in both cases are due to the wrong wave propagation model, which further supports the importance of taking sound speed variations into account in FFD-PAT image reconstruction.

**Figure 6.** RECONSTRUCTION **f**rec (**LEFT**) AND THE DIFFERENCE |**<sup>f</sup>**rec − **f**-| (**RIGHT**) USING THE FOURIER ALGORITHM. Note that the Fourier algorithm assumes constant sound speed **c**0 in the image reconstruction, similar to the results shown in the bottom row of Figure 5.

## *4.4. Quantitative Error Analysis*

For all results shown in Figure 5, we stopped the steepest decent method after five iterations, since the reconstructions did not significantly improve by using more iterations. In fact, even after a single iteration, the 2-reconstruction error and the 2-residual error are quite close to their minimal values. This phenomenon is investigated in a more quantitative manner in Figure 7. The images on the left-hand side show the evaluation of the relative *L*2-reconstruction error **<sup>f</sup>***k* − **f**-/**<sup>f</sup>**- in dependence of the iteration index *k*. The images on the right show the relative residual errors W**<sup>f</sup>***k* − **h**-/**<sup>h</sup>**-. The rapid convergence of both error metrics indicates that the finite time inversion problem is well conditioned, which is also suggested by our theoretical analysis presented in Section 3. The relative data errors, residuals and reconstruction errors for all reconstructions shown above are summarized in Table 1.

**Table 1.** RELATIVE ERROR METRICS FOR THE FFD DATA AND RECONSTRUCTIONS. The first column shows the relative 2-norm of the noise added to the data. The second column shows the relative 2-reconstruction error after step one (not applicable to Fourier reconstruction). The third column shows the relative residual error which are minimized in step 2 of the iterative algorithms. The last column shows the relative 2-error of final reconstruction.


**Figure 7.** RELATIVE RECONSTRUCTION ERROR (**LEFT**) AND RELATIVE RESIDUALS (**RIGHT**). (**Top**): Simulated data and correct sound speed. (**Middle**): Noisy data and correct sound speed. (**Bottom**): Simulated data and wrong sound speed (constant and equal to one).
