1.

Introduction

Photoacoustic tomography (PAT) is a noninvasive biomedical imaging technique, based on the photoacoustic effect, which combines the optical absorption contrast of tissue with the high spatial resolution of ultrasound imaging techniques.^1–4 In PAT, the target is illuminated with a short optical pulse and an acoustic pressure signal is generated via the thermoacoustic effect.^2,3 Wide-band ultrasonic transducers are used to measure the pressure signal at a number of locations outside the object and the measured data serve as the input to numerical image reconstruction algorithms that estimate the total absorbed optical energy density. The utility of PAT has been demonstrated in a number of in vivo studies of biological structure and function^5–7 and in facilitating a number of medically relevant diagnostic tasks.^2,8–13

Functional imaging of mouse brains has been successfully demonstrated with PAT.^5,6 For larger animals with thicker skulls, PAT has been demonstrated to provide qualitatively useful images of brain structures;^14–17 however, the transmission of ultrasonic waves through the skull bone(s) induces strong changes to the photoacoustic signal through the processes of absorption, strong reflections, and longitudinal-to-shear mode conversion. These effects can result in distortions and degraded spatial resolution in the reconstructed images. Time-reversal algorithms for PAT have been shown to compensate for density variations, dispersion and absorption,¹⁸ but they do not account for media that support shear wave propagation.

Under the assumption that the to-be-imaged object is embedded in a planar, layered medium, the effects of absorption, dispersion, and shear-mode conversion can be compensated for using a recently derived PAT reconstruction formula.¹⁹ The reconstruction formula, which is based on the angular spectrum representation of the photoacoustic field, estimates the Fourier components of the absorbed optical energy density. In this paper, the PAT reconstruction formula valid for planar, layered media is used to evaluate the role of modeling shear waves in PAT. Specifically, the transmission function for a three layer system when shear waves are modeled is compared to the transmission function for a three layer system for a system that does not compensate for shear waves. Comparisons are made based on the incident angle of the field in the plane wave representation as well as on the thickness of the solid layer. Object-independent error functions are introduced to quantify the importance of modeling shear waves for PAT in this geometry. We show that, when shear waves are omitted from the imaging model, a large number of Fourier components of the object will be mis-estimated.

2.

Background

The assumed geometry of our system is shown in Fig. 1. The to-be-imaged object is contained in layer 2, characterized by speed of sound $c_2$ and density ${\rho }_2$, where both are consistent with biological tissue (e.g., a fluid). The photoacoustic signal is detected in layer 0, on the plane $z=0$, characterized by thickness $d_0$, speed of sound $c_0$, and density ${\rho }_0$. Layer 1 is assumed to be a solid, thus supporting shear waves, with thickness $d_1$. The solid is assumed to be dispersive and absorptive for both longitudinal and shear waves.

Fig. 1

A diagram of the layered system. Ultrasonic detection occurs in the plane $z=0$. Layers 2 and 0 are assumed to be fluids, with densities ${\rho }_m$ and speed of sound $c_m$. Layer 0 has thickness $d_0$. Layer 1 is assumed to be a solid, bone, with density ${\rho }_1$, longitudinal speed of sound $c_l$, shear speed of sound $c_s$, thickness $d_1$, and longitudinal and shear absorption coefficients ${\overline{\alpha }}_l$ and ${\overline{\alpha }}_s$. The dispersive speeds of sound are referenced to 1 MHz. Layer 2 is assumed to be unbounded in the negative $z$-direction.

In PAT, the measured acoustic pressure is used to infer the absorbed optical energy density, denoted by $A(\mathbf{r})$. The space-frequency domain representation, $\tilde{p}(\mathbf{r},\omega )$, of the acoustic field is related to the space-time representation, $p(\mathbf{r},t)$, by a Fourier transform, viz

Let $\overline{p}(k_x,k_y,\omega )$ denote the two-dimensional (2-D) spatial Fourier transform of the pressure data $\tilde{p}(x,y,z,\omega )$ evaluated on the measurement plane $z=0$:

Dispersion and absorption in the bone are both important effects that must be modeled to develop an accurate imaging model for transcranial PAT. For many tissues of interest, it is known²² that the linear absorption coefficient, as a function of temporal frequency, obeys a power law: ${\alpha }_m(\omega )={\overline{\alpha }}_m{|\omega |}^y$. For $y\ne 1$, the dispersion relationship is given by²³

Fig. 2

(a) The longitudinal speed of sound as a function of frequency, (b) for the solid layer and the shear speed of sound. The medium is assumed to have known speeds of sound $c_l=2900\mathrm{m}/\sec$ and $c_s=1450\mathrm{m}/\sec$ at 1 MHz. The absorption coefficients are assumed to be ${\overline{\alpha }}_l=81.13{(\mathrm{Mrad}/\sec )}^{-0.93}{\mathrm{m}}^{-1}$ and ${\overline{\alpha }}_s=162.3{(\mathrm{Mrad}/\sec )}^{-0.93}{\mathrm{m}}^{-1}$. The speeds of sound averaged across the measurement band are ${\overline{c}}_l=2785\mathrm{m}/\sec$ and ${\overline{c}}_l=1392\mathrm{m}/\sec$. At 1 MHz, ${\alpha }_l(\omega )=170\mathrm{Np}/\mathrm{m}$ and ${\alpha }_s(\omega )=341\mathrm{Np}/\mathrm{m}$.

In this paper, we assumed a three-layer structure (see Fig. 1). The bottom layer was assumed to be soft tissue, unbounded in the $-z$ direction, with speed of sound $1483\mathrm{m}/\sec$ and density $1000\mathrm{kg}/{\mathrm{m}}^3$. The middle layer was assumed to be bone (a solid), with longitudinal speed of sound of $2900\mathrm{m}/\mathrm{s}$, shear speed $1450\mathrm{m}/\sec$ at 1 MHz and density $1900\mathrm{kg}/{\mathrm{m}}^3$. Both shear and longitudinal waves were assumed to obey a power law in the middle layer (bone) with $y=0.93$ and ${\alpha }_l(\omega )=170\mathrm{Np}/\mathrm{m}$ and ${\alpha }_s(\omega )=341\mathrm{Np}/\mathrm{m}$ at 1 MHz. The choice of $y=0.93$ falls within the reported ranges for power-law fitting.²⁶ The choices of ${\alpha }_1$ and ${\alpha }_s$ were chosen so that ${\alpha }_l$ fell within the range of reported absorption coefficients at 1 MHz and that the dispersive speeds of sound (shear and longitudinal) grew at similar rates. The top layer was assumed soft tissue, with speed of sound $1520\mathrm{m}/\sec$ and density $1100\mathrm{kg}/{\mathrm{m}}^3$. The thickness of the solid layer was varied in the simulations. The pressure wavefield was assumed to be recorded in a layer that is matched to the top layer a distance $d_0=1\mathrm{mm}$ away from the tissue/bone edge.

3.

Transfer Function Analysis

The manner in which modeling shear waves in the PAT image reconstruction algorithm affects image quality in transcranial PAT is best understood by analyzing the ATF, $T({\mathbf{k}}_{\parallel },\omega )$, in Eq. (5). Let ${\overline{T}}_s(\theta ,\omega )=T(\omega \sin (\theta )/c_2,\omega )$ denote the ATF for the layered medium when shear waves are modeled in layer 1. Let ${\overline{T}}_l(\theta ,\omega )=T(\omega \sin (\theta )/c_2,\omega )$ denote the ATF for the layered medium when shear waves are not modeled in layer 1.

Plots of the magnitude and phase of ${\overline{T}}_s$ and ${\overline{T}}_l$ as a function of the thickness of the solid layer are shown in Figs. 3 and 4, with each panel representing a different incident planewave angle. The assumed angular frequency is $\omega =2\pi ·1.5\mathrm{MHz}$. In Figs. 3(a) and 4(a), the incident planewave angle is 1 deg. One notes that $T_s$ and $T_l$ are similar for all thicknesses. In the limit that the thickness goes to zero, the magnitude of the ATF tends toward 0.9401, the transmission coefficient for normal incidence for a two-layer medium with properties consistent with layer 0 and layer 2. In Fig. 4(b)–4(d), the ATFs are shown for an incident angle of 15, 30, and 45 deg, respectively. One sees that $T_s$ and $T_l$ differ greatly as the incidence angle increases for any thickness of the bone. For any incidence angle except zero degrees (normal incidence), the propagation of shear waves will have some effect on the transmission of acoustic energy from layer 2 to layer 0, although that effect is best seen for larger angles. Moreover, one notes that ${\overline{T}}_l$ at 30 and 45 deg falls off rapidly with thickness.

Fig. 3

The magnitude of the ATFs, ${\overline{T}}_l(\theta ,\omega )$ (black dashes) and ${\overline{T}}_s(\theta ,\omega )$ (blue line), as a function of the thickness of the solid layer for planewaves at $\omega =2\pi ·1.5\mathrm{MHz}$ at an incidence angle of (a) 1 deg, (b) 15 deg, (c) 30 deg, and (d) 45 deg.

Fig. 4

The unwrapped phase of the ATFs, ${\overline{T}}_l(\theta ,\omega )$ (black dashes) and ${\overline{T}}_s(\theta ,\omega )$ (blue line), as a function of the thickness of the solid layer for planewaves at $\omega =2\pi ·1.5\mathrm{MHz}$ at (a) an incidence angle of 1 deg, 15 deg, 30 deg, and (d) and 45 deg..

Plots of the magnitude and phase of ${\overline{T}}_s$ and ${\overline{T}}_l$ as a function of angle are shown in Figs. 5 and 6 for four different temporal frequencies. In panel a of each figure, the assumed angular frequency is $\omega =2\pi ·0.5\mathrm{MHz}$. In panel b, $\omega =2\pi ·1.0\mathrm{MHz}$, in panel c, $\omega =2\pi ·1.5\mathrm{MHz}$ and in panel d, $\omega =2\pi ·2.0\mathrm{MHz}$. In all four panels, the thickness of the solid is assumed to be 2 mm. One notes that both the magnitude and phase of $T_l$ and $T_s$ are nearly identical for angles less than 25 deg in panels b, c, and d. Beyond that angle, the ATFs differ significantly. The critical angle for the fluid-solid interface found at 1 MHz is ${\theta }_c={\sin }^{-1}\frac{1483}{2900}=30.75\mathrm{deg}$, that is, acoustic planewaves incident on the interface at an angle greater than ${\theta }_c$ correspond to non-propagating longitudinal waves in the solid layer. Note that, due to the assumed dispersion relationship, the critical angle decreases with increased frequency from 1 MHz. Physically, for $\theta >{\theta }_c$, propagating shear waves generated at the fluid/solid interface between layers 2 and 1, traveling at $1450\mathrm{m}/\sec$ at 1 MHz, transmit energy from one side of the solid to the other. The longitudinal waves are evanescent and decay rapidly from the interface between layers 2 and 1. When shear waves are not modeled, there is no mechanism for these high spatial frequency waves to transmit from one fluid medium (layer 2) to the other (layer 0). Accordingly, PAT reconstruction algorithms for planar layered medium that do not account for the propagation of shear waves will yield images with a inferior/lower lateral spatial resolution than algorithms that account for shear waves. These findings generalize beyond the choice of specific system parameters chosen in this paper. Varying the values of ${\alpha }_m$ and $y$ will only affect the rate at which the ATFs fall off as a function of thickness in Fig. 3, but the relationship between $T_s$ and $T_l$ will remain the same. Likewise, using different values for the speeds of sound or densities in the solid and fluid layers will change the critical angle, ${\theta }_c$, but not the fact that $T_s$ and $T_l$ are in good agreement when the incident angle is less than the critical angle.

Fig. 5

The magnitude of the ATFs, ${\overline{T}}_l(\theta ,\omega )$ (black dashes) and ${\overline{T}}_s(\theta ,\omega )$ (blue lines), as a function of incident angle, for a bone layer with properties $d_1=3\mathrm{mm}$ and (a) $\omega =2\pi ·0.5\mathrm{MHz}$, (b) $\omega =2\pi ·1.0\mathrm{MHz}$, (c) $\omega =2\pi ·1.5\mathrm{MHz}$, and (d) $\omega =2\pi ·2.0\mathrm{MHz}$.

Fig. 6

The unwrapped phase of the ATFs, ${\overline{T}}_l(\theta ,\omega )$ (black dashes) and ${\overline{T}}_s(\theta ,\omega )$ (blue lines), as a function of incident angle, for a bone layer with properties $d_1=3\mathrm{mm}$ and (a) $\omega =2\pi ·0.5\mathrm{MHz}$, (b) $\omega =2\pi ·1.0\mathrm{MHz}$, (c) $\omega =2\pi ·1.5\mathrm{MHz}$, and (d) $\omega =2\pi ·2.0\mathrm{MHz}$.

4.

Error Analysis

It is clear that the ATF when shear waves are accounted for differs from the ATF when shear waves are not accounted for. These differences will manifest when attempting to image objects through a solid layer in two ways. Since shear waves travel more slowly than longitudinal waves in the solid layer, misestimating the time-of-flight of the pressure field in the solid layer will lead to phase errors in the estimated Fourier components of the object. Likewise, errors in the estimated amplitude of the Fourier components of the object will occur because of not accounting for longitudinal-to-shear wave conversion and vice-versa at solid-fluid interfaces. Both effects will lead to distortions in the image.

To better understand the way in which not accounting for shear waves affects image fidelity, two error measures are introduced that are object-independent:

In Figs. 7Fig. 8Fig. 9–10, the amplitude and phase errors are shown in the $k_x-k_z$ plane for four different thicknesses of the solid layer. We assumed that the maximum measured frequency of the acoustic field in each case was 7.5 MHz, setting the radius of the Ewald measurement sphere at $k_{\max }=2\pi ·7.5\mathrm{MHz}/c_2=31,776{\mathrm{m}}^{-1}$. Note that $E_a$ approaches $-1$ when $T_l$ is near zero, i.e., the model that does not account for shear waves predicts a null value for the transmission amplitude.

Fig. 7

The amplitude error, $E_a$, and the phase error, $E_p$, inside the Ewald measurement sphere, for a solid with thickness $d=10\mu \mathrm{m}$.

Fig. 8

The amplitude error, $E_a$, and the phase error, $E_p$, inside the Ewald measurement sphere, for a solid with thickness $d=100\mu \mathrm{m}$.

Fig. 9

The amplitude error, $E_a$, and the phase error, $E_p$, inside the Ewald measurement sphere, for a solid with thickness $d=1\mathrm{mm}$. The maximum positive error in amplitude was set to be 3 to improve the dynamic range in the image.

Fig. 10

The amplitude error, $E_a$, and the phase error, $E_p$, inside the Ewald measurement sphere, for a solid with thickness $d=1\mathrm{cm}$. The maximum positive error in amplitude was set to be 3 to improve the dynamic range in the image.

In Fig. 7, the amplitude and phase error are shown when the solid layer is assumed to be 10 μm thick. One notes that the amplitude and phase error are relatively insignificant for small values of $k_x$ for all $k_z$. However, larger errors exist for both the amplitude and phase near the line $k_z=0$, with the size of the error growing for both amplitude and phase as $k_x$ gets larger. In Fig. 8, the amplitude and phase error are shown when the solid layer is assumed to be 100 μm thick. Even for this relatively thin layer, the amplitude distortions are already significant. The boundary between green and purple in the amplitude plot or between light blue and purple in the phase plot corresponds approximately to 28 deg, or near the critical angle, ${\theta }_c$, for the fluid-solid interface. For larger incidence angles, the amplitude errors near $-1$ and the phase errors approach $\pi$.

In Fig. 9, the amplitude and phase error are shown when the solid layer is assumed to be 1 mm thick, and in Fig. 10, the amplitude and phase error are shown when the solid layer is assumed to be 1 cm thick. In both cases, which are relevant to primate brain imaging, one sees that the amplitude error is near $-1$ for all incidence angles above ${\theta }_c$. This severely limits the ability of a non-shear wave based method for resolving large transverse frequency components of the to-be-imaged object in a planar detection geometry. In all figures, the maximum positive error in amplitude was set to be 3. The large positive error locations correspond to points at which $T_s$ approaches 0 [see, for example, Fig. 5(b) at 30 deg].

To demonstrate the manner by which the errors in the ATFs due to neglecting shear waves can result in image artifacts, computer-simulation studies were performed. Because the layered medium problem is inherently 2-D, only 2-D simulations were run. A numerical phantom representing the object, $A(\mathbf{r})$, was considered that consisted of three uniform disks of radii of 2, 3, and 4 mm, located within the bottom layer of the structure with centers along the line $x=0$ [see Fig. 11(a)]. In these simulations, the solid layer was assumed to be the top edge of the image (i.e., the solid layer begins at the line $z=25\mathrm{mm}$). The disks were quasi-bandlimited by convolving each with a two-dimensional Gaussian function of radius 400 μm. Uniform disks were employed in the simulations because their 2-D Fourier transforms are known analytically, and therefore the simulated pressure data could be computed without numerical approximation. The object was assumed to be embedded in a the same three-layer structure as was described at the end of Sec. 2 (see Fig. 1) and was used in generating the error maps in Figs. 7Fig. 8Fig. 9–10.

Fig. 11

(a) An image of the numerical phantom-three band-limited disks. (b) The reconstruction of those disks, found using Eq. (10), is shown when the solid layer is assumed to be $d=100\mu \mathrm{m}$. (c) Reconstruction of those disks, found using Eq. (10), is shown when the solid layer is assumed to be $d=1\mathrm{mm}$.

Simulated pressure data were computed by use of $T_s({\mathbf{k}}_{\parallel },\omega )$, which properly modeled the effects of shear wave propagation in the solid layer. From these data, images were reconstructed by use of an inverse Fourier transform method that neglected shear wave propagation and employed the transfer function $T_l({\mathbf{k}}_{\parallel },\omega )$. A non-uniform sampling scheme in $\omega$ was assumed to avoid interpolation-based errors in the reconstructed images. This prevented the confounding of model-mismatch errors (from neglecting shear waves), with errors associated with the specific reconstruction model used. For example, a Fourier-based inversion scheme, such as the one found in Ref. 19 requires interpolation of the data in the $k_z$-direction. Given that a known subset of Fourier components are misestimated when neglecting shear waves, inclusion of these data will result in an error that is unrelated to the errors associated with the neglect of shear waves. This assumption results in a reconstruction formula of the form

In Fig. 11(b), the reconstruction using Eq. (10) is shown for the case when $T_s$ and $T_l$ correspond to the solid (bone) layer being 100 μm thick. One notes that the reconstructed disks have slight errors, most prominent near the left and right edges of each disk. The reconstructions also contain distortions in the background in the form of vertical ‘bands’. Both the errors in the disks and the bands are a result of the errors in misestimating the Fourier components along the line $k_z=0$ as well as those that correspond to waves incident at an angle larger than ${\theta }_c$.

Figure 11(c) shows the reconstruction using Eq. (10) is shown for the case when $T_s$ and $T_l$ correspond to the solid (bone) layer being 1 mm thick. In this case, the minor distortions along the right and left edges of the disks are much more pronounced, resulting in what appear to be small objects near the disks. Likewise, the vertical banding is significantly more pronounced. One notes that these errors—associated with misestimating components in the $k_z$-direction—are more pronounced because of the larger errors associated with the thicker bone layer (see Figs. 8 and 9). In both Fig. 11(b) and 11(c), the horizontal and near-horizontal boundaries of the disks are estimated much better than the vertical boundaries. This is consistent with the fact that neglecting shear waves results in very little error in the Fourier components corresponding to propagating longitudinal waves in the solid layer (i.e., waves incident at an angle less than ${\theta }_c$). One notes that previous studies in transcranial PAT ^14–17 did not possess these banding artifacts; however, the measurement geometry was not planar in those studies, which can result in different manifestations of the artifacts.

5.

Summary and Discussion

In this work, a recently-introduced PAT image reconstruction algorithm, for use when the optical absorber is embedded in a layered medium that supports shear waves in some of the layers, was employed to determine the role shear waves play in PAT image formation. Using system parameters relevant to transcranial PAT, it was shown that, when shear waves are omitted from the imaging model, a smaller number of spatial frequency components of the object can be estimated compared to situations when shear waves are included. The spatial frequency components that are no longer accurately able to be estimated are components parallel to the layered medium and consistent with planewave components incident beyond the critical angle of the solid-fluid interface. These mis-estimated Fourier components result in distortions in the reconstructed images. The strength of the distortions are directly related to the thickness of the solid (bone) layer.

The Fourier-based reconstruction method for inverting pressure data by disregarding shear waves in PAT imaging algorithms is mathematically equivalent to standard time-reversal (TR) and filtered backprojection (FBP) algorithms for PAT when the object of interest is embedded in a layered medium and acoustic detection is made on a continuous measurement surface. Since standard formulations of TR and FBP algorithms are insensitive to shear waves, images produced by these methods will exhibit the same systematic error as presented here. One notes that the geometry and type of transducers used to approximate a continuous measurement surface will also affect the systematic error introduced by not including shear waves.

The electrical and spatial impulse responses due to transduction and the finite size of each transducer, respectively, can be modeled as filters.^20,28 In the Fourier representation of the photoacoustic field presented in this work, the effects from the transducers would be modeled by a multiplicative factor in the denominator of the right hand side of Eq. (5). The electrical impulse response often serves to diminish very low temporal frequency information, which corresponds to eliminating data from the center of Ewald’s sphere. Depending on the size and shape, the spatial pass-band of the transducer may diminish data from Ewald’s sphere that correspond to high spatial frequencies, that is, data for which the largest errors occur by not accounting for shear waves. Much wider pass-bands in both spatial and temporal frequencies occur when using planar, polymer-film ultrasound sensors²⁹ to acquire the photoacoustic signal.

For the reconstruction formula that compensates for shear waves, estimates of $\mathcal{A}$ in the $k_z$-direction are non-uniformly sampled due to the relationship found in Eq. (6). To attain estimates of $A(\mathbf{r})$ from the non-uniformly sampled data, either interpolation of the data to a uniform grid in $k_z$ followed by an FFT is required, or both steps can be performed at once using a non-uniform FFT (NUFFT).³⁰ In either case, the interpolation is performed along vertical lines in the Ewald sphere. One notes that the error maps shown in Figs. 7Fig. 8Fig. 9–10 exhibit errors that are non uniform in the $k_z$-direction.

There are a number of open problems that remain to be addressed in regards to accounting for shear wave physics in transcranial PAT. The manner in which neglecting shear waves affects image quality when different image reconstruction algorithms and non-planar imaging geometries are employed has not been investigated systematically. The non-uniform shape of the skull will certainly induce different types of artifacts than the ‘bands’ presented in this work. The interplay between transducer placement/geometry and the errors described in this manuscript associated with neglecting shear waves is also an important topic of future investigation.

In some cases, the properties of the solid layer may not be determined accurately. One notes that mis-estimation of the acoustic absorption properties in the bone layer will affect both the phase (through using the wrong dispersion model) and the amplitude of the Fourier components of the object. Methods for estimating the acoustic absorption properties of bone that are sufficiently accurate to facilitate transcranial PAT applications remains an ongoing research topic.

When the shear wave properties are not known, it may be possible to develop an iterative reconstruction algorithm which uses the error maps found in Sec. 4 to create reconstructions using only the ‘minimal error’ components of the measured data. One notes that the regions of the Ewald sphere that possess minimal error are bounded by the critical angle of the solid-fluid interface for longitudinal waves when the solid layer is on the order of a millimeter or larger. That is, the region of estimable Fourier components depends only on the properties of the longitudinal waves in the medium. The development of such an algorithm remains as future work.