1. Introduction

Non-invasive laser-induced photoacoustic tomography (PAT), as a hybrid imaging method, overcomes the resolution disadvantage of pure optical imaging and contrast/speckle disadvantages of pure ultrasound imaging [1–7]. It has been shown that PAT is a powerful imaging method for visualizing both the structural and functional information of biological tissues with excellent spatial resolution and satisfactory imaging depth [8, 9]. PAT is concerned with an inverse problem where a single short-pulsed light beam illuminates an object and the photoacoustic waves excited by thermo elastic expansion are measured using wideband ultrasound transducers in multiple locations around the object. The geometry of the object and spatial distribution of the absorbed optical energy density and acoustic property are recovered from the measured scattered fields using a reconstruction algorithm. To date various PAT reconstruction algorithms have been developed and applied to the detection of breast cancer, skin cancer, vascular diseases and brain tumors in small animals [10–14]. These algorithms, however, almost all rely on analytical solutions to the photoacoustic wave equation in a regularly shaped measurement configuration without appropriate boundary conditions applied. The major assumption made in these linear algorithms is that biological tissues are acoustically homogeneous which is not true in reality. To overcome these limitations, we have developed an iterative nonlinear reconstruction algorithm that is based on the finite element solution to the full Helmoholtz-like photoacoustic wave equation without the homogeneous acoustic property approximations [15]. Together with an iterative Newton method, our algorithm is able to precisely solve the Helmholtz wave equation and fulfill reliable inverse computation for an arbitrary measurement configuration [16–18]. Further, using simulated data [15] we show that our algorithm has an additional advantage compared to the existing ones: the images of both acoustic speed and absorbed optical energy density can be recovered simultaneously. It is known that there exist significant differences in acoustic properties between normal and tumor tissues and that tissue acoustic property images have the potential to differentiate benign from malignant breast lesions [19, 20]. Thus our PAT method offers a unique ability of reconstructing both tissue optical and acoustical properties in a single modality for a more complete diagnosis of diseases. In this paper, we demonstrate experimental evidence that the optical and acoustic images of heterogeneous media can indeed be obtained simultaneously using our reconstruction method.

2. Methods and materials

Since our reconstruction algorithm has been described in detail elsewhere [15–18], we only give a brief outline here for context. The frequency-domain wave equation for PAT imaging in an acoustically heterogeneous medium is written as [15],

where p is the acoustic pressure; k ₀ = ω/c ₀ is the wave number described by the angular frequency, ω and the speed of acoustic wave in a reference or coupling medium, c ₀; β is the thermal expansion coefficient; C_p is the specific heat; Φ is the absorbed optical energy density that is the product of optical absorption coefficient and incident optical fluence; O is a coefficient that depends on both acoustic speed and attenuation as follows [15]:

in which c is the speed of acoustic wave in the medium/tissue; α is the acoustic attenuation coefficient. The finite element discretization of Eq. (1) is stated as,

where the elements of the matrix A and b are expressed as

in which Φ , O and p are expressed by their real and imaginary components, ψ is the basis function, <> indicates integration over the problem domain and the second-order absorption boundary conditions are employed here [16],

where $\eta =\frac{-\mathit{ik}-\frac{3}{2\rho }+\frac{i3}{{8\mathit{k\rho }}^2}}{\frac{1-i}{\mathit{k\rho }}},\gamma =\frac{\frac{-i}{{2\mathit{k\rho }}^2}}{\frac{1-i}{\mathit{k\rho }}}$ and θ is the angular coordinate at a radial position ρ. To obtain a matrix equation capable of inverse solution, a combined Marquardt and Tikhonov regularization scheme is used for the inverse calculation [15],

where p o = ($p_1^o$,$p_2^o$,…,$p_M^o$ )^T, p^c = ($p_1^c$,$p_2^c$,…,$p_M^c$ )^T, and $p_i^o$ and $p_i^c$ are observed and computed complex acoustic field data for i=1, 2…, M boundary measurement locations, Δχ is the update vector for the optical and acoustic properties , ℑ is the Jacobian matrix formed at the boundary measurement sites, λ is a scalar and I is the identity matrix. Thus here the image formation task is to update optical and acoustic property distributions via iterative solution of Eqs. (3) and (6) so that a weighted sum of the squared difference between computed and measured acoustic pressure can be minimized.

The experimental system used was based on a mechanical scanning mechanism [16–18]. In this system, pulsed light with incident fluence 10 mJ/cm² well below the safety standard from a Nd: YAG laser (wavelength: 532nm, pulse duration: 3–6ns) were coupled into the phantom via an optical subsystem and generated acoustic pressure wave. A wide-bandwidth, 1 MHz transducer was used to receive the acoustic signals. The transducer and the phantom were immersed in a water tank. A rotary stage rotated the receiver relative to the center of the tank. One set of data was taken at 120 positions when the receiver was scanned circularly over 360°. The complex wavefield signal was first amplified by a preamplifier (Gain: 17dB, 0.03KHz~3MHz, Onda, Corporation, Sunnyvale, CA), and then amplified further by a Pulser/Receiver (GE Panametrics, Waltham, MA). A data acquisition board converted it into digital one that was fed to a PC. The radius of the receiver motion path was 25.4mm in the phantom experiments conducted.

3. Results and discussion

In the first two experiments, we embedded one glycerol powder-containing target (10mm in diameter) in a 30 mm-diameter solid cylindrical phantom. We then immersed the object-bearing solid phantom in a 50.8mm-diameter water background, as displayed in Fig. 1(a). The acoustic speed in the region of the target containing 10% (phantom 1) and 28% (phantom 2) glycerol was 1540 and 1620m/s, respectively, while the acoustic speed for the background was 1480m/s. For the third experiment, a 5mm- and a 10mm-diameter target (each contained 28% glycerol) were embedded into a 30 mm-diameter solid cylindrical phantom. Again the entire solid phantom was immersed in a 50.8mm-diameter water background, as shown in Fig. 1(b). In the final experiment, we embedded two nanoparticle-containing objects (3mm in diameter) in a 50 mm-diameter solid cylindrical phantom, as plotted in Fig. 1(c). Each target had 0.5nM nanoparticles. We then immersed the object-bearing solid phantom in a 110.6 mm-diameter water background. For all the experiments, the phantom materials used consisted of Intralipid as scatterer and India ink as absorber with Agar powder for solidifying the Intralipid and India ink solution. The optical absorption and scattering coefficients for the phantom were 0.004mm^-1 and 1.0mm^-1, respectively. The optical absorption contrast between the target and phantom is about 10 for experiments 1–3, and 20 for experiment 4.

 

Fig. 1. Experiment geometry for (a) phantoms 1 and 2, (b) phantom 3 and (c) phantom 4. The dimension units are millimeters

Download Full Size | PDF

Figures 2(a)–2(b) show the reconstructed acoustic and optical images for experiment 1, while Figs. 2(c)–2(d) present the recovered acoustic and optical images for experiment 2. We see that the object in each case is clearly detected. By estimating the full width at half maximum (FWHM) of the acoustic property, as displayed in Fig. 3, we found that the recovered size of these objects ranged from 9.1 to 9.3 mm, in good agreement with the actual object size of 10mm. In our experiments, the 10 and 28% glycerol containing target used gave an acoustic speed of 1540 and 1620m/s, respectively. As can be seen from Fig. 3, we are able to obtain quantitatively resolved images in terms of the acoustic property as well as the position, shape and size of the objects.

The optical images obtained (i.e., absorbed optical energy density) are dependent on the incident optical fluence, thus giving just a relative measure of the optical absorption property of the media. It is noted the existing PAT reconstruction methods can image only the distribution of absorbed optical energy density that is the product of both the local optical absorption coefficient and fluence distribution within the irradiated medium. It is well known that it is the tissue absorption coefficient that directly correlates with tissue structural and functional information such as blood oxygenation. We have developed a reconstruction algorithm that is based on the finite element solution to the photoacoustic wave equation and photon diffusion equation for recovering the optical properties of heterogeneous media. For the developed algorithm, the absorption coefficient can be iteratively reconstructed from the recovered absorbed energy density and the optical fluence that is computed by the diffuse optical tomography model using the radiant exposure as the source [18].

Figures 4(a) and 4(b) display the reconstructed acoustic speed and optical images for the third experiment having two targets. By estimating the FWHM of the acoustic profiles plotted in Fig. 5, we found that the recovered size of the two objects was 4.6 and 9.8mm, respectively. Again, both size values were in good agreement with the exact values of 5 and 10mm. As shown in Fig. 4, the reconstructed acoustic properties for these cases are also accurate relative to the exact value, which is further confirmed by the quantitative analysis shown in Fig. 5.

Figures 6(a)–6(b) present the reconstructed acoustic and optical properties of two objects having 0.5nM nanoparticles within each of the objects. By estimating the full width half maximum (FWHM) of the acoustic speed plotted in Fig. 7, we found that the recovered diameter size of these objects ranges from 2.8mm to 3.3mm, in good agreement with the actual object size of 3mm. In addition, it is noted the imaging resolution in Fig. 4 is not as good as that in Fig. 6 due to the low optical contrast between the targets and phantoms for experiments 1–3. Further, from Figs. 2, 4 and 6, we can also see that the image quality of both optical and acoustic images appears to be similar when different acoustic contrasts between the target and the background were used.

In summary, we have demonstrated experimental evidence that both optical and acoustic property images can be obtained using our photoacoustic imaging method. This capability of reconstructing two entirely distinct tissue parameters in a single modality may bring in significant clinical value for decision-making.

 

Fig. 2. Reconstructed acoustic speed (left column) and absorbed optical energy density (right column) images with different acoustic contrast levels with respect to the background. (a) and (b) for experiment 1. (c) and (d) for experiment 2. The axes (left and bottom) illustrate the spatial scale, in millimeters, whereas the color scale (right) records the acoustic speed, in mm/s, or absorbed optical energy density, in a relative unit.

Download Full Size | PDF

 

Fig. 3. Recovered acoustic speed profiles plotted along x=0mm (a), y=0mm (b) from the images shown in Figs. 2(a) and 2(c), respectively.

Download Full Size | PDF

 

Fig. 4. Reconstructed acoustic speed (left) and absorbed optical energy density (right) images for experiment 3. The axes (left and bottom) illustrate the spatial scale, in millimeters, whereas the color scale (right) records the acoustic speed, in mm/s, or absorbed optical energy density, in a relative unit.

Download Full Size | PDF

 

Fig. 5. Recovered acoustic speed profiles plotted along y=0mm (a), y=-8mm (b) from the image shown in Fig. 4(a).

Download Full Size | PDF

 

Fig. 6. Reconstructed acoustic speed (left) and absorbed optical energy density (right) images for experiment 4. The axes (left and bottom) illustrate the spatial scale, in millimeters, whereas the color scale (right) records the acoustic speed, in mm/s, or absorbed optical energy density, in a relative unit.

Download Full Size | PDF

 

Fig. 7. Recovered acoustic speed profile plotted along x=0mm from the image shown in Fig. 6(a).

Download Full Size | PDF