Photoacoustic image reconstruction from full field data in heterogeneous media
Abstract
We consider image reconstruction in fullfield photoacoustic tomography, where 2D projections of the full 3D acoustic pressure distribution at a given time are collected. We discuss existing results on the stability and uniqueness of the resulting image reconstruction problem and review existing reconstruction algorithms. Open challenges are also mentioned. Additionally, we introduce novel onestep reconstruction methods allowing for a variable speed of sound. We apply preconditioned iterative and variational regularization methods to the onestep formulation. Numerical results using the onestep formulation are presented, together with a comparison with the previous twostep approach for fullfield photoacoustic tomography.
Keywords: Photoacoustic Tomography, image reconstruction, full field detection, forwardbackward splitting, onestep reconstruction, heterogeneous medium.
topsep=0em \setitemizetopsep=0em,label=
1 Introduction
Photoacoustic tomography (PAT) is a hybrid imaging modality that beneficially combines the high spatial resolution of ultrasound imaging with the good contrast of optical tomography [1, 2]. In PAT, a semitransparent sample is illuminated by a short laser pulse which induces an acoustic pressure wave depending on the light absorbing structures inside the sample. The induced pressure waves propagate in space, are detected outside of the sample, and measurements are used to recover the photoacoustic (PA) source. The standard approach in PAT is to record timeresolved acoustic signals on a detection surface partially or fully enclosing the investigated object. Plenty of reconstruction methods have been developed for the standard setting. This includes analytic inversion methods [3, 4, 5, 6, 7, 8, 9], time reversal [10, 11, 12], continuous iterative methods [13, 14, 15] and discrete iterative methods [16, 17, 18].
In this paper, we consider full field detection PAT (FFDPAT), which collects different measurement data. In FFDPAT, 2D linear projections of the 3D pressure field at a fixed time instant are captured [19]. This can be implemented by using a special phase contrast method and a CCDcamera that records the full field projections [20]. Similar to the integrating line detector approach. [21, 22, 23], a 3D data set is measured by collecting 2D full field projections from a 1D set of projection directions. Previous reconstructions techniques for FFDPAT where mostly based on a constant sound speed assumption [19, 20]. Recently, in [24] we investigated the case of a variable speed of sound and established the stability and uniqueness in a complete data situation. Moreover, we proposed a twostep reconstruction procedure (including the partial data case), where in the first step the 3D pressure field is recovered using the measured projection data. In the second step, the PA source is recovered from by solving a finite time wave inversion problem.
We study image reconstruction in FFDPAT allowing for a spatially variable speed of sound. In section 2 we formulate the reconstruction problem and recall existing results. In Section 3, we present the new onestep formulation where we recover the PA source directly from the full field data. To stabilize the inversion, we use variational regularization (generalized Tikhonov regularization) including a preconditioning strategy and minimize the Tikhonov functional by proximal forwardbackward splitting. Numerical results for the onestep and the twostep method are presented in Section 4. The paper ends with some conclusions given in Section 5.
2 Full field detection in photoacoustic tomography
In this section, we formulate the image reconstruction problem in FFDPAT and review known uniqueness and stability results as well as existing image reconstruction methods.
2.1 Modeling and problem formulation
We allow a variable sound speed and model the acoustic wave propagation in PAT by the wave equation
(1)  
(2) 
Here is the sound speed at location and is the initial pressure distribution (the PA source) that encodes the inner structure of the sample. We assume that the PA source vanishes outside a bounded volume and that the sound speed takes the constant value outside a bounded volume containing .
To formulate the FFDPAT reconstruction problem we introduce some notations. We denote by the solution of the wave equation (1), (2) with initial data at the given measurement time , and the initialtofinal time wave operator by
(3) 
Here and below is the set of all smooth functions in that have compact support in ; the same notation is used for , where is sufficiently large. We define the Xray transform
(4) 
Then equals the 2D Radon transform of applied in horizontal planes . In particular, can be inverted using any of the existing reconstruction formulas for the 2D Radon transform [25, 26, 27].
In the complete data situation, the challenge in FFDPAT is to reconstruct the PA source from projection data . In practice, one cannot measure the function on the whole projection plane. In particular, at least integrals over lines intersecting the investigated object are missing. Due to practical constraints, projection data might be unknown for even more lines [19, 20]. Thus, in FFDPT we face with the following reconstruction problem.
Problem 1 (Image reconstruction problem in FFDPAT). Let be a given final measurement time and for any , let be the set in the projection plane where the projection data are available. The goal in FFDPAT is to recover the PA source from data , where
(5) 
(with for and for ) is the FFDPAT forward operator.
Problem 1 is closely related to the following finite time wave inversion problem, where in particular the spatial dimensions appear in FFDPAT.
Problem 2 (Finite time wave inversion problem). Let denote the spatial dimension. For given and a region where pressure measurements are made, recover the source from data where is the solution of the dimensional wave equation with initial conditions .
In the remainder of this section we review known results and open problems related to Problems 1 and 2. Especially in the case of a variable sound speed, not so much is known for these image reconstruction problems.
2.2 Review for the constant sound speed case
Problem 1 has first been studied in [19] for a constant speed of sound , where we proposed a noniterative reconstruction method that is outlined in the following. The commutation relation between the 2D Radon transform and the wave equation states . This shows that the FFD data, for any , are given by where solves the 2D wave equation
(6)  
(7) 
Hence the FFDPAT Problem 1 amounts to the solution of the 2D instance of Problem 2 for every . If we have uniqueness, stability or an inversion procedure for Problem 2, we also have corresponding results for Problem 1.
In the general case, uniqueness and stability for Problem 2 are unknown. For some special cases, in [19] we derived an inversion method for the 2D instance based on the reduction to the 1D instance of Problem 2. Below we formulate this procedure for arbitrary dimension. For that purpose, consider the case of full data and denote by
(8) 
the dimensional Radon transform that maps a function defined on and vanishing outside to the integrals of over all hyperplanes (dimensional affine subspaces) of . Let denote the solution of the dimensional wave equation with initial conditions . The commutation relation between the Radon transform and the wave equation shows that satisfies the 1D wave equation with initial conditions . Evaluating the solution of the 1D wave equation at gives . If vanishes outside the ball , then for . For this implies
(9) 
In particular, consists of two separated and translated copies of , where each of the copies can be used to recover the 1D initial source . This implies the following results.
Theorem 3. If , vanishes outside and , then the final time wave inversion Problem 2 is uniquely solvable. Moreover, the initial source can be reconstructed from by the following procedure

Compute the Radon transform .

Invert the 1D wave equation: Choose , set for .

Compute where is any reconstruction method for the Radon transform.
Uniqueness in Theorem 3 refers to the fact that whenever are two distinct PA sources. Moreover, using the twosided stability estimate for the Radon transform, the reconstruction approach of Theorem 3 implies the twosided stability estimate with respect to the norms
(10) 
Here and elsewhere, the notion means that the inequalities hold for some constants . Combining this with the considerations at the beginning of this subsection we obtain the following result for FFDPAT.
Theorem 4. If , and , then the FFDPAT Problem 1 is uniquely solvable via . Moreover the twosided stability estimate holds.
In [19] we also apply the Radon transform approach to the limited data case .
2.3 Review for the variable sound speed case
In the case of variable speed of sound where is not constant, the Xray projections do not satisfy the 2D wave equation. In [24] we therefore proposed a different approach where we first invert the Xray transform and subsequently solve the 3D finite time wave inversion Problem 2.
The following nontrivial fact for Problem 2 has been proven in [24].
Theorem 5. For every , the finite time wave inversion Problem 2 with is uniquely solvable. Moreover, the twosided stability estimate holds.
As a consequence we have the following result for Problem 1 in the variable sound speed case.
Theorem 6. For every , the FFDPAT Problem 1 with is uniquely solvable via . Moreover, the twosided stability estimate holds.
For the limited data case, in [24] we proposed the twostep procedure . In the case that is sufficiently large we observed accurate reconstruction results using the twostep procedure. However, does not exactly invert and therefore the twostep method introduces a systematic (albeit small) error. Moreover, let us mention that in the limiteddata case, no theoretical results on the uniqueness and stability for Problem 1 are known.
3 Preconditioned onestep inversion methods
The application of the inverse Xray transform in the twostep method for limited data introduces a systematic error. Therefore, in this paper we propose a onestep method where we recover directly from data instead of first applying the inverse Xray transform. For implementing the onestep strategy we use iterative and variational reconstruction methods. We include a preconditioned technique accounting for the smoothing of by degree in order to accelerate the iteration.
3.1 Preconditioned onestep Landweber method
Denote by the preconditioning operator defined by , and define the backprojection for . Here is the Fourier transform in second variable. The backprojection operator is the formal adjoint of with the respect to the standard inner product. The composition is the standard filtered backprojection inversion formula for the 2D Radon transform [27] applied with fixed . We write for the formal adjoint of the initialtofinal time wave operator. In [24] we have shown that where is the solution of the backwards wave equation on with and denoting the indicator function of .
Algorithm 7 (Preconditioned onestep Landweber algorithm for FFDPAT).

Initialize: , ,

While (not stop) do:

.

.

for step size

Note that is the adjoint of the forward operator with respect the weighted inner product . Hence the above algorithm undoes the implicit smoothing of and . In particular, assuming the stability estimate , Algorithm 7 is linearly convergent:
(11) 
where . While we expect the stability estimate to be satisfied for important special cases, it is expected not to be satisfied for severe limited data cases or a trapping sound speed. In this case, we expect Problems 1 and 2 to be severely illposed, and that no estimate of the form holds. Variational methods that include an additional regularization term are a reasonable alternative in the illposed and the wellposed case.
3.2 Preconditioned onestep variational regularization
Instead of looking for a theoretically exact solution of the onestep formulation, in variational regularization we minimize the generalized Tikhonov functional
(12) 
where is a regularization term, that additionally stabilizes the iteration and is the regularization parameter. For its numerical solution, we use the forward backward splitting [28], which alternates between a gradient step for the data fitting term and an implicit step with respect to the regularization term .
Algorithm 8 (Preconditioned proximal onestep algorithm of FFDPAT).

Initialize: , ,

While (not stop) do

.

for step size .

The proximal mapping in Algorithm 8 is defined by
(13) 
This implicit treatment of the regularizer allows efficient treatment of nondifferentiable regularizers, where the plain gradient method is not applicable. Also in the differentiable case, the forwardbackward splitting is useful. For example, if , then the Tikhonov functional is illconditioned. The plain gradient method therefore becomes inefficient whereas the forwardbackward splitting treats implicitly and is not affected by the illconditioning. An additional benefit of Algorithm 8 is the explicit smoothing step, which we observed to overall stabilize the iterative procedure.
4 Numerical results
For the presented numerical results, we consider a 2D version of the FFDPAT Problem 1. This arises when we have translational symmetry in the direction meaning that the sound speed and the PA source are independent of the direction. We assume that the PA source is contained in the 2D ball . The 2D full field data are given by where is the solution of the 2D wave equation at time and reduces to the 2D Radon transform. The measurement region is taken as for any .
4.1 Implementation details
All numerical results are obtained using Matlab. In the numerical implementation we solve the 2D wave equation and its adjoint using the kspace method [29, 30, 31] as described in [15]. The Radon transform, its adjoint and inverse we computed with the Matlab build in functions radon and iradon using linear interpolation. Data simulation consist in first computing the solution of the wave equation , then computing the Xray transform and finally multiplying with where with . We discretize the PA source using a Cartesian grid on the domain and numerically compute at time on an Cartesian grid on the domain . The Xray transform is evaluated for directions equidistantly distributed between and .
For image reconstruction we use the proximal gradient version of the onestep method (Algorithm 8). For comparison purpose we also apply the twostep method . For stably inverting we use the proximal gradient method for minimizing the Tikhonov functional . This results in the following twostep variant of Algorithm 8.
Algorithm 9 (Proximal twostep algorithm for FFDPAT).

Compute .

Initialize: , ,

While (not stop) do: for step size .
For the onestep and the twostep algorithm we use the regularizer , the regularization parameter , and a constant step size .
4.2 Simulation and reconstruction results
Figure 1 (left) shows the PA source to be reconstructed. We used a trapping sound speed [11] shown in Figure 1 (right). The simulated data are shown in Figure 2. The top left picture shows the pressure and the top right picture shows its Xray transform to which we have added Gaussian white noise with a standard deviation of of the mean value of , resulting in a relative data error . The FFD data are shown in the bottom right image in Figure 2. The bottom left image shows which is the input for the second step of the twostep method.
Reconstruction results with the onestep (top row) and the twostep (center row) method are shown in Figure 3. The images on the left hand side show the evaluation of the relative reconstruction error in dependence of the iteration index . The images on the right show the reconstructions results of the onestep method (top) and twostep method (center) after 60 iterations. The relative reconstruction error is for the onestep and for the twostep method. Visually as well as in terms of the error, the onestep method shows slightly improved results compared to the twostep method. However, as already observed in [24] the twostep method works remarkably well. The bottom row in Figure 3 shows the reconstruction results assuming a constant sound speed for image reconstruction from the data shown in Figure 2. The constant sound speed has been taken as the average value of the actual sound speed used for data generation. One observes that the reconstruction results using a constant sound speed are inferior, demonstrating that it is crucial to take sound speed variations into account.
5 Conclusion
In this paper, we have studied image reconstruction using full field detection PAT. We reviewed existing result for the inversion problem and introduced new preconditioned onestep reconstruction methods with and without explicit regularization. We presented numerical results for the 2D limited data setting using the preconditioned forward backward onestep splitting Algorithm 8. We compared the results with the previous twostep reconstruction approach, where in the first step the final wave data are approximately recovered by applying the inverse Xray transform to the full field data . For a fair comparison, we used the same regularization term and minimization algorithm for the onestep and the twostep method. As shown in Figure 3 the onestep as well as the twostep method produce accurate results. The onestep method, however, yields slightly fewer visual artefacts and reduces the relative reconstruction error compared to the twostep method.
There are several open problems related to the FFDPAT inversion Problem 1. Uniqueness and stability of reconstruction are known for the full data case, but in the case of limited case neither stability nor uniqueness of reconstruction are known. This is also the case for the related final time wave inversion Problem 2. Our numerical investigations suggest uniqueness as well as stability in the case is taken sufficiently large and the measurement domain is sufficiently large. Theoretically investigating these issues will be subject of future research. Moreover, numerical and experimental investigations in 3D will be performed, especially for the case where the measurement domains in the field data are small and where limited data artifact are expected.
Acknowledgments
M.H. and G.Z acknowledge support of the Austrian Science Fund (FWF), project P 30747N32. The research of L.N. is supported by the National Science Foundation (NSF) Grants DMS 1212125 and DMS 1616904. The work of R.N. has been supported by the FWF, project P 28032.
References
 [1] K. Wang and M. A. Anastasio, “Photoacoustic and thermoacoustic tomography: image formation principles,” in Handbook of Mathematical Methods in Imaging, pp. 781–815, Springer, 2011.
 [2] L. V. Wang, “Multiscale photoacoustic microscopy and computed tomography,” Nature Phot. 3(9), pp. 503–509, 2009.
 [3] L. V. Nguyen, “A family of inversion formulas in thermoacoustic tomography,” Inverse Probl. Imaging 3(4), pp. 649–675, 2009.
 [4] D. Finch, M. Haltmeier, and Rakesh, “Inversion of spherical means and the wave equation in even dimensions,” SIAM J. Appl. Math. 68(2), pp. 392–412, 2007.
 [5] D. Finch, S. K. Patch, and Rakesh, “Determining a function from its mean values over a family of spheres,” SIAM J. Math. Anal. 35(5), pp. 1213–1240, 2004.
 [6] M. Haltmeier, “Universal inversion formulas for recovering a function from spherical means,” SIAM J. Math. Anal. 46(1), pp. 214–232, 2014.
 [7] L. A. Kunyansky, “Explicit inversion formulae for the spherical mean Radon transform,” Inverse Probl. 23(1), pp. 373–383, 2007.
 [8] M. Xu and L. V. Wang, “Universal backprojection algorithm for photoacoustic computed tomography,” Phys. Rev. E 71, 2005.
 [9] P. Burgholzer, G. J. Matt, M. Haltmeier, and G. Paltauf, “Exact and approximate imaging methods for photoacoustic tomography using an arbitrary detection surface,” Phys. Rev. E 75(4), p. 046706, 2007.
 [10] P. Stefanov and G. Uhlmann, “Thermoacoustic tomography with variable sound speed,” Inverse Problems 25(7), pp. 075011, 16, 2009.
 [11] J. Qian, P. Stefanov, G. Uhlmann, and H. Zhao, “An efficient neumann seriesbased algorithm for thermoacoustic and photoacoustic tomography with variable sound speed,” SIAM J. Imaging Sci. 4(3), pp. 850–883, 2011.
 [12] Y. Hristova, P. Kuchment, and L. Nguyen, “Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media,” Inverse Problems 24(5), pp. 055006, 25, 2008.
 [13] S. R. Arridge, M. M. Betcke, B. T. Cox, F. Lucka, and B. E. Treeby, “On the adjoint operator in photoacoustic tomography,” Inverse Problems 32(11), p. 115012 (19pp), 2016.
 [14] Z. Belhachmi, T. Glatz, and O. Scherzer, “A direct method for photoacoustic tomography with inhomogeneous sound speed,” Inverse Problems 32(4), p. 045005, 2016.
 [15] M. Haltmeier and L. V. Nguyen, “Analysis of iterative methods in photoacoustic tomography with variable sound speed,” SIAM J. Imaging Sci. 10(2), pp. 751–781, 2017.
 [16] K. Wang, R. Su, A. A. Oraevsky, and M. A. Anastasio, “Investigation of iterative image reconstruction in threedimensional optoacoustic tomography,” Phys. Med. Biol. 57(17), p. 5399, 2012.
 [17] X. L. DeanBen, A. Buehler, V. Ntziachristos, and D. Razansky, “Accurate modelbased reconstruction algorithm for threedimensional optoacoustic tomography,” IEEE Trans. Med. Imag. 31(10), pp. 1922–1928, 2012.
 [18] C. Huang, K. Wang, L. Nie, and M. A. Wang, L. V.and Anastasio, “Fullwave iterative image reconstruction in photoacoustic tomography with acoustically inhomogeneous media,” IEEE Trans. Med. Imag. 32(6), pp. 1097–1110, 2013.
 [19] R. Nuster, G. Zangerl, M. Haltmeier, and G. Paltauf, “Full field detection in photoacoustic tomography,” Opt. Express 18(6), pp. 6288–6299, 2010.
 [20] R. Nuster, P. Slezak, and G. Paltauf, “High resolution threedimensional photoacoustic tomography with CCDcamera based ultrasound detection,” Biomed. Opt. Express 5(8), pp. 2635–2647, 2014.
 [21] P. Burgholzer, J. BauerMarschallinger, H. Grün, M. Haltmeier, and G. Paltauf, “Temporal backprojection algorithms for photoacoustic tomography with integrating line detectors,” Inverse Probl. 23(6), pp. S65–S80, 2007.
 [22] H. Grün, T. Berer, P. Burgholzer, R. Nuster, and G. Paltauf, “Threedimensional photoacoustic imaging using fiberbased line detectors,” J. Biomed. Optics 15(2), pp. 021306–021306–8, 2010.
 [23] G. Paltauf, R. Nuster, M. Haltmeier, and P. Burgholzer, “Photoacoustic tomography using a MachZehnder interferometer as an acoustic line detector,” Appl. Opt. 46(16), pp. 3352–3358, 2007.
 [24] G. Zangerl, M. Haltmeier, L. V. Nguyen, and R. Nuster, “Full field inversion in photoacoustic tomography with variable sound speed,” 2018. arXiv:1808.00816.
 [25] S. R. Deans, The Radon transform and some of its Applications, John Wiley & Sons, New York, 1983.
 [26] S. Helgason, The Radon Transform, vol. 5 of Progress in Mathematics, Birkhäuser, Boston, second ed., 1999.
 [27] F. Natterer, The Mathematics of Computerized Tomography, vol. 32 of Classics in Applied Mathematics, SIAM, Philadelphia, 2001.
 [28] P. L. Combettes and V. R. Wajs, “Signal recovery by proximal forwardbackward splitting,” Multiscale Model. Sim. 4(4), pp. 1168–1200 (electronic), 2005.
 [29] T. D. Mast, L. P. Souriau, D. D. Liu, M. Tabei, A. I. Nachman, and R. C. Waag, “A kspace method for largescale models of wave propagation in tissue,” IEEE Trans. Acoust. Speech Signal Process. 48(2), pp. 341–354, 2001.
 [30] B. Cox, S. Kara, S. Arridge, and P. Beard, “kspace propagation models for acoustically heterogeneous media: Application to biomedical photoacoustics,” J. Acoust. Soc. Am. 121(6), pp. 3453–3464, 2007.
 [31] M. Tabei, T. D. Mast, and R. C. Waag, “A kspace method for coupled firstorder acoustic propagation equations,” J. Acoust. Soc. Am. 111(1), pp. 53–63, 2002.