Translator Disclaimer
7 May 2012 Photoacoustic computed tomography correcting for heterogeneity and attenuation
Author Affiliations +
We report an investigation of image reconstruction in photoacoustic tomography for objects that possess heterogeneous material and acoustic attenuation distributions. When the object contains materials, such as bone and soft-tissue, that are modeled using power law attenuation models with distinct exponents, we demonstrate that the effects of acoustic attenuation due to the most strongly attenuating material can be compensated for if the attenuation of the other less attenuating material(s) are neglected. Experiments with phantom objects are presented to validated our findings.



Photoacoustic computed tomography (PCT) is a rapidly emerging biophotonic imaging modality that combines optical image contrast with ultrasound detection.14 In PCT, the object is illuminated with a pulsed optical wavefield, resulting in the generation of acoustic wavefields via the thermoacoustic effect. From measurements of the acoustic wavefields, a tomographic reconstruction algorithm is employed to obtain an image that depicts the spatially variant absorbed optical energy density distribution A(r) within the object. Because the optical absorption properties of tissue are highly related to its molecular constitution, biomedical applications of PCT can reveal the pathological condition of the tissue5,6 and therefore facilitate a wide-range of diagnostic tasks.2,711

The thermoacoustically induced pressure signals measured in PCT are broadband and acoustic attenuation is frequency-dependent. It has been demonstrated12 that the fidelity of reconstructed images can degrade if acoustic attenuation is not compensated for in the PCT reconstruction algorithm. However, relatively few tomographic reconstruction algorithms are available for such compensation for acoustic attenuation.1217 Moreover, all of the previously investigated methods have assumed that the acoustic attenuation properties of the object are homogeneous. An important biomedical application in which that assumption will be grossly violated is transcranial PCT,18 in which the models of acoustic attenuation in soft-tissue and skull bone have distinct forms.

In this work, we report an investigation of PCT reconstruction of optical absorbers embedded in a heterogeneous, lossy medium. A time-reversal-based reconstruction algorithm proposed by Treeby et al.,14 which was previously demonstrated for media possessing homogeneous acoustic absorption properties, is modified for acoustically heterogeneous and lossy acoustic media obeying a power law attenuation model. As described below, in general the attenuation coefficient component of the power law is permitted to be spatially variant, while the power law exponent is required to be constant. When the object contains materials, such as bone and soft-tissue, that are modeled using power law attenuation models with distinct exponents, we demonstrate that the effects of acoustic attenuation due to the most strongly attenuating material (e.g., bone) can be compensated for if the attenuation due to the other less attenuating material(s) (e.g., soft-tissue) is neglected. Experiments with phantom objects are conducted to corroborate our findings.


Time Reversal with Heterogeneous Absorption

Let p(r,t) denote the thermoacoustically induced pressure wavefield at location rR3 and time t0. Upon propagating out of the object, the pressure wavefields are temporally sampled by use of point-like ultrasonic transducers located at a collection of positions on a closed measurement aperture that surrounds the object. We assume that the effect of the transducer’s electrical impulse response on the measured data has been appropriately deconvolved.10 The sought after absorbed energy density A(r) is related to the induced pressure wavefield p(r,t) at t=0 as

Eq. (1)

where Γ is the dimensionless Grueneisen parameter. From knowledge of the measured data p(rS,t), an estimate of A(r), or equivalently, p(r,t=0), can be obtained by use of a tomographic reconstruction algorithm.

The measured data p(rS,t) contain the effects of acoustic attenuation, which is experienced by the pressure wavefield upon propagating to the ultrasonic transducer locations. For a wide variety of materials, including biological tissues, the acoustic attenuation coefficient α can be described by a frequency power law of the form19

Eq. (2)

where f is the temporal frequency in MHz, α0 is the frequency-independent attenuation coefficient in dBMHzycm1, and y is the power law exponent that is typically in the range of 0 to 3.0.20

To compensate for acoustic attenuation corresponding to a power law model, Treeby et al.14 proposed a method based on time-reversal image reconstruction principles for a lossy acoustic medium. The reconstruction method operates by iteratively solving the following three coupled acoustic equations backward in time:

Eq. (3)


Eq. (4)


Eq. (5)

subject to the initial and boundary conditions:

Eq. (6)

Here, u(r,t) is the acoustic particle velocity, c0(r) is the speed of sound (SOS), ρ(r,t) and ρ0(r) are the acoustic and ambient density, respectively, pm(t) denotes the pressure data measured at transducer location rm, and T denotes the maximum time at which the pm(t) were recorded. The quantities τ(r) and η(r) describe the acoustic absorption and dispersion proportionality coefficients that are defined as

Eq. (7)


Acoustic attenuation is modeled by the absorbing equation of state in Eq. (5), which employes two lossy derivative operators based on the fractional Laplacian to separately accounts for the acoustic absorption and dispersion that are consistent with Eq. (2). Attenuation compensation is achieved by reversing the absorption proportionality coefficient in sign but leaving the equivalent dispersion parameter unchanged during the time-reversal reconstruction. As described in Ref. 14, the solution of these equations yields an image that depicts p(r,t=0) in which the effects of acoustic attenuation have been compensated for. Image reconstruction based on Eqs. (3)(4)–(5) can be accomplished by use of the k-space pseudospectral method.21 Compared with other numerical methods, like the finite-difference time-domain method, k-space pseudospectral methods can significantly improve the computational efficiency by employing the fast Fourier transform (FFT) algorithm to compute the spatial derivatives. The accuracy of the temporal derivatives can also be improved by using k-space adjustments.14 However, to date, this method has only been investigated for reconstructing objects that possess homogeneous acoustic attenuation parameters, i.e. objects whose acoustic attenuation is described by a fixed power law with a constant attenuation coefficient α0 and power law exponent y.

We modified the original implementation of the k-space model21 for use with heterogeneous lossy media. Specifically, two modifications were implemented: (1) The k-space adjustment parameter κ in Eq. (15) in Ref. 14 was removed. This parameter is not required because the equation of state [Eq. (5)] does not involve temporal derivatives, and k-space adjustment is only used to improve the stability and accuracy of the computation of temporal derivatives in the k-space method; and (2) The implementation was modified to permit α0 in Eq. (7) to be a spatially varying quantity α0(r).

Note that although α0(r) can be spatially variant, the power law exponent y is required to be a constant in the k-space time-reversal method. When the object is composed of soft tissues, the assumption of a constant power law exponent is justified. However, when the object contains regions corresponding to distinct power law exponents, which occurs, for example, in the presence of both bone and soft-tissue, the reconstruction method must be modified to avoid image blurring and distortions due to use of a fixed power law exponent. When acoustic attenuation due to a single power law exponent is dominant, e.g., the skull attenuation in transcranial PCT, we propose a simple strategy for circumventing this problem. Namely, the acoustic attenuation effects due to the most strongly attenuating component (e.g., bone) can be compensated for by use of the correct power law parameters, while the less important attenuation effects due to the other component(s) are neglected. Let Vs denote the region of support of the most strongly attenuating object component and let α0,s(r) and ys denote the quantities that specify the power law in Eq. (2) for this component. If one specifies y=ys and α0(r)=α0,s(r) for rVs and α0(r)=0 otherwise, the k-space time-reversal reconstruction method described above will compensate for acoustic attenuation resulting from the most strongly attenuating component.


Computer Simulations

To corroborate the correctness of the modified wave solver code for use with acoustically heterogeneous, lossy media, a computer-simulation study was conducted. The modified wave solver was employed to simulate the propagation of a monopolar pulsed acoustic plane-wave through a one-dimensional heterogeneous lossy medium. The assumed propagation medium consisted of an acoustically absorbing structure of length L=10mm that was embedded in an infinite homogeneous lossless medium with a SOS and density corresponding to water at room temperature. The SOS and density of the absorbing structure were 3000ms1 and 2000kgm3, and its acoustic attenuation was assumed to be described by the power law αs(f)=α0,sfys with α0,s=1dBMHzyscm1 and ys=1.5. When solving Eqs. (3)–(5) the k-space wave solver employed a computational grid of dimension 1×512pixels (51.2 mm), a time step of 1 ns, and a total simulation time of 40 μs.

The pressure wavefield that was propagated through the acoustically inhomogeneous medium was computed as a function of time at the edge of the computational grid. Samples of the magnitude of its one-dimensional Fourier transform As(f) were computed by use of the discrete Fourier transform. The pressure wavefield was also computed at the same location for the case when the acoustic heterogeneity was absent, with the corresponding Fourier magnitude spectrum being denoted as Aw(f). The frequency-dependent attenuation coefficient was estimated from the simulated measurements as:22

Eq. (8)

where f0 is a reference frequency. The estimated u(f) and the corresponding true values given by utrue(f)=α0(fyf0y) were plotted together in Fig. 1(a) and found to be nearly identical. The mean square error between the curves was 0.48.

Fig. 1

A plot of the simulated data function u(f) (blue circles) and the expected function u*(f) (solid line) from 0 to 5 MHz (panel a). This curve corresponds to α0=1dBMHzycm1 and y=1.5. A plot of the experimentally determined data function u(f) (blue circles) and the best fit curve (solid line) (panel b). This curve corresponds to α0=1.3dBMHzycm1 and y=0.9.



Experimental Validation

An experimental PCT study was conducted to demonstrate the use of the modified time-reversal image reconstruction method for use with acoustically heterogeneous, lossy media. A well-characterized phantom, displayed in Fig. 2, represented the to-be-imaged object. The phantom contained 6 optically absorbing structures (pencil leads) embedded in agar. These structures were surrounded by an acrylic cylinder that had inner and outer radii of 7.1 and 7.6 cm, respectively, and a height of 3 cm. The density and SOS of the acrylic were measured and found to be 1200kgm3 and 3100ms1.

Fig. 2

A photograph of the pencil leads held in agar and surrounded by an acrylic cylindrical shell.


The optical absorbers and acrylic cylinder were immersed in water and irradiated by a laser beam from the top. Light from a tunable dye laser (NS, Sirah), pumped by a Q-switched Nd:YAG laser (PRO-350-10, Newport) at 630 nm was used as the illumination source, and the incident laser fluence on the target surface was 8mJ/cm2 with a 10 Hz pulse repetition rate. The photoacoustic (PA) signals were detected by use of a single ultrasound transducer that was scanned along a circular trajectory of radius 9.5 cm. The transducer was cylindrically focused and therefore, the reconstruction problem was treated as a two-dimensional (2-D) one. The photoacoustic signals were recorded with 20 MHz sampling rate at 1000 equally spaced locations on the scanning circle and were amplified by a 50-dB amplifier (5072 PR, Panametrics, Waltham, MA). It has been demonstrated that the 2-D time reversal algorithm can yield accurate reconstructed images if the maximum time of signal recording time T is sufficiently large.23 Therefore, 20,000 temporal samples were acquired at each recording location to ensure that the magnitudes of the PA signals at the cut-off time T were sufficient small (approximately at the noise level).

In the image reconstruction procedure we sought to compensate for acoustic attenuation of the PA signals due to the acrylic cylinder, which represented the dominant acoustic absorber in the object. To determine the absorption parameters α0 and y of acrylic, a transmission experiment was conducted by use of a modified broadband through-transmission technique proposed by He.22 A flat acrylic specimen of thickness 11 mm was employed, whose composition was identical to the acrylic cylinder. The transmitting and receiving transducers employed were both Panametrics V306, having a central frequency of 2.25 MHz with a bandwidth of 70%. From transmission measurements with and without the acrylic specimen present, the corresponding amplitude spectra Aw(f) and As(f) were computed and used to calculate the measured values u(f) in Eq. (8). A nonlinear least squares method was used to fit the measured data to the frequency power law. Figure 1(b) displays the measured values u(f) (blue circles) and the fitted curve u*(f) (solid line). The estimated absorption parameters were found to be α0=1.3dBMHzycm1 and y=0.9.

For use in the time-reversal reconstruction code, the 2-D SOS map c0(r), density map ρ0(r), and attenuation coefficient α0(r) were constructed. The maps c0(r) and ρ0(r) were assigned the values for acrylic within the annular region occupied by that material and assigned the values 1480ms1 and 1000kgm3 elsewhere. Similarly, the map α0(r) was assigned the value α0=1.3dBMHzycm1 within the annular region occupied by the acrylic and was set to zero elsewhere, reflecting that we neglected the relatively weak acoustic attenuation due to the water bath and agar. The power law exponent was set at y=0.9, as determined above.

The measured PA signals were pre-processed by a curvelet denoising technique prior to application of the image reconstruction algorithm. The images were reconstructed on a grid of 500×500pixels of dimension 0.5 mm. To mitigate noise amplification in the reconstructed images, the time-reversed pressure signals were subjected to a low-pass filter specified by a tapered cosine window. The filter cutoff frequency corresponded to the frequency at which the value of average power spectrum of PA signals matched the noise level.

Two additional images were reconstructed to demonstrate the relative importance of compensating for the SOS and density heterogeneities vs. acoustic attenuation. One image was reconstructed by employing a constant SOS value of 1520ms1 and constant density value of 1000kgm3 in the reconstruction algorithm, but properly compensated for the attenuation in the acrylic cylinder. The second image was reconstructed by properly incorporating the spatially variant SOS and density distributions in the reconstruction algorithm, but ignored acoustic attenuation.

The reconstructed images are displayed in Fig. 3. Figure 3(a) displays the reference image corresponding to the case where the acrylic cylinder was absent. Figure 3(b)3(d) displays images of the phantom when the acrylic cylinder was present: Fig. 3(b) displays the image obtained by assuming the constant SOS and mass density values described above but compensating for the acoustic attenuation due to the acrylic; Fig. 3(c) displays the image reconstructed by properly compensating for the spatially variant SOS and density distributions but neglecting acoustic attenuation; the image in Fig. 3(d) was reconstructed by properly compensating for both the spatially variant SOS and density distributions and acoustic attenuation due to the acrylic cylinder.

Fig. 3

The reconstructed image depicting the six optical absorbers in four cases: (a) the acrylic shell was absent during imaging (reference image); (b) the acrylic shell was present and the SOS and density heterogeneities were ignored in the reconstruction method but acoustic attenuation was compensated for; (c) the acrylic shell was present and the SOS and density heterogeneities were compensated for in the reconstruction method but the acoustic attenuation was ignored; and (d) the acrylic shell was present and the SOS and density heterogeneities and acoustic attenuation were compensated for in the reconstruction method.


We found that compensation for SOS and density heterogeneities without attenuation compensation [Fig. 3(c)] yields better spatial resolution than compensation for only attenuation [Fig. 3(b)]. Specifically, the average full-width-at-half-maximum (FWHM) of the optical absorbers in Fig. 3(c) was 25% less than the average FWHM in Fig. 3(b). These results confirm that, for the object studied, the SOS and mass density heterogeneities can influence the reconstructed image more strongly than acoustic attenuation.17 However, as expected, Fig. 3(d) reveals that compensation for both SOS and density heterogeneities along with acoustic attenuation yields the image that possesses the best spatial resolution. This is most noticeable for the optical absorber closest to the acrylic cylinder (far left absorber in the reconstructed images). For that structure, the FWHM in the vertical direction was further reduced by 40% over the FWHM corresponding to Fig. 3(c). This can be explained by the fact that the average acoustic path length through the acrylic cylinder for PA waves generated from the optical absorber closest to the cylinder is longer than for PA waves generated from the other optical absorbers.

Profiles through the centers of the reconstructed images are displayed in Fig. 4. The profiles denoted by solid red, solid green, dotted black, and dashed blue lines correspond to the images in Fig. 3(a)3(d), respectively. The averaged peak magnitude of the six optical absorbers in the reconstructed image with compensation of both SOS and density heterogeneities along with acoustic attenuation (dashed blue line) is 92% of that corresponding to the reference image (solid red line). The averaged peak magnitude in the reconstructed image that compensated only for SOS and density heterogeneities and neglected acoustic attenuation (dotted black line) was 64% of the averaged peak magnitude in the reference image (solid red line), while the reconstructed image that only compensated for acoustic attenuation (solid green) was 57% of that corresponding to the reference image. One notes that in the reconstructed image that only compensates for attenuation (solid green), not only is the peak magnitude underestimated, but the peak positions are also shifted as compared to the reference image. These shifts are larger for the optical absorbers closer to the acrylic cylinder. This demonstrates that, even for relatively simple heterogeneous SOS distributions, using a constant effective SOS value in the reconstruction algorithm can result in image distortions.

Fig. 4

Profiles through the centers of the reconstructed images. The profiles depicted as solid red, solid green, dotted black, and dashed blue lines correspond to the images in Fig. 3(a)3(d), respectively.




We have investigated the use of a time-reversal algorithm for PCT image reconstruction that can compensate for acoustic attenuation in heterogeneous lossy acoustic media. For applications in which acoustic attenuation in a multi-component object is described by frequency power laws having distinct exponents, we demonstrated that the acoustic attenuation due to the most strongly attenuating component can be effectively compensated for. The transmission experiment outlined in this paper to estimate the acoustic attenuation properties of the cylinder is impractical for in-vivo imaging applications. In that case, adjunct imaging data, such as a CT image of the skull24,25 may provide a means of estimating α(r,f), as well as information about the skull geometry, for use with the time-reversal algorithm. We believe that our findings will facilitate the further development of PCT for important applications including transcranial brain imaging.


The authors acknowledge Brad Treeby and Ben Cox for helpful discussions regarding the k-Wave toolbox and Zijian Guo for useful discussions regarding the design of the experiment. Chao Huang, Robert W. Schoonover, and Mark A. Anastasio acknowledge support in part by NIH awards EB010049 and EB009715. Liming Nie and Lihong V. Wang acknowledge support by the NIH awards EB010049, CA134539, EB000712, CA136398, and EB008085. L.V.W. has a financial interest in Microphotoacoustics, Inc. and Endra, Inc., which, however, did not support this work.



L. V. Wang, “Prospects of photoacoustic tomography,” Med. Phys., 35 (12), 5758 –5767 (2008). MPHYA6 0094-2405 Google Scholar


M. XuL. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum., 77 (4), 041101 (2006). RSINAK 0034-6748 Google Scholar


A. A. OraevskyA. A. Karabutov, “Optoacoustic tomography,” Biomedical Photonics Handbook, CRC Press, Boca Raton, FL (2003). Google Scholar


L. V. Wang, “Photoacoustic imaging and spectroscopy,” Photoacoustic Imaging and Spectroscopy, CRC Press, Boca Raton, FL (2009). Google Scholar


W. Joineset al., “Microwave power absorption differences between normal and malignant tissue,” Int. J. Radiat. Oncol. Biol. Phys., 6 (6), 681 –687 (1980). IOBPD3 0360-3016 Google Scholar


W. CheongS. PrahlA. Welch, “A review of the optical properties of biological tissues,” IEEE J. Quant. Electron., 26 (12), 2166 –2185 (1990). IEJQA7 0018-9197 Google Scholar


R. KrugerD. ReineckeG. Kruger, “Thermoacoustic computed tomography-technical considerations,” Med. Phys., 26 (9), 1832 –1837 (1999). MPHYA6 0094-2405 Google Scholar


M. Haltmeieret al., “Thermoacoustic computed tomography with large planar receivers,” Inverse Probl., 20 (5), 1663 –1673 (2004). INPEEY 0266-5611 Google Scholar


B. T. Coxet al., “Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method,” Appl. Opt., 45 (8), 1866 –1875 (2006). APOPAI 0003-6935 Google Scholar


K. Wanget al., “An imaging model incorporating ultrasonic transducer properties for three-dimensional optoacoustic tomography,” IEEE Trans. Med. Imaging, 30 (2), 203 –214 (2011). ITMID4 0278-0062 Google Scholar


Z. XuQ. ZhuL. Wang, “In vivo photoacoustic tomography of mouse cerebral edema induced by cold injury,” J. Biomed. Opt., 16 (6), 066020 (2011). JBOPFO 1083-3668 Google Scholar


P. J. La RivièreJ. ZhangM. A. Anastasio, “Image reconstruction in optoacoustic tomography for dispersive acoustic media,” Opt. Lett., 31 (6), 781 –783 (2006). OPLEDP 0146-9592 Google Scholar


Y. XuL. V. Wang, “Time reversal and its application to tomography with diffracting sources,” Phys. Rev. Lett., 92 (3), 033902 (2004). PRLTAO 0031-9007 Google Scholar


B. E. TreebyE. Z. ZhangB. T. Cox, “Photoacoustic tomography in absorbing acoustic media using time reversal,” Inverse Probl., 26 (11), 115003 (2010). INPEEY 0266-5611 Google Scholar


P. Burgholzeret al., “Compensation of acoustic attenuation for high-resolution photoacoustic imaging with line detectors,” Proc. SPIE, 6437 643724 (2007). PSISDG 0277-786X Google Scholar


D. ModgilM. A. AnastasioP. J. La Rivière, “Photoacoustic image reconstruction in an attenuating medium using singular value decomposition,” in IEEE Nucl. Sci. Conf. R., 4489 –4493 (2008). Google Scholar


X. L. Den-BenD. RazanskyV. Ntziachristos, “The effects of acoustic attenuation in optoacoustic signals,” Phys. Med. Biol., 56 (18), 6129 –6148 (2011). PHMBA7 0031-9155 Google Scholar


L. NieZ. GuoL. V. Wang, “Photoacoustic tomography of monkey brain using virtual point ultrasonic transducers,” J. Biomed. Opt., 16 (7), 076005 (2011). JBOPFO 1083-3668 Google Scholar


T. L. Szabo, “Time domain wave equations for lossy media obeying a frequency power law,” J. Acoust. Soc. Am., 96 (1), 491 –500 (1994). JASMAN 0001-4966 Google Scholar


T. L. Szabo, “Diagnostic ultrasound imaging,” Diagnostic Ultrasound Imaging: Inside Out, Elsevier, Oxford, UK (2004). Google Scholar


B. TreebyB. Cox, “k-wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt., 15 (2), 021314 (2010). JBOPFO 1083-3668 Google Scholar


P. He, “Determination of ultrasonic parameters based on attenuation and dispersion measurements,” Ultrason. Imag., 20 (4), 275 –287 (1998). ULIMD4 0161-7346 Google Scholar


Y. HristovaP. KuchmentL. Nguyen, “Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media,” Inverse Probl., 24 (5), 055006 (2008). INPEEY 0266-5611 Google Scholar


J. F. Aubryet al., “Experimental demonstration of noninvasive transskull adaptive focusing based on prior computed tomography scans,” J. Acoust. Soc. Am., 113 (1), 84 –93 (2003). JASMAN 0001-4966 Google Scholar


S. PichardoV. SinK. Hynynen, “Multi-frequency characterization of the speed of sound and attenuation coefficient for longitudinal transmission of freshly excised human skulls,” Phys. Med. Biol., 56 (1), 219 –250 (2011). PHMBA7 0031-9155 Google Scholar
© 2012 Society of Photo-Optical Instrumentation Engineers (SPIE) 0091-3286/2012/$25.00 © 2012 SPIE
Chao Huang, Liming Nie, Robert W. Schoonover, Lihong V. Wang, and Mark A. Anastasio "Photoacoustic computed tomography correcting for heterogeneity and attenuation," Journal of Biomedical Optics 17(6), 061211 (7 May 2012).
Published: 7 May 2012


Back to Top