Open Access
Add to BibSonomyAbstract
In this study, we adopt a model-based correction method to reduce the finite aperture effect in photoacoustic tomography (PAT) – the tangential resolution deteriorates as the imaging point moves away from the circular scanning center. Such degradation in resolution originates from the spatial impulse responses (SIRs) of the used finite-sized unfocused transducer. Based on a linear, discrete PAT imaging model, the proposed method employs a spatiotemporal optimal filter designed in minimum mean square error sense to compensate the SIRs associated with an unfocused transducer at every imaging point; thus retrospective restoration of the tangential resolution can be achieved. Simulation and experimental results demonstrate that this method can substantially improve the degraded tangential resolution for PAT with finite-sized unfocused transducers while retaining the radial resolution.
©2010 Optical Society of America
1. Introduction
Photoacoustic imaging is a promising hybrid imaging modality that combines light and ultrasound to detect absorbed photons ultrasonically through the photoacoustic effect. It overcomes the resolution drawback of pure optical imaging due to overwhelming light scattering in biological tissues, and owns the most compelling features of both optics and ultrasound – high optical absorption contrast and sub-millimeter ultrasound resolution up to an imaging depth of centimeters [1,2]. When biological tissues, being illuminated by laser pulses, absorb the pulsed laser energy, photoacoustic signals are induced due to transient thermo-elastic expansion. These photoacoustic waves are then detected by ultrasonic transducers and mapped back to form images representative of optical absorption distribution. Photoacoustic imaging has shown promise in a wide spectrum of applications such as imaging of blood vessels [3], detection of melanoma and breast tumor [4,5], and monitoring of blood oxygenation [6,7].
There are three imaging modes – forward-, backward-, and tomography-mode – in photoacoustic imaging. Here, we intend to solve the problems associated with photoacoustic tomography (PAT) mode. In PAT, a planar circular scanning geometry is commonly used [2]. A short pulse is used to illuminate the imaging object from the top, and an ultrasonic transducer is circularly scanned around the object to collect photoacoustic signals. Then a back-projection algorithm is applied to reconstruct distribution of laser energy deposition. So far most rigorous PAT back-projection algorithms have been developed based on “ideal point detectors” which are capable of detecting photoacoustic waves with a very large acceptance angle and accurate time delay, and thus offering uniform spatial resolution in reconstructed images. Nonetheless, the small active area of point-like detectors leads to poor signal sensitivity. Hence, flat unfocused transducers with large detection area are commonly used in practice to address such an issue, sacrificing the receiving angle and signal delay accuracy due to the finite-length spatial impulse responses (SIRs, i.e., diffraction effect) associated with the transducers [8]. These compromises cause the “finite aperture effect” – the tangential resolution deteriorates as the imaging point moves away from the circular scanning center [9]. Such spatial-variant tangential resolution is not proper for imaging big objects, such as breast.
A virtual point detector idea with a positive or negative focused transducer has been proposed to solve such a problem [10,11]. By using this idea, the tangential resolution has been experimentally demonstrated to be improved more than three fold in PAT when the object is far away enough from the circular scanning center [10]. However, the tangential resolution is still spatially variant, and the improvement comes with the price of the degradation in radial resolution, as shown in the numerical study by Li [11].
In this study, we adopt a model-based correction method to improve the degraded tangential resolution for PAT with finite-sized unfocused transducers, with no need to change the transducer type. This method is based on a linear, discrete imaging model of PAT. Using this model, a spatiotemporal optimal filter designed in minimum mean square error sense is used to deconvolve the SIRs associated with an unfocused transducer at every imaging point; thus retrospectively restoring the degraded tangential resolution. Performance of our proposed method is first evaluated by simulations and further demonstrated with experimental PAT data.
2. Materials and methods
2.1 Spatial impulse response of a finite sized transducer
The spatial impulse response (SIR) concept is a time domain approach to model acoustic radiation and reception, and stems from linear acoustics where wave propagation is treated as a linear time-invariant system. Without consideration of the electrical impulse response of the transducer, for a Dirac photoacoustic source located at r, the detected signal by a finite-sized transducer at rT is defined as SIR. Based on linear acoustics, this SIR can be expressed as the integral of detected Dirac pulses with propagation delays over the transducer surface [12,13]:
where ST represents the active area of the transducer.According to Eq. (1), the SIR takes the form of a finite-length pulse, depending on the position of the source and the shape and size of the active area of the transducer. It accounts for the diffraction effect over the transducer’s active area; thus determines the directivity patterns of the transducer, i.e., detection angle. Figure 1 compares the SIRs of a point detector and a finite-sized unfocused transducer (i.e., a circular transducer which is commonly used in PAT) to an off-axis Dirac photoacoustic source. Note that in this study, the SIRs are all calculated numerically with the DREAM toolbox which is available from the website http://www.signal.uu.se/Toolbox/dream/ [8,13]. Compared with that of a point detector, the SIR of a circular transducer introduces limited detection angle, waveform distortion, and time delay error to the back-projection algorithm established on point detectors; thus forming the finite aperture effect.
Fig. 1 Illustration of normalized SIRs of a point detector and a finite-sized unfocused transducer (i.e., circular transducer) to an off-axis point source, and the problems caused by the SIR of the finite-sized unfocused transducer.
2.1 Two-dimensional imaging model of PAT
In this section, a two-dimensional (2D) discrete and linear imaging model of PAT with finite-sized unfocused transducer is built using the SIR concept and matrix formalism. Using this model, a model-based correction method is developed to compensate the finite aperture effect. In PAT, with the SIR representation, the output signal produced by a finite-sized transducer at a given scanning angle θl that arises from a point optical absorber at position r can be expressed as:
where ⊗ denotes a one-dimensional temporal convolution, he(t) is the electrical impulse response of the transducer, p(t) is the unit photoacoustic waveform generated by the point absorber which can be approximated with the derivative of a unit delta function, and o(r) is the optical absorption or optical energy deposition on the point absorber determining the amplitude of p(t).Figure 2 illustrates sampling in PAT data acquisition and sonogram used for matrix formalism of the 2D PAT imaging model along with Eq. (2). The region of interest (ROI), O, is divided into M by N image elements, and each element o(xn, ym) is a scalar representing the optical absorption of the absorber or the optical energy deposition on the absorber at position (xn, ym). The transducer is circularly scanned to acquire photoacoustic signals at different scanning angles, and the measurements are taken from L angles in 360-degree full view. Therefore, the full data set (i.e., PAT sonogram) Y of one complete circular scan consists of L discrete-time A-lines. The received A-line signal by the finite-sized transducer at scanning angle θl is denoted as a K by 1 vector yl , containing K samples. Re-writing Eq. (2) using superposition principle, yl can be expressed as a sum of photoacoustic signals from all absorbers in the ROI:
where hSIR(θl, xn, ym) is the SIR vector from the transducer position θl to o(xn, ym), h e and p are vectors of the electrical impulse response of the transducer and the unit photoacoustic waveform, respectively. After stacking the Y and O into column vectors of y and o, respectively, and using the notation in Eq. (3), the 2D discrete imaging model of PAT with finite-sized transducers can be expressed in the following matrix notation [14]:where e is a vector denoting additive noise, and S denotes the SIR matrix which is constructed with the K by M matrices for all combinations of scanning angles θl of the transducer and positions xn in ROI O:2.2 Spatiotemporal optimal filter: compensating for the SIRs
Based on Eq. (4), a spatiotemporal optimal filter K can be derived to reconstruct the PAT image where the finite aperture effect caused by the SIRs is compensated. In this study, this spatiotemporal optimal filter K is designed to minimize the reconstruction error betweenoand the reconstructed image in the mean square error sense. With the assumption that both oand eare Gaussian, mutually uncorrelated random vectors with covariance matrices Co and Ce, respectively, the spatiotemporal optimal filter K can be found [14]:
where T denotes matrix transpose.PAT image reconstruction with this spatiotemporal filter K can be viewed as a two-step reconstruction process. Without taking the electrical impulse response of the transducer into account, the filter K deconvolves the SIRs associated with the scanning unfocused transducer at every imaging point in the acquired PAT sonogram Y and also performs back-projection with the processed sonogram so that retrospective restoration of the tangential resolution can be achieved.
3. Results
3.1 Simulations
Simulations were performed to evaluate the efficacy of the proposed method on the improvement of tangential resolution for PAT with finite-sized unfocused transducers. The comparison between the proposed method and the back-projection (i.e., delay and sum) algorithm were carried out. Note that in the following, all the results from our proposed method were done without compensation for the electrical impulse response of the transducer in order to facilitate the comparison with the back-projection algorithm that does not compensate the electrical impulse response. All the simulations were done using the DREAM toolbox [8,13]. Six point absorbers located at different distances away from the scanning center were simulated to obtain the point spread functions for resolution examination. The point absorbers were located at 0 mm, 5 mm, 10 mm, 15 mm, 20 mm, and 25 mm away from the scanning center. The circular scanning radius was set to 50 mm. The simulated unfocused circular transducer had a diameter of 6 mm and a center frequency of 5 MHz with 80% fractional bandwidth. The grid size of the ROI (i.e., the reconstruction area) was 0.2 mm, the signal sampling rate was 20 MHz, and the speed of sound was set to 1500 m/sec. 720 scanning angles in 360-degree full view were used to mitigate the back-projection streak artifacts at the imaging points close to the scanning transducer. Note that there was no noise added in the simulations.
Figures 3(a) and 3(b) show the reconstructed PAT images (i.e., point spread functions) of the six simulated point absorbers with the model-based correction method and the back-projection algorithm, respectively. The right panels in Fig. 3 (a) and 3 (b) show the close-up reconstructed images of each point absorber. Figure 3(a) clearly shows all the six point absorbers – all with similar point spread functions whereas Fig. 3(b) fails to show the point shape of the absorbers except for the ones at and near the scanning center. It is evident in Fig. 3(b) that with the back-projection algorithm, the object is blurred and elongated in the tangential direction when the point absorber is far from the scanning center. Note that the sidelobe level in the model-based corrected images (Fig. 3(a)) becomes relatively higher in tangential direction typically when the point absorber falls within the near field of the transducer, which is 30 mm from the simulated 5-MHz transducer surface. Within the near field of the transducer, the detection angle is more limited as illustrated in Fig. 1, and hence fewer A-lines in the PAT sonogram can contribute to the model-based correction method, consequently ending with higher sidelobes. However, the sidelobe level in the model-based close-up images is still lower than that in the back-projection reconstructed ones.
Fig. 3 Reconstructed PAT images of simulated point absorbers where (x,y) = (0,0) is the scanning center. (a) Model-based correction method (b) Back-projection.
The tangential resolution, radial resolution, and peak amplitude of the six point spread functions in Figs. 3(a) and 3(b) as a function of the distance away from the scanning center are quantified and shown in Fig. 4 . Note that the resolution is defined as the full width half maximum of the projected tangential and radial profiles of each point spread function. Figure 4(a) shows that our proposed method significantly improves the tangential resolution – space-invariant tangential resolution is retrieved. Our proposed method provides four-fold improvement in tangential resolution over the back-projection algorithm when the point absorber is 10 mm away from the scanning center. Comparing Figs. 4(a) and 4(b), we can find that the model-based method provides the same tangential and radial resolution over all the imaging area. For both of the methods, the radial resolution, which is mainly determined by the transducer bandwidth, is almost space-invariant as described in reference [9]. In addition to the resolution improvement, as shown in Fig. 4(c), the proposed method also offers uniform peak amplitude of the point spread functions over all the imaging area because the diffraction effect of the unfocused transducer is compensated. This feature is important for quantitative and molecular PAT [15].
Fig. 4 (a) Tangential resolution, (b) radial resolution, and (c) peak amplitude as a function of the distance away from the scanning center
3.2 Experiments
Performance of our proposed method was then further demonstrated with experimental PAT data. The PAT system used here was similar to those in the references [7] and [10]. A porcine gel phantom embedded with five human black hairs (~70-μm in diameter) was made. Cross-sectional views of the five human black hairs were imaged to evaluate the point spread functions of the PAT system. Laser pulses at 800 nm were applied to excite the photoacoustic signals. The delivered laser energy was homogenized within a circular region with ~12-mm diameter and centered at the scanning center. The photoacoustic signals were detected by a 5-MHz unfocused circular transducer (V310, Panametrics-NDT, USA), collected at 240 positions along a complete circle, and digitized at a 25-MHz sampling rate with 8-bit resolution. The active area of the transducer was 6 mm in diameter. The circular scanning radius, grid size of the reconstructed ROI, and the speed of sound used in the both reconstruction methods were the same as those in the simulations. Note that the signals were averaged 40 times to improve the signal-to-noise ratio (SNR).
Figures 5(a) and 5(b) show the reconstructed cross-sections of the five hairs with the model-based method and the back-projection algorithm, respectively. Figure 5 is presented in the same format as Fig. 3. The weak image amplitude of the two hairs located at 16 mm and 22 mm away from the scanning center was caused by our limited light illumination area. The experimental results of the proposed method are generally in good agreement with the simulation ones. Unlike the simulations, the back-projection reconstructed image suffers the interferences from hair sections out of image plane, as indicated by the arrows in Fig. 5(b), which are more pronounced in the close-up images. Such out-of-image-plane interferences hamper the visibility of the cross-sections of the hairs in Fig. 5(b) while they are suppressed by the proposed method since our spatiotemporal optimal filter is designed to reconstruct the designated 2D image plane only. Note that currently, the proposed method takes 10 minutes to reconstruct Fig. 5(a) on a PC with a 3.33-GHz 6-core CPU and 4-GB memory.
Fig. 5 Reconstructed PAT images using experimental data where cross-sections of five human black hairs are imaged and the scanning center is at (x,y) = (0,0). (a) Model-based correction method (b) Back-projection. The arrows in (b) indicate the interferences from hair sections out of image plane.
4. Conclusions
The proposed method effectively improves the degraded tangential resolution for PAT with finite-sized unfocused transducers while retaining the radial resolution. Uniform image quality in terms of resolution and image intensity over all the imaging area is achievable, which is important for breast imaging, quantitative and molecular PAT. In addition, by taking the electrical impulse response of the transducer into account, the proposed method can potentially improve the radial resolution. It is worth noting that similar and even 3D models and spatiotemporal filter design based on other metrics, instead of the minimum mean square error sense used here, can also be developed for different photoacoustic imaging modes although we concentrate on the tomography mode in this study [16,17]. Future work will focus on studying the effect of SNR and tissue inhomogeneity of speed of sound on the proposed model-based method.
Acknowledgments
Funding supports from National Science Council, Taiwan (NSC 97-2221-E-007-084-MY3), the Electronics and Optoelectronics Research Laboratories, ITRI, Taiwan (9351AA5110), and National Tsing Hua University, Taiwan (Boost program 98N2531E1 and 99N2905E1) are greatly appreciated. The authors would also like to thank Prof. Ray-Kuang Lee, Institute of Photonics Technologies, National Tsing Hua University, Taiwan, for providing workstations to run our simulations in the preliminary study phase.
References and links
1. X. Wang, Y. Pang, G. Ku, X. Xie, G. Stoica, and L. V. Wang, “Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain,” Nat. Biotechnol. 21(7), 803–806 (2003). [CrossRef] [PubMed]
2. L. V. Wang, “Tutorial on photoacoustic microscopy and computed tomography,” IEEE J. Sel. Top. Quantum Electron. 14(1), 171–179 (2008). [CrossRef]
3. M.-L. Li, J. C. Wang, J. A. Schwartz, K. L. Gill-Sharp, G. Stoica, and L. V. Wang, “In-vivo photoacoustic microscopy of nanoshell extravasation from solid tumor vasculature,” J. Biomed. Opt. 14(1), 010507 (2009). [CrossRef] [PubMed]
4. J. T. Oh, M.-L. Li, H. F. Zhang, K. Maslov, G. Stoica, and L. V. Wang, “Three-dimensional imaging of melanoma skin cancer in-vivo by dual wavelength photoacoustic microscopy,” J. Biomed. Opt. 11, 034032 (2006). [CrossRef]
5. S. A. Ermilov, T. Khamapirad, A. Conjusteau, M. H. Leonard, R. Lacewell, K. Mehta, T. Miller, and A. A. Oraevsky, “Laser optoacoustic imaging system for detection of breast cancer,” J. Biomed. Opt. 14(2), 024007 (2009). [CrossRef] [PubMed]
6. M.-L. Li, J.-T. Oh, X. Y. Xie, G. Ku, W. Wang, C. Li, G. Lungu, G. Stoica, and L. V. Wang, “Simultaneous molecular and hypoxia imaging of brain tumors in vivo using spectroscopic photoacoustic tomography,” Proc. IEEE 96(3), 481–489 (2008). [CrossRef]
7. G. F. Lungu, M.-L. Li, X. Xie, L. V. Wang, and G. Stoica, “In vivo imaging and characterization of hypoxia-induced neovascularization and tumor invasion,” Int. J. Oncol. 30(1), 45–54 (2007).
8. F. Lingvall, “Time domain reconstruction methods for ultrasonic array imaging,” Ph.D. dissertation, Signals and Systems, Uppsala University, Uppsala, Sweden (2004).
9. M. Xu and L. V. Wang, “Analytic explanation of spatial resolution related to bandwidth and detector aperture size in thermoacoustic or photoacoustic reconstruction,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 67(5), 056605 (2003). [CrossRef] [PubMed]
10. M. Pramanik, G. Ku, and L. V. Wang, “Tangential resolution improvement in thermoacoustic and photoacoustic tomography using a negative acoustic lens,” J. Biomed. Opt. 14(2), 024028 (2009). [CrossRef] [PubMed]
11. M.-L. Li and L. V. Wang, “A study of reconstruction in photoacoustic tomography with a focused transducer,” Proc. SPIE 6437, 64371E (2007). [CrossRef]
12. P. R. Stepanishen, “Transient radiation from pistons in an infinite planar baffle,” J. Acoust. Soc. Am. 49(5B), 1629–1638 (1971). [CrossRef]
13. B. Piwakowski and B. Delannoy, “Method for computing spatial pulse response: Time-domain approach,” J. Acoust. Soc. Am. 86(6), 2422–2432 (1989). [CrossRef]
14. F. Lingvall, T. Olofsson, and T. Stepinski, “Synthetic aperture imaging using sources with finite aperture: deconvolution of the spatial impulse response,” J. Acoust. Soc. Am. 114(1), 225–234 (2003). [CrossRef] [PubMed]
15. A. Rosenthal, D. Razansky, and V. Ntziachristos, “Fast semi-analytical model-based acoustic inversion for quantitative optoacoustic tomography,” IEEE Trans. Med. Imaging 29(6), 1275–1285 (2010). [CrossRef] [PubMed]
16. M. Roumeliotis, P. Ephrat, J. Patrick, and J. J. L. Carson, “Development and characterization of an omnidirectional photoacoustic point source for calibration of a staring 3D photoacoustic imaging system,” Opt. Express 17(17), 15228–15238 (2009). [CrossRef] [PubMed]
17. M. Roumeliotis, R. Z. Stodilka, M. A. Anastasio, G. Chaudhary, H. Al-Aabed, E. Ng, A. Immucci, and J. J. L. Carson, “Analysis of a photoacoustic imaging system by the crosstalk matrix and singular value decomposition,” Opt. Express 18(11), 11406–11417 (2010). [CrossRef] [PubMed]
