Open Access
Add to BibSonomy
Abstract
Almost all reconstruction methods in photoacoustic tomography (PAT) have been developed by assuming that sound propagation is linear, which is valid for ordinary PAT applications but would become inappropriate when the sound amplitude is higher than a certain threshold level. In the current study, we investigate the effect of nonlinear sound propagation on PAT by using a numerical method which utilizes the time-reversal (TR) technique. In the forward stage, the Euler equations are solved to simulate nonlinear sound propagation, and the flow variables (pressure, velocity and density) are recorded by an array of virtual sensors. The recorded data are used to reconstruct the initial fields within the computational domain by using both linear and nonlinear TR techniques. Furthermore, TR results constructed with and without the recorded flow velocity field, which is difficult to measure for practical applications, have also been compared. The current results show that nonlinear reconstructions produce images with superior clarity, resolution and contrast compared to those reconstructed by the linear method, particularly when the recorded velocity field is used in the reconstruction.
© 2017 Optical Society of America
1. Introduction
Photoacoustic tomography for biomedical imaging has attracted growing interests in the last decade since it combines the advantages of both optical and acoustical methods [1–4]. When a tissue sample is irradiated by a short laser pulse, part of the laser energy is absorbed and converted into heat. This local transient heat transfer gives rise to a high pressure distribution due to thermoelastic expansion of the nearby medium. The resulting pressure wave propagates ultrasonically before being picked up by the transducers [5]. Various algorithms have been proposed to reconstruct the initial source location based on linear sound propagation theory [6–9]. More details of these algorithms can be found in the review papers by Li and Wang [10], Wang and Yao [11].
However, to the best knowledge of the authors, previous studies have rarely considered the effect of nonlinear sound propagation on PAT imaging. The nonlinearity could potentially be important for two reasons. First, stronger photoacoustic signals are often desired to improve the contrast, sensor sensitivity and detection specificity of PAT imaging [12,13]. Typically, the initial pressure distribution can be as high as as 2.2 × 104 Pa for a biological tissue [10], and a temperature rise of 1mK can lead to an increase in pressure by approximately 800 Pa [14]. This high resultant sound amplitude leads to a more significant nonlinear effect, and the use of linear wave assumption may lead to inaccurate reconstructions. Second, the sound pressure time signal at a fixed monitoring point is often found to have an N-type waveform [15, 16]. From an analytical aspect, Sun et al. [17] studied the shock wave generation from a decaying photoacoustic sphere irradiated by a laser beam, and identified tooth-like waves that are formed by a nonlinear process. This is similar to the generation of the N-type waveform by a wave with an initial high-amplitude pressure distribution. In this case, for PAT applications, linear methods might not be able to produce detailed and accurate images and nonlinear methods should be used instead.
In this study, the effect of nonlinear sound propagation on reconstructing the initial pressure distribution is demonstrated. The study is conducted by using high fidelity numerical simulations for the nonlinear propagation of photoacoustic sound. For simplicity, the surrounding medium considered is air, which is still relevant in photoacoustic applications [18,19]. In addition, the initial pressure distribution is prescribed directly and the propagation process is studied by the numerical simulations. To highlight the effect of nonlinearity on PAT imagining, TR technique is adopted to reconstruct the initial source location. A benchmark Shepp-Logan phantom image is used to demonstrate the capability of the current method.
2. Methods
The current study considers the nonlinear ultrasonic propagation of sound induced by an initial in-homogeneous pressure distribution in the flow field which has been induced by the photoacoustic effect. We further specify the application in air and assume that the effects of viscosity and heat transfer are negligible. Therefore, the sound propagation is governed by the Euler equations [20]
where γ = 1.4 is the specific heat ratio, p, u and ρ are the flow pressure, velocity and density respectively. Usually, each flow variable can be decomposed into the ambient mean value and its fluctuation about the mean, which are denoted by the subscript (·)0 and the superscript (·)′ respectively, i.e. To use TR to reconstruct the initial pressure distribution, flow variables (p, u and ρ) are recorded by an array of surrounding virtual sensors. Being similar to the TR applications to PAT in the context of linear wave propagation [21–23], the recorded information are introduced to the computational domain to solve Eq. (1) with a time-reversed sequence, i.e. p(t) = pr (T − t), ρ(t) = ρr (T − t) and u(t) = −ur (T − t), where T is the simulation time and the subscript (·)r means the recorded data. It can be easily verified that the time-reversed variables also satisfy the Euler equations, provided that there is no shock wave in the flow field [24]. The presence of shocks can induce irreversible entropy increase, which violates the exact invariance in TR [25,26]. For comparison, we also perform TR computation based on the linearized Euler equations (LEE), which are derived from Eq. (1) under the assumption that the fluctuation variables are of small amplitude Usually, both the pressure p and velocity u are required to obtain an accurate reconstruction. However, it is hard to measure the velocity field in practice. To overcome this issue, we also study both linear and nonlinear TR by using the recorded p as the sole input. In the simulations, the fourth-order low-dissipation and low-dispersion Runge-Kutta scheme [27] is employed for time marching. A tenth-order explicit filter is used to suppress the non-physical spurious waves. Furthermore, absorbing zones are implemented to enforce non-reflective boundary conditions. The high-order weighted essential non-oscillation (WENO) scheme [28] based on local characteristic decomposition and flux splitting is implemented to capture the potential discontinuities when solving the Euler equations. When the linear equations are solved, an advanced compact scheme is employed for spatial discretization [29]. In this study, all the flow variables are non-dimensionalized by using the reference length l* = 1m, reference density ρ* = 1.225kg/m3 is the ambient density, reference speed u* = c0 = 340m/s is the speed of sound, reference time t* = l*/u* and reference pressure . The ambient pressure is assumed as 1 atmosphere pressure (1atm ≈ 105Pa), suggesting that the dimensionless density ρ0 = 1.0, pressure p0 = 1/γ ≈ 0.714.3. Results and discussion
In this section, reconstructions of the standard phantom image created by Shepp and Logan [30] are investigated. The initial pressure distribution has complicated spatial structures, and discontinuities are presented. Two values of (maximum) pressure fluctuation amplitude Ap = 0.1p0 and 0.5p0 are studied. Nonlinearity is less significant to the former case due to its smaller sound amplitude.
For all simulations, the range of the computational domain is x × y = (−2, 2) × (−2, 2) and the discretization sizes are Δx = Δy = 0.005 that is sufficiently small to resolve the acoustic waves, leading to a total number of points of 640, 000. The time step size is Δt = 0.25Δx to maintain the simulation stability. Figure 1 shows the distributions of reconstructed pressure fluctuation p′/p0 for the case with Ap = 0.1p0. Both linear and nonlinear methods can identify the main structures due to the smaller effect of nonlinearity. It should be noted that the levels for different reconstructions shown in the figure are adaptively adjusted to identify the structures clearly. Among the four images, Fig. 1(a) that represents the nonlinear prediction based on both recorded p and u has the best quality and gives the most clearly-resolved small structures. The second best image is Fig. 1(b) that is obtained using the nonlinear method based on the recorded pressure p solely. It can be seen the reconstructed the edges at the interfaces are not as sharp as in Fig. 1(a). The two images reconstructed by the linear method, i.e. Figs. 1(c) and (d), are much blurrier than those by the nonlinear method, with wider edges and poorly-defined small structures.
Fig. 1 The reconstructed pressure fluctuation fields of the Shepp-Logan phantom with the amplitude Ap = 0.1p0. (a) nonlinear result with p and u; (b) nonlinear result with p; (c) linear result with p and u; (d) linear result only p.
Figure 2 shows the instantaneous pressure and density distributions at time t = 1.5 in the forward stage. The pressure field shown in Fig. 2(a) appears to be diffused due to the interference of the nonlinear waves. The sharp-edged wave fronts indicate the presence of shock waves in the flow field caused by the large discontinuities in the initial pressure field. For the density distribution, the profile in the region external to x × y = (−0.5, 0.5) × (−0.5, 0.5) is similar to that of the pressure field. In the middle of the computational domain, a small-amplitude structure with a similar shape to the Shepp-Logan phantom image is present. This interesting structure corresponds to the entropy fluctuation arises from the inhomogeneity of the pressure and the density distributions at the beginning of the simulation. It is motionless because the background mean flow is stationary for this case [31].
Fig. 2 Instantaneous distributions of the pressure (a) and density (b) fields at t = 1.5 for the Ap = 0.1p0 case. The level of the pressure is p ∈ (0.71, 0.75); and the level of the density field is ρ ∈ (0.85, 1.05).
To further demonstrate the nonlinear effect of sound propagation on PAT imaging, the reconstructions for the higher-amplitude case Ap = 0.5p0 by the four methods are shown in Fig. 3. The best result is obtained by using the nonlinear method based on both recorded p and u, as shown in Fig. 3(a). The image reconstructed by the nonlinear method based on recorded pressure p only has lower quality and smaller amplitude since a part of the energy is not introduced to the TR computation (see Fig. 3(b)). Although the main structures can still be identified, the reconstruction with p only shows wider edges at the interfaces and may provide misleading information for tissue diagnosis. The two images produced by the linear methods, shown in Figs. 3(c) and (d), have the worst quality. Nearly all the structures are diffused. This shows that the linear method is not suitable to be used for large-amplitude PAT cases. To quantitatively evaluate the reconstruction quality, we define an error function that
where Ω is the computational domain, pc is the reconstructed results by the TR computations, pi is the initial pressure distribution and Ap is introduced to scale the results by the varying initial amplitude. Figure 4 shows the tendency of ℰ of the reconstruction results by the 4 methods (nonlinear method with/without u; linear method with/without u) with the initial pressure amplitude Ap ranges from 0.001p0 to 0.5p0. For the nonlinear methods, reconstruction error exists and increase with Ap slowly. Especially for the nonlinear method with both p and u as input, the error is maintained at a relatively low level (the red solid line). In addition, below the threshold value , the difference between the linear and nonlinear results are viewed small, suggesting that the nonlinearity maybe unimportant for these cases. For cases with , the value of ℰ increase with Ap drastically, and the reconstruction results only use pressure information are even worse (the blue dashed line). For all cases, the initial pressure field is not exactly reconstructed and ℰ increases with Ap. This is caused by the occurrence of the shocks induced by the discontinuous initial pressure distribution as shown in Fig. 2. In this case, some information is lost during the forward propagation stage and the initial field is therefore not exactly reconstructed [25]. Nevertheless, by solving the nonlinear equations, less information is lost and the result seems to be better, and the main structures in the initial field are fairly well detected.Fig. 3 The reconstructed pressure fields of the Shepp-Logan phantom with the amplitude Ap = 0.5p0. The description for the figures is the same as Fig. 1.
Fig. 4 The tendency of the reconstruction error ℰ with the initial pressure amplitude for the 4 reconstruction methods.
4. Summary
Sound propagation in PAT has always been modelled by linear wave equation. The current study investigates the potential effect of nonlinear sound propagation on PAT image quality. In the current method, forward sound propagation is simulated based on the nonlinear Euler equations. TR computations based on linear and nonlinear wave equations have been conducted. The test results clearly show that nonlinear reconstructions have much higher image quality with better clarity and contrast. It is difficult to measure the velocity field in many practical PAT applications. Therefore, TR simulations based on the recorded pressure only have also been investigated. It is found that with only the recorded pressure, nonlinear reconstructions still have higher image quality than those by the linear method. The current study shows that taking account of nonlinear sound propagation may do help to enhance the quality of PAT imaging. However, it should be noted that only the sound propagation processes is considered in this numerical simulation. Other nonlinear processes such as the bubble oscillation that causes sound generation, are not accounted for.
Acknowledgment
National Natural Science Foundation of China (Grant Nos. 11322222 and 11561130148). The third author received the Newton Advanced Fellowship from the Royal Society (Ref: NA14081) to conduct the research collaboration at the University of Cambridge.
References and links
1. M. H. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum. 77, 041101 (2006). [CrossRef]
2. H. F. Zhang, K. Maslov, G. Stoica, and L. V. Wang, “Functional photoacoustic microscopy for high-resolution and noninvasive in vivo imaging,” Nature Biotechnology 24, 848–851 (2006). [CrossRef] [PubMed]
3. C. Tao and X. J. Liu, “Reconstruction of high quality photoacoustic tomography with a limited-view scanning,” Opt. Express 18, 2760–2766 (2010). [CrossRef] [PubMed]
4. L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science 335, 1458–1462 (2012). [CrossRef] [PubMed]
5. J. J. Yao and L. V. Wang, “Photoacoustic tomography: fundamentals, advances and prospects,” Contrast Media Molecular Imaging 6332–345 (2011). [CrossRef] [PubMed]
6. M. H. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E 71, 016706 (2005). [CrossRef]
7. P. Burgholzer and G. J. Matt, “Exact and approximative imaging methods for photoacoustic tomography using an arbitrary detection surface,” Phys. Rev. E 75, 046706 (2007). [CrossRef]
8. N. N. Cao and A. Nehorai, “Tumor localization using diffuse optical tomography and linearly constrained minimum variance beamforming,” Opt. Express 15, 897–909 (2007). [CrossRef]
9. J. Jose, R. G. H. Willemink, S. Resink, D. Piras, J. C. G. van Hespen, C. H. Slump, W. Steenbergen, T. G. van Leeuwen, and S. Manohar, “Passive element enriched photoacoustic computed tomography (per pact) for simultaneous imaging of acoustic propagation properties and light absorption,” Opt. Express 19, 2093–2104 (2011). [CrossRef] [PubMed]
10. C. H. Li and L. V. Wang, “Photoacoustic tomography and sensing in biomedicine,” Phys. Med. Biol. 54, R59 (2009). [CrossRef] [PubMed]
11. L. V. Wang and J. J. Yao, “A practical guide to photoacoustic tomography in the life sciences,” Nature Methods 13, 627 (2016). [CrossRef] [PubMed]
12. C. W. Wei, M. Lombardo, K. Larson-Smith, I. Pelivanov, C. Perez, J. Xia, T. Matula, D. Pozzo, and M. O’Donnel, “Nonlinear contrast enhancement in photoacoustic molecular imaging with gold nanosphere encapsulated nanoemulsions,” J. Appl. Phys. 104, 033701 (2014).
13. K. Wilson, K. Homan, and S. Emelianov, “Biomedical photoacoustics beyond thermal expansion using triggered nanodroplet vaporization for contrast-enhanced imaging,” Nature Communications 3, 618 (2012). [CrossRef] [PubMed]
14. L. V. Wang and S. Hu, “Photoacoustic Tomography: In Vivo Imaging from Organelles to Organs,” Science 335, 1458–1462 (2012). [CrossRef] [PubMed]
15. L. V. Wang, “Tutorial on photoacoustic microscopy and computed tomography,” IEEE J. Sel. Top. Quantum Electron. 14, 171–179 (2008). [CrossRef]
16. G. J. Diebold, T. Sun, and M. I. Khan, “Photoacoustic monopole radiation in one, two, and three dimensions,” Phys. Rev. Lett. 67, 3384–3387 (1991). [CrossRef] [PubMed]
17. J. M. Sun, B. S. Gerstman, and B. Li, “Bubble dynamics and shock waves generated by laser absorption of a photoacoustic sphere,” J. Appl. Phys. 88, 2352–2362 (2000). [CrossRef]
18. R. G. Stearns and G. S. Kino, “High-frequency photoacoustics in air,” IEEE Transactions on Ultrasonics, Ferro-electrics, and Frequency Control 2, 179–190 (1987). [CrossRef]
19. A. Elia, P. M. Lugara, C. Di Franco, and V. Spagnolo, “Photoacoustic techniques for trace gas sensing based on semiconductor laser sources,” Sensors 9, 9616–9628, (2009). [CrossRef] [PubMed]
20. S. W. Rienstra and A. Hirschberg, “An Introduction to Acoustics” (Eindhoven University of Technology, 2017).
21. B. E. Treeby, E. Z. Zhang, and B. T. Cox, “Photoacoustic tomography in absorbing acoustic media using time reversal,” Inverse Problems 26, 115003 (2010). [CrossRef]
22. Y. Xu and L. V. Wang, “Time reversal and its application to tomography with diffracting sources,” Phys. Rev. Lett. 92, 033902 (2004). [CrossRef] [PubMed]
23. E. Bossy, K. Daoudi, A. Boccara, M. Tanter, J. Aubry, G. Montaldo, and M. Fink, “Time reversal of photoacoustic waves,” Appl. Phys. Lett. 89, 184108 (2006). [CrossRef]
24. K. B. Cunningham, M. F. Hamilton, A. P. Brysev, and L. M. Krutyansky, “Time-reversed sound beams of finite amplitude,” J. Acoustical Soc. Am. 109, 2668–2674 (2001). [CrossRef]
25. M. Tanter, J. L. Thomas, F. Coulouvrat, and M. Fink, “Breaking of time reversal invariance in nonlinear acoustics,” Phys. Rev. E 64, 016602 (2000). [CrossRef]
26. C. Hedberg, “Basics of nonlinear time reversal acoustics,” AIP Conference Proceedings 1106, 164–172 (2009). [CrossRef]
27. F. Q. Hu, M. Y. Hussaini, and J. Manthey, “Low-dissipation and low-dispersion runge-kutta schemes for computational acoustics,” J. Comput. Phys. 124, 177–191 (1996). [CrossRef]
28. G. S. Jiang and C. W. Shu, “Efficient implementation of weighted ENO schemes,” J. Comput.l Phys. 126, 202–228 (1996). [CrossRef]
29. G. Ashcroft and X. Zhang, “Optimized prefactored compact schemes,” J. Comput. Phys. 190, 459–477 (2003). [CrossRef]
30. L. A. Shepp and B. F. Logan, “The Fourier Reconstruction of a Head Section,” IEEE Trans. Nuclear Science 3, 21–43 (1974). [CrossRef]
31. L. S. G. Kovasznay, “Turbulence in supersonic flow,” J. Aeronautical Sci. 20, 657–682 (1953). [CrossRef]
