1.

Introduction

Photoacoustic computed tomography (PCT) is a rapidly emerging biophotonic imaging modality that combines optical image contrast with ultrasound detection.^1–4 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(\mathbf{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 tissue^5,6 and therefore facilitate a wide-range of diagnostic tasks.^2,7–11

The thermoacoustically induced pressure signals measured in PCT are broadband and acoustic attenuation is frequency-dependent. It has been demonstrated¹² 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.^12–17 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,¹⁸ 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.,¹⁴ 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.

2.

Time Reversal with Heterogeneous Absorption

Let $p(\mathbf{r},t)$ denote the thermoacoustically induced pressure wavefield at location $\mathbf{r}\in {\mathbb{R}}^3$ and time $t\ge 0$. 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.¹⁰ The sought after absorbed energy density $A(\mathbf{r})$ is related to the induced pressure wavefield $p(\mathbf{r},t)$ at $t=0$ as

The measured data $p({\mathbf{r}}_S,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 $\alpha$ can be described by a frequency power law of the form¹⁹

To compensate for acoustic attenuation corresponding to a power law model, Treeby et al.¹⁴ 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:

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(\mathbf{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.²¹ 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.¹⁴ 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 ${\alpha }_0$ and power law exponent $y$.

We modified the original implementation of the k-space model²¹ for use with heterogeneous lossy media. Specifically, two modifications were implemented: (1) The k-space adjustment parameter $\kappa$ 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 ${\alpha }_0$ in Eq. (7) to be a spatially varying quantity ${\alpha }_0(\mathbf{r})$.

Note that although ${\alpha }_0(\mathbf{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 $V_s$ denote the region of support of the most strongly attenuating object component and let ${\alpha }_{0,s}(\mathbf{r})$ and $y_s$ denote the quantities that specify the power law in Eq. (2) for this component. If one specifies $y=y_s$ and ${\alpha }_0(\mathbf{r})={\alpha }_{0,s}(\mathbf{r})$ for $\mathbf{r}\in V_s$ and ${\alpha }_0(\mathbf{r})=0$ otherwise, the k-space time-reversal reconstruction method described above will compensate for acoustic attenuation resulting from the most strongly attenuating component.

3.

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=10\mathrm{mm}$ 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 $3000\mathrm{m}{\mathrm{s}}^{-1}$ and $2000\mathrm{kg}{\mathrm{m}}^{-3}$, and its acoustic attenuation was assumed to be described by the power law ${\alpha }_s(f)={\alpha }_{0,s}{f}^{y_s}$ with ${\alpha }_{0,s}=1\mathrm{dB}{\mathrm{MHz}}^{-y_s}{\mathrm{cm}}^{-1}$ and $y_s=1.5$. When solving Eqs. (3)–(5) the k-space wave solver employed a computational grid of dimension $1\times 512\text{pixels}$ (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 $A_s(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 $A_w(f)$. The frequency-dependent attenuation coefficient was estimated from the simulated measurements as:²²

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 ${\alpha }_0=1\mathrm{dB}{\mathrm{MHz}}^{-y}{\mathrm{cm}}^{-1}$ 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 ${\alpha }_0=1.3\mathrm{dB}{\mathrm{MHz}}^{-y}{\mathrm{cm}}^{-1}$ and $y=0.9$.

4.

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 $1200\mathrm{kg}{\mathrm{m}}^{-3}$ and $3100\mathrm{m}{\mathrm{s}}^{-1}$.

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 $8\mathrm{mJ}/{\mathrm{cm}}^2$ 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.²³ 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 ${\alpha }_0$ and $y$ of acrylic, a transmission experiment was conducted by use of a modified broadband through-transmission technique proposed by He.²² 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 $A_w(f)$ and $A_s(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 ${\alpha }_0=1.3\mathrm{dB}{\mathrm{MHz}}^{-y}{\mathrm{cm}}^{-1}$ and $y=0.9$.

For use in the time-reversal reconstruction code, the 2-D SOS map $c_0(\mathbf{r})$, density map ${\rho }_0(\mathbf{r})$, and attenuation coefficient ${\alpha }_0(\mathbf{r})$ were constructed. The maps $c_0(\mathbf{r})$ and ${\rho }_0(\mathbf{r})$ were assigned the values for acrylic within the annular region occupied by that material and assigned the values $1480\mathrm{m}{\mathrm{s}}^{-1}$ and $1000\mathrm{kg}{\mathrm{m}}^{-3}$ elsewhere. Similarly, the map ${\alpha }_0(\mathbf{r})$ was assigned the value ${\alpha }_0=1.3\mathrm{dB}{\mathrm{MHz}}^{-y}{\mathrm{cm}}^{-1}$ 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\times 500\text{pixels}$ 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 $1520\mathrm{m}{\mathrm{s}}^{-1}$ and constant density value of $1000\mathrm{kg}{\mathrm{m}}^{-3}$ 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.¹⁷ 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.

5.

Summary

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 skull^24,25 may provide a means of estimating $\alpha (\mathbf{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.