## Abstract

Photoacoustic (PA) imaging is based upon the generation of an ultrasound pulse arising from subsurface tissue absorption due to pulsed laser excitation, and measurement of its surface time-of-arrival. Expensive and bulky pulsed lasers with high peak fluence powers may provide shortcomings for applications of PA imaging in medicine and biology. These limitations may be overcome with the frequency-domain PA measurements, which employ modulated rather than pulsed light to generate the acoustic wave. In this contribution, we model the single modulation frequency based PA pressures on the measurement plane through the diffraction approximation and then employ a convolution approach to reconstruct the sectional image slices. The results demonstrate that the proposed method with appropriate data post-processing is capable of recovering sectional images while suppressing the defocused noise resulting from the other sections.

© 2011 OSA

## 1. Introduction

In recent years photoacoustic (PA) imaging has experienced rapid development in the biophotonics community [1–7]. PA takes the advantage of both the optical and ultrasonic imaging modalities, promising high-resolution structural and functional imaging of tissue [8]. Pulsed lasers are generally used to induce ultrasonic waves that originate from transient absorption of the propagating light pulse creating a local tissue expansion. The local expansion generates an ultrasonic pulse that travels to the tissue interface for facile detection. The pressure of PA wave is detected by ultrasound transducers positioned outside or on the surface of tissue and can be used to reconstruct the origin of the initial PA pulse. The technique is thus suitable for quantifying tissue properties or identifying pathological structures. PA has been demonstrated to image the blood vasculature in both small animals and humans via high resolution [9]. However, time-domain PA requires a short laser pulses in order to fulfill the assumptions upon which PA imaging is based [10]. In addition, expensive and bulky pulsed lasers with high peak powers can limit time-domain PA applications in medicine and biology.

Frequency-domain (FD) based PA methodologies may offer an attractive alternative measurement scheme although its signal-to-noise ratio (SNR) is 20dB - 40 dB worse than that of pulse system [11, 12]. In this approach, instead of a pulse of laser light, the tissue or object surface is illuminated with intensity-modulated continuous-wave (CW) laser or laser diode as an excitation source. Generally the modulation frequency in FD-based PA imaging system is scanned over the frequency range of interest to preserve full axial (depth) resolution, but this method may be time-consuming and can complicate instrumentation. Recently, a single-frequency detection using amplitude-based axial resolution has been demonstrated for imaging subsurface blood vessels in tissues [12]. However, in order to localize the PA sources one needs to sweep the modulation frequency within the bandwidth of ultrasonic transducer.

Herein we propose a convolution approach for the reconstruction of sectional images in FD based PA imaging using the complex information (i.e. both amplitude and phase). In our method, the PA pressure on the measurement plane is modeled using diffraction approximation and the inversion is implemented by performing the convolution operation and appropriate post-processing. This approach is similar to digital holography [13–15], and is the first time, to the best our knowledge, to be introduced to the FD based photoacosutic tomography. In comparison to digital holography PA based methods have a significant advantage in that the ultrasound amplitude attenuation is two orders of magnitude less than that of a propagating light imaging, allowing the acoustic waves to travel long distances in tissue without significant distortion or attenuation. Thus, the technique may be especially suitable to image deep biological tissues if incident fluence rates can be reasonably controlled and interference from other light absorbing structures not within the focus plane can be accounted for. To remove the defocused noise from the other sections, only real part of the reconstructed sectional image was used. Then the finding of maximum focus measure instead of sweeping modulation frequency was applied to determine the depth information. Combined with the light propagation model in tissue, the absorption coefficients of blood vessels can be quantitatively retrieved, making it possible for functional imaging.

## 2. Methods

#### 2.1 Photoacoustic wave equation

In the frequency domain, the generation and propagation of laser-induced acoustic wave through highly scattering media (such as biological tissue) can be described by the following Helmholtz-like equation [10, 16, 17]:

#### 2.2 Light propagation in scattering media

The gold standard for describing light propagation in turbid media is the radiative transport equation (RTE). However, a direct solution of RTE is a challenging task even with high performance computers. Instead, an approximation to the radiative transport equation is used and has been accomplished in biological tissues [18]. In this approximation, the multiply scattered light intensity is described by the diffusion equation, which is given in the frequency domain over a 3D bounded domain Ω by [19, 20]

#### 2.3 Reconstruction of sectional images by the convolution approach

The PA wave on the measurement plane $\left(x,y,z=0\right)$ can be expressed as the assembly of the spherical waves which converge at the different points $(\zeta ,\eta ,z)$ within the tissue. The measured acoustic wave $p\left(x,y,\omega \right)$ can be computed by the diffraction formula:

#### 2.4 Data post-processing

Suppose the focused plane is located at ${z}_{1}$_{,} as shown in Fig. 1
. A reconstruction of this sectional image requires recovering $A\text{\Phi}\left(\zeta ,\eta ,{z}_{1}\right)$from the PA wave of $p\left(x,y,\omega \right).$ The conventional approach to obtain sectional images involves computing the convolution of PA wave with the conjugate of the impulse response $({g}^{*})$at the focused position ${z}_{1}$:

The first term above represents the focused section at ${z}_{1}$ and the second term refers to as the defocused noise contributed by acoustic waves generated at the other sections. Hence, the reconstructed sectional image is a mixed distribution of the focused image points and the defocused noise. The reconstructed sectional image will be blurred by different degrees depending on the strength and distance of contributing sources from the focused section. Note also that the first term in Eq. (7) can be regarded as real due to the negligible phase shift of exciting laser light (see the numerical simulation section), while the second term is complex. By retaining the real part of Eq. (7), the imaginary part of the defocused noise will be completely removed since it contains no meaningful signal content at the focus plane. Suppose a section located at ${z}_{1}$depth contains blood vessels with orders of magnitude larger absorption than other tissue sections, and the maximum absorbed optical energy density would be deposited at the ${z}_{1}$ section. With the proposed convolution approach, the maximum absorbed optical energy density is reconstructed at the ${z}_{1}$section. By projecting the maximum reconstructed absorbed light energy density along the $z$ direction, the ${z}_{1}$ depth information can be obtained. Combined with the light propagation model, the map of optical absorption coefficient may be recovered. Figure 2 shows a block diagram of reconstruction of sectional images in the FD based photoacoustic imaging process.

## 3. Numerical simulations

As for the evaluation of the excitation fluence distribution, we assume a simple case of a slab tissue shown in Fig. 1 that can later be modified for different geometries. Simulated data of fluence was acquired on a 1.28 x 1.28 x 2.0 cm^{3} domain. The slab tissue was modeled as homogeneous with the following values: absorption coefficient ${\mu}_{a}$of 0.02 cm^{−1}; reduced scattering coefficient ${\mu}_{s}\text{'}$ of 10.0 cm^{−1}; reflection coefficient of 0.433 on the top surface and 0 on all other surfaces. While more complicated and informative sources such as patterned modulated incident light could be used, a homogeneously distributed area source modulated at 50MHz was simulated to illuminate the top surface of the chosen domain. Figure 3
shows the calculated complex photon fluence at the excitation wavelength of 785 nm as a function of depth, computed using finite element method. We can see from Fig. 3 that the amplitude of the complex photon fluence decreases exponentially and its phase increases slowly with respect to the depth. When compared to the acoustic wave phase delay shown below, the initial phase value of the induced PA wave can be assumed to be zero. For the pressure wave measurements, the slab tissue and the acoustic transducer are immersed in water and the distance between the top surface of the slab tissue and the measurement plane is 4.0 cm, as shown in Fig. 1.

Reconstruction was performed on a domain of interest divided into 20 sections$(M=20)$, where two distinct sections, ${z}_{16}=4.8$cm and ${z}_{18}=4.9$cm away from $z=0$measurement plane, include absorptive targets mimicking blood vessels with optical absorption coefficient of 3.0 (cm^{−1}). For these numerical simulations, the coefficient $A=i{k}_{0}\frac{{c}_{0}\beta}{{C}_{p}}$ in Eq. (1) is set to 1 for simplicity. Figure 4(a)
and 4(b) present the normalized amplitude and phase of the simulated acoustic pressure field on the measurement plane, respectively, from which it is difficult to decipher the structure of blood vessels. Gaussian noise of 20 dB was first added to the simulated PA pressure to simulate experimentally measured data. The reconstructed 3D map was sampled in intervals of $\Delta \zeta =\Delta \eta =\text{1}00\mu \text{m}$, which defines transverse resolution in the $\zeta $and $\eta $directions, respectively, and a longitudinal sampling interval of $\text{5}00\mu \text{m}$, which defines resolution in the $z$ direction.

Figure 5(a) and 5(b) show the reconstructed sectional images at ${z}_{16}=4.8$cm and ${z}_{18}=4.9$cm, respectively, without data post-processing. Their corresponding absorptive coefficient profiles through $\zeta =0.8$ cm are shown in Fig. 5(c) and 5(d), respectively. The figures show that the absorption coefficients of simulated blood vessels are clearly reconstructed in the corresponding focused sections. We can also observe the strong defocused noise originating from the interference in sections adjacent to the focused sections.

With the data post-processing, the reconstructed focused sections are weakly influenced by the defocused noise, as shown in Fig. 6
. Figure 7
shows the reconstructed sectional image at depth ${z}_{11}=4.55$ cm which contains no blood vessels and the similar resutls are obtained for other sections. The C-scan or *en face* images (a gray-scale plot of adding all the reconstructed sectional images) are shown in Fig. 8
, from which we can see that optical absorption coefficients of blood vessels have been reconstructed properly.

As a complementary of the blood vasculature, the lymphatic circulatory system returns water and proteins from various tissues back to the bloodstream. The lymphatics can be imaged using fluorescence optical microscopy [21], and near-infrared fluorescence imaging [22], and photoacoustic imaging [23], through intradermal injection of vital dye such as indocyanine green (ICG). With the administration of dye, the lymphatic vessels may have high optical absorption coefficients compared to the blood vessels. Figure 9(a)
and 9(b) show the reconstructed optical coefficient coefficients of blood vessels with ${\mu}_{a}$of 3.0 cm^{−1} located at 4.8 cm and lymphatic vessels with${\mu}_{a}$of 30.0 cm^{−1} located at 4.9 cm, respectively. Their corresponding profiles through $\zeta =0.8$cm are shown in Fig. 9(c) and 9(d), respectively. The corresponding C-scan or *en face* images are shown in Fig. 10
, from which we can see that both the blood vessels and lymphatic vessels were reconstructed simultaneously. Similar simulation results can be obtained under other modulation frequencies.

## 4. Discussion

As for the acoustic diffraction theory, it is reasonable by employing a small size transducer with relative large distance from the object and immersing the transducer and the imaging domain in water. In addition, the acoustic pressure, as given in Eq. (4), does not consider the photonic and thermal contributions to the system. This implies that the thermal and acoustic confinements are satisfied. It should be noted that both the thermal and stress confinements can be measured by the ratios of the squared photonic to the thermal and acoustic wave numbers, respectively [10]. For the parameters given in the numerical simulations as well as the typical thermal diffusivity of soft tissue as 0.11mm^{2}/s and the speed of light in tissue is on the order of 2 x10^{8} m/s, the calculated ratios of the photonic to the thermal and acoustic wave numbers are small and it follows that the thermal subsystem has little effects on the measurements [10]. If acoustic attenuation in tissue is not negligible, one may need to modify the photoacoutic wave Eq. (1). As photons experience attenuation in tissue, the photonic contribution becomes negligible and only the acoustic contribution remains at a large distance. Otherwise, one needs to consider the photonic contribution by modifying the impulse response,$\text{g}\left(\zeta ,\eta ,x,y,z\right)$. To obtain more accuracy distribution of photon fluence in biological tissues, high order spherical harmonics and simplified spherical harmonics models may be used to solve the RTE [24, 25]. On the other hand, using the Monte Carlo simulation to solve the RTE, one can obtain any desired accuracy [26].

In numerical simulations, we assume that the blood vessels at different sections are not overlapped along the *z* direction and thus this technique is suitable for imaging the shape of tissue of object. If the blood vessels overlap along $z$direction, based on the proposed convolution method with data post-processing techniques, the blood vessel with the maximum absorbed optical energy density in $z$direction would be constructed. The proposed approach is direct and fast. If one wants to reconstruct the overlapping structures, the Wiener filter approach, the Wigner distribution function approach and the inversion methods used in the field of optical scanning holography [14-15, 27] can be adapted and applied to the photoacoustic imaging. Both the Wiener filter and Wigner distribution function approaches work well for only two sections. In addition, the Wiener filter approach requires no noise virtually introduced in the imaging system, and the imaging reconstruction is a time-consuming process for the Wigner distribution approach. Fortunately, inverse imaging approaches, such as the conjugate-gradient based algorithm with Tikhonov regularization, are more effective in suppressing noise and retrieving the multiple sections. However, the sharp edges become blurred in the reconstructed sections, due to minimizing the total energy in the reconstruction algorithm. With the total-variation regularization and a nonnegativity constraint, the sharp edges can be well preserved in the reconstructed sections using a primal-dual Newton’s method combined with gradient projection. The minimum separation between two sections depends on the modulation frequency, SNR and the sampling size in the $z$ direction. For our experimental conditions, two sections with minimum separation of $\text{5}00\mu \text{m}$can be reconstructed.

In summary, we demonstrate that the distributions of the optical absorption are possibly to be reconstructed using the convolution approach with appropriate data post processing coupled with optical diffusion equation in the FD based PA imaging.

## References and Links

**1. **X. Wang, Y. Pang, G. Ku, X. Xie, G. Stoica, and L. V. Wang, “Non-invasive laser-induced photoacoustic tomography for structural and functional imaging of the brain in vivo,” Nat. Biotechnol. **21**(7), 803–806 (2003). [CrossRef]

**2. **H. P. Brecht, R. Su, M. Fronheiser, S. A. Ermilov, A. Conjusteau, and A. A. Oraevsky, “Whole-body three-dimensional opto-acoustic tomography system for small animals,” J. Biomed. Opt. **14**(6), 064007 (2009). [CrossRef]

**3. **D. Razansky, M. Distel, C. Vinegoni, R. Ma, N. Perrimon, R. W. Köster, and V. Ntziachristos, “Multispectral opto-acoustic tomography of deep-seated fluorescent proteins in vivo,” Nat. Photonics **3**(7), 412–417 (2009). [CrossRef]

**4. **Y. Sun and H. Jiang, “Quantitative three-dimensional photoacoustic tomography of the finger joints: phantom studies in a spherical scanning geometry,” Phys. Med. Biol. **54**(18), 5457–5467 (2009). [CrossRef]

**5. **L. V. Wang, “Prospects of photoacoustic tomography,” Med. Phys. **35**(12), 5758–5767 (2008). [CrossRef]

**6. **L. V. Wang, “Multiscale photoacoustic microscopy and computed tomography,” Nat. Photonics **3**(9), 503–509 (2009). [CrossRef]

**7. **B. T. Cox, S. R. Arridge, and P. C. Beard, “Estimating chromophore distributions from multiwavelength photoacoustic images,” J. Opt. Soc. Am. A **26**(2), 443–455 (2009). [CrossRef]

**8. **L. V. Wang, “Tutorial on photoacoustic microscopy and computed tomography,” IEEE J. Sel. Top. Quantum Electron. **14**(1), 171–179 (2008). [CrossRef]

**9. **B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt. **15**(2), 021314 (2010). [CrossRef]

**10. **N. Baddour, “Theory and analysis of frequency-domain photoacoustic tomography,” J. Acoust. Soc. Am. **123**(5), 2577–2590 (2008). [CrossRef]

**11. **A. Petschke and P. J. La Rivière, “Comparison of intensity-modulated continuous-wave lasers with a chirped modulation frequency to pulsed lasers for photoacoustic imaging applications,” Biomed. Opt. Express **1**(4), 1188–1195 (2010). [CrossRef]

**12. **K. Maslov and L. V. Wang, “Photoacoustic imaging of biological tissue with intensity-modulated continuous-wave laser,” J. Biomed. Opt. **13**(2), 024006 (2008). [CrossRef]

**13. **U. Schnars and W. P. O. Jutner, “Digital recording and numerical reconstruction of holograms,” Meas. Sci. Technol. **13**(9), R85–R101 (2002). [CrossRef]

**14. **X. Zhang, E. Y. Lam, and T. C. Poon, “Reconstruction of sectional images in holography using inverse imaging,” Opt. Express **16**(22), 17215–17226 (2008). [CrossRef]

**15. **E. Y. Lam, X. Zhang, H. Vo, T. C. Poon, and G. Indebetouw, “Three-dimensional microscopy and sectional image reconstruction using optical scanning holography,” Appl. Opt. **48**(34), H113–H119 (2009). [CrossRef]

**16. **Z. Yuan, C. Wu, H. Zhao, and H. Jiang, “Imaging of small nanoparticle-containing objects by finite-element-based photoacoustic tomography,” Opt. Lett. **30**(22), 3054–3056 (2005). [CrossRef]

**17. **Z. Yuan, Q. Wang, and H. Jiang, “Reconstruction of optical absorption coefficient maps of heterogeneous media by photoacoustic tomography coupled with diffusion equation based regularized Newton method,” Opt. Express **15**(26), 18076–18081 (2007). [CrossRef]

**18. **M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt. **28**(12), 2331–2336 (1989). [CrossRef]

**19. **E. M. Sevick, B. Chance, J. Leigh, S. Nioka, and M. Maris, “Quantitation of time- and frequency-resolved optical spectra for the determination of tissue oxygenation,” Anal. Biochem. **195**(2), 330–351 (1991). [CrossRef]

**20. **F. Fedele, J. P. Laible, and M. J. Eppstein, “Coupled complex adjoint sensitivities for frequency-domain fluorescence tomography: theory and vectorized implementation,” J. Comput. Phys. **187**(2), 597–619 (2003). [CrossRef]

**21. **P. Baluk, J. Fuxe, H. Hashizume, T. Romano, E. Lashnits, S. Butz, D. Vestweber, M. Corada, C. Molendini, E. Dejana, and D. M. McDonald, “Functionally specialized junctions between endothelial cells of lymphatic vessels,” J. Exp. Med. **204**(10), 2349–2362 (2007). [CrossRef]

**22. **J. C. Rasmussen, I. C. Tan, M. V. Marshall, C. E. Fife, and E. M. Sevick-Muraca, “Lymphatic imaging in humans with near-infrared fluorescence,” Curr. Opin. Biotechnol. **20**(1), 74–82 (2009). [CrossRef]

**23. **C. Kim, K. H. Song, F. Gao, and L. V. Wang, “Sentinel lymph nodes and lymphatic vessels: noninvasive dual-modality in vivo mapping by using indocyanine green in rats--volumetric spectroscopic photoacoustic imaging and planar fluorescence imaging,” Radiology **255**(2), 442–450 (2010). [CrossRef]

**24. **A. D. Klose and E. W. Larsen, “Light transport in biological tissue based on the simplified spherical harmonics equations,” J. Comput. Phys. **220**(1), 441–470 (2006). [CrossRef]

**25. **Y. Lu, B. Zhu, H. Shen, J. C. Rasmussen, G. Wang, and E. M. Sevick-Muraca, “A parallel adaptive finite element simplified spherical harmonics approximation solver for frequency domain fluorescence molecular imaging,” Phys. Med. Biol. **55**(16), 4625–4645 (2010). [CrossRef]

**26. **H. Shen and G. Wang, “A tetrahedron-based inhomogeneous Monte Carlo optical simulator,” Phys. Med. Biol. **55**(4), 947–962 (2010). [CrossRef]

**27. **X. Zhang and E. Y. Lam, “Edge-preserving sectional image reconstruction in optical scanning holography,” J. Opt. Soc. Am. A **27**(7), 1630–1637 (2010). [CrossRef]