Photoacoustic image reconstruction from full field data in heterogeneous media
We consider image reconstruction in full-field 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 one-step reconstruction methods allowing for a variable speed of sound. We apply preconditioned iterative and variational regularization methods to the one-step formulation. Numerical results using the one-step formulation are presented, together with a comparison with the previous two-step approach for full-field photoacoustic tomography.
Keywords: Photoacoustic Tomography, image reconstruction, full field detection, forward-backward splitting, one-step reconstruction, heterogeneous medium.
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 semi-transparent 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 time-resolved 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 (FFD-PAT), which collects different measurement data. In FFD-PAT, 2D linear projections of the 3D pressure field at a fixed time instant are captured . This can be implemented by using a special phase contrast method and a CCD-camera that records the full field projections . 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 FFD-PAT where mostly based on a constant sound speed assumption [19, 20]. Recently, in  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 two-step 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 FFD-PAT 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 one-step 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 forward-backward splitting. Numerical results for the one-step and the two-step 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 FFD-PAT 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
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 FFD-PAT 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 initial-to-final time wave operator by
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 X-ray transform
In the complete data situation, the challenge in FFD-PAT 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 FFD-PT we face with the following reconstruction problem.
Problem 1 (Image reconstruction problem in FFD-PAT). 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 FFD-PAT is to recover the PA source from data , where
(with for and for ) is the FFD-PAT forward operator.
Problem 1 is closely related to the following finite time wave inversion problem, where in particular the spatial dimensions appear in FFD-PAT.
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  for a constant speed of sound , where we proposed a non-iterative 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
Hence the FFD-PAT 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  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
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
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 two-sided stability estimate for the Radon transform, the reconstruction approach of Theorem 3 implies the two-sided stability estimate with respect to the -norms
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 FFD-PAT.
Theorem 4. If , and , then the FFD-PAT Problem 1 is uniquely solvable via . Moreover the two-sided stability estimate holds.
In  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 X-ray projections do not satisfy the 2D wave equation. In  we therefore proposed a different approach where we first invert the X-ray transform and subsequently solve the 3D finite time wave inversion Problem 2.
The following non-trivial fact for Problem 2 has been proven in .
Theorem 5. For every , the finite time wave inversion Problem 2 with is uniquely solvable. Moreover, the two-sided 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 FFD-PAT Problem 1 with is uniquely solvable via . Moreover, the two-sided stability estimate holds.
For the limited data case, in  we proposed the two-step procedure . In the case that is sufficiently large we observed accurate reconstruction results using the two-step procedure. However, does not exactly invert and therefore the two-step method introduces a systematic (albeit small) error. Moreover, let us mention that in the limited-data case, no theoretical results on the uniqueness and stability for Problem 1 are known.
3 Preconditioned one-step inversion methods
The application of the inverse X-ray transform in the two-step method for limited data introduces a systematic error. Therefore, in this paper we propose a one-step method where we recover directly from data instead of first applying the inverse X-ray transform. For implementing the one-step 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 one-step 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  applied with fixed . We write for the formal -adjoint of the initial-to-final time wave operator. In  we have shown that where is the solution of the backwards wave equation on with and denoting the indicator function of .
Algorithm 7 (Preconditioned one-step Landweber algorithm for FFD-PAT).
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:
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 ill-posed, and that no estimate of the form holds. Variational methods that include an additional regularization term are a reasonable alternative in the ill-posed and the well-posed case.
3.2 Preconditioned one-step variational regularization
Instead of looking for a theoretically exact solution of the one-step formulation, in variational regularization we minimize the generalized Tikhonov functional
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 , 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 one-step algorithm of FFD-PAT).
Initialize: , ,
While (not stop) do
for step size .
The proximal mapping in Algorithm 8 is defined by
This implicit treatment of the regularizer allows efficient treatment of non-differentiable regularizers, where the plain gradient method is not applicable. Also in the differentiable case, the forward-backward splitting is useful. For example, if , then the Tikhonov functional is ill-conditioned. The plain gradient method therefore becomes inefficient whereas the forward-backward splitting treats implicitly and is not affected by the ill-conditioning. 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 FFD-PAT 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 k-space method [29, 30, 31] as described in . 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 X-ray 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 X-ray transform is evaluated for directions equidistantly distributed between and .
For image reconstruction we use the proximal gradient version of the one-step method (Algorithm 8). For comparison purpose we also apply the two-step method . For stably inverting we use the proximal gradient method for minimizing the Tikhonov functional . This results in the following two-step variant of Algorithm 8.
Algorithm 9 (Proximal two-step algorithm for FFD-PAT).
Initialize: , ,
While (not stop) do: for step size .
For the one-step and the two-step 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  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 X-ray 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 two-step method.
Reconstruction results with the one-step (top row) and the two-step (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 one-step method (top) and two-step method (center) after 60 iterations. The relative -reconstruction error is for the one-step and for the two-step method. Visually as well as in terms of the -error, the one-step method shows slightly improved results compared to the two-step method. However, as already observed in  the two-step 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.
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 one-step reconstruction methods with and without explicit regularization. We presented numerical results for the 2D limited data setting using the preconditioned forward backward one-step splitting Algorithm 8. We compared the results with the previous two-step reconstruction approach, where in the first step the final wave data are approximately recovered by applying the inverse X-ray transform to the full field data . For a fair comparison, we used the same regularization term and minimization algorithm for the one-step and the two-step method. As shown in Figure 3 the one-step as well as the two-step method produce accurate results. The one-step method, however, yields slightly fewer visual artefacts and reduces the relative -reconstruction error compared to the two-step method.
There are several open problems related to the FFD-PAT 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.
M.H. and G.Z acknowledge support of the Austrian Science Fund (FWF), project P 30747-N32. 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.
-  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.
-  L. V. Wang, “Multiscale photoacoustic microscopy and computed tomography,” Nature Phot. 3(9), pp. 503–509, 2009.
-  L. V. Nguyen, “A family of inversion formulas in thermoacoustic tomography,” Inverse Probl. Imaging 3(4), pp. 649–675, 2009.
-  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.
-  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.
-  M. Haltmeier, “Universal inversion formulas for recovering a function from spherical means,” SIAM J. Math. Anal. 46(1), pp. 214–232, 2014.
-  L. A. Kunyansky, “Explicit inversion formulae for the spherical mean Radon transform,” Inverse Probl. 23(1), pp. 373–383, 2007.
-  M. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E 71, 2005.
-  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.
-  P. Stefanov and G. Uhlmann, “Thermoacoustic tomography with variable sound speed,” Inverse Problems 25(7), pp. 075011, 16, 2009.
-  J. Qian, P. Stefanov, G. Uhlmann, and H. Zhao, “An efficient neumann series-based algorithm for thermoacoustic and photoacoustic tomography with variable sound speed,” SIAM J. Imaging Sci. 4(3), pp. 850–883, 2011.
-  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.
-  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.
-  Z. Belhachmi, T. Glatz, and O. Scherzer, “A direct method for photoacoustic tomography with inhomogeneous sound speed,” Inverse Problems 32(4), p. 045005, 2016.
-  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.
-  K. Wang, R. Su, A. A. Oraevsky, and M. A. Anastasio, “Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography,” Phys. Med. Biol. 57(17), p. 5399, 2012.
-  X. L. Dean-Ben, A. Buehler, V. Ntziachristos, and D. Razansky, “Accurate model-based reconstruction algorithm for three-dimensional optoacoustic tomography,” IEEE Trans. Med. Imag. 31(10), pp. 1922–1928, 2012.
-  C. Huang, K. Wang, L. Nie, and M. A. Wang, L. V.and Anastasio, “Full-wave iterative image reconstruction in photoacoustic tomography with acoustically inhomogeneous media,” IEEE Trans. Med. Imag. 32(6), pp. 1097–1110, 2013.
-  R. Nuster, G. Zangerl, M. Haltmeier, and G. Paltauf, “Full field detection in photoacoustic tomography,” Opt. Express 18(6), pp. 6288–6299, 2010.
-  R. Nuster, P. Slezak, and G. Paltauf, “High resolution three-dimensional photoacoustic tomography with CCD-camera based ultrasound detection,” Biomed. Opt. Express 5(8), pp. 2635–2647, 2014.
-  P. Burgholzer, J. Bauer-Marschallinger, H. Grün, M. Haltmeier, and G. Paltauf, “Temporal back-projection algorithms for photoacoustic tomography with integrating line detectors,” Inverse Probl. 23(6), pp. S65–S80, 2007.
-  H. Grün, T. Berer, P. Burgholzer, R. Nuster, and G. Paltauf, “Three-dimensional photoacoustic imaging using fiber-based line detectors,” J. Biomed. Optics 15(2), pp. 021306–021306–8, 2010.
-  G. Paltauf, R. Nuster, M. Haltmeier, and P. Burgholzer, “Photoacoustic tomography using a Mach-Zehnder interferometer as an acoustic line detector,” Appl. Opt. 46(16), pp. 3352–3358, 2007.
-  G. Zangerl, M. Haltmeier, L. V. Nguyen, and R. Nuster, “Full field inversion in photoacoustic tomography with variable sound speed,” 2018. arXiv:1808.00816.
-  S. R. Deans, The Radon transform and some of its Applications, John Wiley & Sons, New York, 1983.
-  S. Helgason, The Radon Transform, vol. 5 of Progress in Mathematics, Birkhäuser, Boston, second ed., 1999.
-  F. Natterer, The Mathematics of Computerized Tomography, vol. 32 of Classics in Applied Mathematics, SIAM, Philadelphia, 2001.
-  P. L. Combettes and V. R. Wajs, “Signal recovery by proximal forward-backward splitting,” Multiscale Model. Sim. 4(4), pp. 1168–1200 (electronic), 2005.
-  T. D. Mast, L. P. Souriau, D. D. Liu, M. Tabei, A. I. Nachman, and R. C. Waag, “A k-space method for large-scale models of wave propagation in tissue,” IEEE Trans. Acoust. Speech Signal Process. 48(2), pp. 341–354, 2001.
-  B. Cox, S. Kara, S. Arridge, and P. Beard, “k-space propagation models for acoustically heterogeneous media: Application to biomedical photoacoustics,” J. Acoust. Soc. Am. 121(6), pp. 3453–3464, 2007.
-  M. Tabei, T. D. Mast, and R. C. Waag, “A k-space method for coupled first-order acoustic propagation equations,” J. Acoust. Soc. Am. 111(1), pp. 53–63, 2002.