# Full Field Inversion in Photoacoustic Tomography with Variable Sound Speed
## Full Field Detection Photoacoustic Tomography
### Two-Step Reconstruction  
One possible approach to solve the inverse problem of FFD-PAT is to directly recover \(f\) from data in (8) via iterative methods. Typically, each iteration step requires the evaluation of the forward operator \(\mathbf{X}_{R}\mathbf{W}_{T}\) and its adjoint \((\mathbf{X}_{R}\mathbf{W}_{T})^{*}=\mathbf{W}_{T}^{*}\mathbf{X}_{R}^{*}\). In this paper, we consider a two-step approach where we first invert \(\mathbf{X}_{R}\) via a direct method and then use an iterative method to invert \(\mathbf{W}_{T}\). This avoids repeated and time consuming evaluation of \(\mathbf{X}_{R}\) and its adjoint. The proposed reconstruction method consists of the following two steps.  
* Inverse Radon transform: Assume that projection data \(g_{R}=\mathbf{X}_{R}\mathbf{W}_{T}f\) are given and consider the zero extension \(g\colon[0,\pi]\times\mathbb{R}^{2}\to\mathbb{R}\) by \(g(\theta,x,z)=g_{R}(\theta,x,z)\) for \((\theta,x,z)\in M_{R}\) and \(g(\theta,x,z)=0\) otherwise. We then define an approximation to \(\mathbf{W}_{T}f\) by applying an inversion formula of the Radon transform in planes \(\mathbb{R}^{2}\times\{z\}\). Using the well-known filtered backprojection inversion formula for the Radon transform [42] yields \[\mathbf{W}_{T}f(x,y,z)\simeq\mathbf{X}^{\sharp}g(x,y,z):=\frac{1}{2\pi^{2}} \int_{0}^{\pi}\mathrm{P.V.}\!\!\int_{\mathbb{R}}\frac{(\partial_{\delta}g)\;( \theta,s,z)ds}{(x\cos(\theta)+y\sin(\theta))-s}d\theta\,,\] where P.V. denotes the principal value integral.
* Final time wave inversion: For the second step, assume that an approximation \(h\simeq\mathbf{W}_{T}f\) to the 3D acoustic field at time \(T\) is given. Recovering the initial pressure then yields the final time wave inversion problem \[\text{Recover $f$ from data }\quad h=\mathbf{W}_{T}f+\text{noise}\,.\] (9) To the best of our knowledge, the final time wave inversion problem (9) has not been considered so far and its investigation will be the main theoretical focus of this work.  
For solving the wave inversion problem (9), we propose iterative solution methods that are described in Section 3. As main theoretical results, we derive its uniqueness and stability.  
Another possible approach for solving the first step would be to work with the exterior Radon transform [43; 44; 45]. Instead, we invert the standard Radon transform after setting the missing values of \(\mathbf{X}\mathbf{W}_{T}f\) to zero. We observe numerically that the missing values are indeed close to zero for large enough measurement time \(T\). Although we do not have a rigorous mathematical proof for this claim, it is supported by at least two facts. First, in the case of a non-trapping sound speed, the known decay estimate [46] for the solution \(p\) of (1)-(3) with initial data \(f\) supported in \(B(0,R)\) states  
\[\left|\frac{\partial^{k+|m|}p(\mathbf{x},t)}{\partial_{t}^{k}\partial_{ \mathbf{x}}^{|m|}}\right|\leq\gamma\,e^{-\beta t}\left\|f\right\|_{2}\quad \text{ for }(\mathbf{x},t)\in B(0,R)\times(T,\infty) \tag{10}\]  
for all \((k,m)\in\mathbb{N}^{2}\). Here \(\beta>0\) is a constant only depending on \(c(\mathbf{x})\) and \(T\), and \(\gamma\) is a constant depending on \(B(0,R)\). Second, in the case of constant sound speed, the X-ray transform reduces (1)-(3) to a 2D wave equation with initial data \(\mathbf{X}f\) supported in a disc of radius \(R\). Thus, in the constant sound speed case, \(\mathbf{X}(p(\,\cdot\,,t))\) describes 2D propagation which rapidly decays inside the disc. Formulating and proving precise statements in such a direction, however, is an interesting unsolved problem.