Open Access
Add to BibSonomyAbstract
Photoacoustic tomography is a promising and rapidly developed methodology of biomedical imaging. It confronts an increasing urgent problem to reconstruct the image from weak and noisy photoacoustic signals, owing to its high benefit in extending the imaging depth and decreasing the dose of laser exposure. Based on the time-domain characteristics of photoacoustic signals, a pulse decomposition algorithm is proposed to reconstruct a photoacoustic image from signals with low signal-to-noise ratio. In this method, a photoacoustic signal is decomposed as the weighted summation of a set of pulses in the time-domain. Images are reconstructed from the weight factors, which are directly related to the optical absorption coefficient. Both simulation and experiment are conducted to test the performance of the method. Numerical simulations show that when the signal-to-noise ratio is −4 dB, the proposed method decreases the reconstruction error to about 17%, in comparison with the conventional back-projection method. Moreover, it can produce acceptable images even when the signal-to-noise ratio is decreased to −10 dB. Experiments show that, when the laser influence level is low, the proposed method achieves a relatively clean image of a hair phantom with some well preserved pattern details. The proposed method demonstrates imaging potential of photoacoustic tomography in expanding applications.
© 2015 Optical Society of America
1. Introduction
Photoacoustic tomography (PAT), also referred to as thermoacoustic tomography or optoacoustic tomography, is an emerging methodology of biomedical imaging, combining high spatial resolution of acoustic measurement in deep tissue and rich optical contrast [1–4 ]. A conventional PAT system projects a laser/microwave beam on the region of interest (ROI). Broadband acoustic waves, i.e., photoacoustic (PA) signals are excited and recorded by the detectors placed around the ROI. Finally, a reconstruction algorithm [1,4–6 ] is employed to map the absorption coefficient distribution from the detected PA signals. Owing to the capability of structural, functional and molecular imaging in vivo at a new path, PAT has shown its potentials in mapping human breast [6], brain [7], gene expression [8], sentinel lymph nodes [1], etc.
With the expanding of PAT applications, it has become an increasing urgent problem to reconstruct image from weak and noisy PA signals. On the one hand, it is highly beneficial to extend the imaging depth to the hard limit for optical penetration [2,4,9], where the optical energy is low and the quantity of diffuse photons is small. In the case of deep penetration, less illumination energy can be delivered into absorbers and then the induced acoustic pressure is relatively weak. Furthermore, the pressure will experience more decay in intensity before recorded by detectors due to increasing distance between absorbers and detectors. Consequently, PA signals from the deep tissue will be weak and the signal-to-noise ratio (SNR) is low. Therefore, the image reconstruction method focusing on the low SNR conditions is especially useful to improve imaging quality of PAT in the deep tissue. On the other hand, the ability of image reconstruction from weak and noisy signal could decrease the dose of laser exposure used for PAT. Low illumination energy guarantee that PAT could be safe enough to be expanded to biomedical applications where the maximum permissible exposure (MPE) [9–11 ] is low, e.g., ocular imaging [12]. In short, image reconstruction from weak and noisy signals has significant value in increasing the imaging depth and decreasing the exposure dose.
However, since the frequency nature of nanosecond pulse laser [1,7,13–15 ], the induced PA signals usually have an wide frequency-band, and such a wide-band signal is easy masked by broadband noise when SNR is low. As a result, the conventional methods are often inefficient to recover images from weak and noisy signals [14,16].
In order to overcome this limitation, we propose an algorithm based on pulse decomposition in the time-domain. In this method, PA signals are decomposed as a set of weighted pulses according to the impulse response of the imaging system. Then the image is reconstructed from the weight factors, which is directly related to the absorption coefficient. Numerical simulations are conducted to compare the proposed method with the conventional back projection method [5,17] under a weak and noisy signal conditions. Finally, hair pattern imaging experiments are carried out to validate the proposed reconstruction method in practical situations.
2. The pulse decomposition method and its implementation
In response to an illumination, a PA pulse, p(r,t), at position r and time t obeys the following wave equation [4],
where c is the speed of sound, Γ is the Gruneisen parameter, and Ie(t) is the temporal profile of the illumination function, A(r) is the absorbed energy density and is related to the optical absorption coefficient µa(r) by A(r) = µa(r)·H′(r), H′(r) is the radiant exposure. Conventionally, the laser is projected evenly on the absorber, and then µa(r) is proportional to A(r), which is approximately regarded as the optical absorption coefficient in the paper to reconstruct the absorber distribution.Solving Eq. (1) we have PA fields,
where r = |r – r′|, and * is used to denote convolution. Given that the real ultrasound detection system has a limit bandwidth, the PA signals recorded by detectors can be expressed as [5], where Id is the impulse response of the ultrasound detection system. The task of PAT image reconstruction is to approximately recover A(r) from pd. In order to solving the inverse problem of image reconstruction under the low SNR conditions, a pulse decomposition algorithm (PDA) in the time-domain is proposed as follows.Introduce the solid angle Ω to rewrite the triple integral in Eq. (2), and ∫r∫∫ dΩdr = ∫∫∫ d 3 r. Substituting Eq. (2) into Eq. (3) obtains the form of the spherical coordinate system with the origin point at r
where . Integrating items in the square brackets converts p as a convolutionChanging Eq. (5) in an integral form, we have
The physical meaning of Φ in Eq. (5) is easily explained as shown in Fig. 1(a). Absorbers located at the same ring can be considered as a group. It means that, considering A(r) is proportional to µa(r), [tΦ(r,t)] is proportional to the sum of the optical absorption coefficients of absorbers located at the subregion, which is a thin ring with a center of r and a radius of [tc]. According to Eq. (5), as long as the thickness of the ring is small enough, a PA signal can be decomposed as a set of pulses with the weight factors of Φ(r,t). As shown in Fig. 1(a), since all absorbers are located at the three rings, the PA signal is the sum of the three pulses with the weight factors of Φ(r,ti− 1), Φ(r,ti) and Φ(r, ti +1).
Fig. 1 (a)Absorbers in the ROI are distributed in three thin rings with the center of r. ta and tb are the propagation time from r to the nearest and furthest rings in the ROI, respectively. The PA signal recorded by the detector is indeed the sum of the three pulses with weight factors Φ(r, ti −1), Φ(r, ti) and Φ(r, ti +1), respectively. (b)There are M detectors placed circularly and evenly around the ROI. Any detector, such as r α or r β, records a PA signal, which will be decomposed to obtain Φ. And Φ is proportional to the integral of the optical absorption coefficient over the arcs with the center of r α or r β. Taking the far field conditions into account, the arcs within the ROI can be approximately regarded as straight lines.
Next, a discrete formula will be provided to calculate Φ in the time-domain, by decomposing pd as a set of PA pulses. The discrete form of Eq. (6) will be expressed simply as follows
where pd(n) is the discrete expression of pd(r,nts), t = nts, n is discrete time and ts is the sampling interval, na and nb are the lower and upper boundaries of , as the discrete forms of ta and tb shown in Fig. 1(a).Eq. (7) can be presented in the style of matrix
With the known H and pd, Eq. (8) can be solved to obtain the solution of Φ by using the least square method [3,18].
After repeating the above decomposition, every PA signal recorded by different detectors will be decomposed as the sum of pulses, and the weight factors of these pulses are proportional to Φ. As shown in Fig. 1(b), if there are M detectors placed circularly and evenly around the ROI, we can get , where r m is the position of the m-th detector, and m ∈ [1,M]. There are direct and strong correlation between A and Φ. From the perspective of geometry, is proportional to the integral of A(r) over the arc with the center of r m, and the radius of |r − r m|, as shown in the left part of Fig. 1(b). Furthermore, the arcs can be approximately considered as the straight line as shown in the exactly enlarged graph at the right part of Fig. 1(b), because the distance between the detector and the ROI is much greater than the physical size of the ROI. Naturally, the reconstruction theory of the parallel projection widely employed in the X-ray tomography [19,20] is available and yields the total formula of the PDA as follows
3. Numerical comparisons
Numerical simulations are conducted to validate performance of the PDA under the low SNR conditions. The measurement geometry is shown in Fig. 2. The ROI is a box, of which the height is HR, and both the length and width are LR. The detector array, surrounding the ROI, consists of M acoustic detectors, distributed circularly and evenly. Each detector has a width of Ld, and its surface is an inward arc with the height Hd and the radius Rs. The center points of all detectors are at the same height as that of the center point of the ROI, and the distance between them is Rd.
Fig. 2 The ROI is a box (a) with the length LR, the width LR and the height HR, which is surrounded by the detectors distributed circularly and evenly. Rd is the distance between the center point of the ROI and the center point of the detector surface (b), and the surface of any detector is an inward arc with the radius Rs, the height Hd and the width Ld.
As shown in Fig. 3(a), there are three absorbers (a cuboid, a cylinder and a triangular prism) in the ROI, with the same absorbed energy density (i.e., the same optical absorption coefficient) and the same height of HR. In the simulations, the speed of sound is c = 1455 m/s, the sampling frequency is fs= 60 MHz, and LR = 6 mm, HR =0.5 mm, Rd= 50 mm, M = 120, Rs= 50 mm, Hd =6 mm, Ld=6 mm.
Fig. 3 The true distribution (a) of absorption coefficients in the ROI contains a rectangular, a circle and a triangular. When SNR is 40 dB, the reconstructed image (b) and (c) are produced by the BP method and the PDA, respectively. When SNR decreases to −2 dB, clean and noisy PA wave are compared in the time-domain (d) and the amplitude spectrum (e). Moreover, for the pixel profile located at the straight line shown in (a), its true value and reconstructed ones by the PDA and BP method with SNR changing from 40 dB to −2 dB are shown in (f), respectively.
In order to simulate the response characteristics of the ultrasound system, the generated PA signals are filtered by a band-pass filter with the low cutoff frequency of 5 MHz and the high cutoff frequency of 15 MHz. The waveform and spectrum of a typical PA signal are given in Figs. 3(d) and 3(e), respectively.
The PDA is proposed to recover the image from the PA signals. Besides, for the sake of comparison, the back projection (BP) method [5], widely used for PAT reconstruction, is employed to recover the image from the same signals. When PA signals are clean, i.e., the signal-to-ratio (SNR) is high (over 20 dB), both methods reconstruct almost the same images and recover the true distribution well with smearing artifacts, which would tend to be eliminated if the frequency band of detectors is wide enough and the detector array is spherical in three dimensions. As shown in both Figs. 3(b) and 3(c), the absorbers can be recognized easily with correct shape and size.
However, when the PA signal become weak and noisy (as shown in Fig. 3(d)), there are strong disturbances at the amplitude spectrum (as shown in Fig. 3(e)). Furthermore, the profile comparison is given in Fig. 3(f) to demonstrate the noise influence directly, in which the profile ”BP(−2dB)” has more obvious deviation from the profile ”true”. To evaluate the impact of noise on PAT imaging, white noise with various SNR is generated and then mixed into single PA signal one by one to simulate noisy signals, where , Pclean is the power of a PA signal, and Pnoise is the power of the added white noise. As analyzed before, the bandwidth of pd is large, and it is hard to depress noise if Eq. (5) is solved by multiplication over such a broad frequency band, and then the noise mixed in A(r) will gradually become the main cause of serious artifacts in the reconstructed image (as shown in Figs. 4(e), 4(g) and 4(i)).
Fig. 4 when SNR decreases from 8 dB to −6 dB, the reconstructed images of the BP method in the first row become unrecognizable gradually, while the images of the PDA in the second row keep recognizable. Additionally, even though SNR decreases from −6 dB to −10 dB further, the images of the PDA still have a certain reference value.
Furthermore, as shown in Fig. 4 with decreasing SNR, there are more and more blur in examples of the images reconstructed by the BP method (as shown in the first row of Fig. 4), and they become unrecognizable when SNR is less than −2 dB. However, the PDA can produce relatively sharp images (as shown in the second row of Fig. 4) even when SNR is −10 dB.
To display reconstruction performance more directly and statistically, a comparison chart, as shown in Fig. 5, exhibits the mean square error (MSE) of two methods, where error is defined as pixel value difference between a reconstructed image and the true one. Additionally, each mean-value dot shows statistical error of 20 repetition simulations with the same SNR.
Fig. 5 A mean plot of MSE vs. SNR is given. The red and blue dots exhibit statistical performance of the BP method and the PDA, respectively. And any dot represents the mean MSE value of 20 repetition simulations with the same SNR.
As shown in Fig. 5, when the PA signals are clean or SNR is increasing, the MSE of the PDA is approaching to zero as the BP method. However, when SNR is decreasing, The PDA performs better obviously. For example, when the SNR is −4 dB, the mean MSE of the BP method and the PDA are 0.0743 and 0.0126, respectively. It indicates that the PDA decreases the MSE to about 17% of the conventional reconstruction method. Furthermore, the mean MSE of the PDA is 0.0283 when SNR is −10 dB, while the mean MSE of the FPB method is 0.0317 when SNR is 4dB. It implies that, if we can accept an image reconstructed by the BP method with a not bad SNR (4 dB), then it is easier for us to accept the image reconstructed by the PDA even though the signals are much more weak and noisy (SNR is −10 dB).
4. Experiments
Finally, experiments are employed to examine the performance of the PDA in practices. The schematic diagram of the experimental setup is shown in Fig. 6(a). The illumination function is produced by a Q-switched Nd:YAG laser with a wavelength of 532 nm and an energy of 80mJ. An agar phantom with a diameter of about 35 mm, as shown in Fig. 6(b), is vertically fixed on the base and sunk in water. A laser beam is projected on the upper surface of the phantom, and the induced PA waves are recorded by an ultrasound detector (V310, Panametrics) with the center frequency of 4.39 MHz and a bandwidth of 100.2% at −6dB. The detector is driven by a computer controlled stepper to scan circularly and evenly. For the rotating detector, the center point of its surface is always kept at the same height as that of the hair pattern ”AI”. PA waves are amplified (SA-230F5,NF) and received by the signal acquisition card (PCI-5105,NI) with a sampling frequency of 60MHz. To test the reconstruction performances with various SNR, we change the laser influence level and get two sets of PA signals.
Fig. 6 (a) Schematic diagram of the PAT imaging system to map an agar phantom surrounded by water. (b) the picture of the agar phantom with hairs placed at the same height and constituting the pattern of ”AI”. (c) the reconstructed image of the PDA when the laser influence level is high. (d) the reconstructed image of the BP method when the laser influence level is high. (e) the reconstructed image of the PDA when the laser influence level is low. (f) the reconstructed image of the BP method when the laser influence level is low.
As shown in Figs. 6(c) and 6(d), when the laser influence level is high, both the BP method and PDA can provide acceptable images, i.e., the pattern is recognizable with correct shapes and sizes. Nevertheless, when the laser influence level is low, as shown in Figs. 6(e) and 6(f), the image recovered by the BP method is totally polluted and even the basic outline of the pattern is masked by the noise, while the PDA still provides a relatively clean image, where the noise is depressed well and some details about the pattern are maintained, such as in the character ”A”, the right downward line is a little longer than the left downward one.
It is noteworthy that, in the process of recovering images from pd, the PDA in the time-domain is selected, other than the conventional methods, mainly because the former achieves better anti-noise performance with acceptable computational cost. On the one hand, caused by the illumination function lasting at the nanosecond level, the energy of PA signals is spread over a wide frequency band. Consequently, if we solve Eq. (6) in the frequency-domain conventionally, a lot of noise will be brought into the spectrum of A(r), and then give rise to serious artifacts in the reconstructed image. On the other hand, owing to the short-pulsed illumination too, most elements in H(n) are zeroes. Therefore, sparse matrix techniques can be adopted to solve Eq. (8), and then reduce the amount of computation significantly. For example, to reconstruct the Figs. 3(b) and 3(c), a Windows Workstation with an Intel Core i7-2600 3.40 GHz processor having 16 GB memory is used in all computations carried out in this work. In the identical reconstruction task, the computing time required by the BP method and PDA are 36.0475 s and 143.3621 s, respectively. It shows that the PDA increases the reconstruction time by 300% approximately.
5. Conclusions
With expanding of PAT applications, reconstructing image from weak and noisy signals deserves more attention since it is highly beneficial to extend the imaging depth and able to decrease the dose of laser exposure. Usually, the frequency bands of PA signals are so wide that the signals are easily masked by noise. Therefore, a pulse decomposition algorithm in the time-domain is proposed to achieve better anti-noise performance. In the algorithm, PA signals are decomposed as a set of PA pulses, whose weight factors are directly related to the optical absorption coefficient.
Simulations and experiments are conducted to test the performance of the proposed algorithm in comparison with the conventional back projection method. They prove that, at the cost of about quadrupled reconstruction time, the proposed algorithm performs better when SNR is low. For example, the simulations show that, when SNR is −4 dB, the proposed method decreases the reconstruction error to about 17%; the experiments show that, even though the laser influence level is decreased much and consequently the PA signals are weak, the proposed method still reconstructs a relatively clean image of the hair phantom with some well-preserved pattern details. The pulse decomposition algorithm demonstrates PAT is acceptable under weak and noisy signal conditions.
Acknowledgments
Project supported by the National Basic Research Program of China (Grant No. 2012CB921504), the National Natural Science Foundation of China (Grant Nos. 11422439, 11274167, 11274171, 61201450, 61201495, 61302175), and Chongqing Science and Technology Commission of China(Grant Nos. 2012jjA40058 and 2012jjA40006).
References and links
1. C. Li and L. V. Wang, “Photoacoustic tomography and sensing in biomedicine,” Phys. Med. Biol. 54, R59–R97 (2009). [CrossRef] [PubMed]
2. L. V. Wang, “Prospects of photoacoustic tomography,” Med. Phys. 35, 5758–5767 (2008). [CrossRef]
3. J. Rudzki-Small, L. J. Libertini, and E. W. Small, “Analysis of photoacoustic waveforms using the nonlinear least square method,” Biophys. Chem. 41, 29–48 (1992). [CrossRef]
4. M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum. 77, 1–22 (2006). [CrossRef]
5. M. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E 71, 1–7 (2005). [CrossRef]
6. C. Tao and X. Liu, “Reconstruction of high quality photoacoustic tomography with a limited-view scanning,” Opt. Express 18, 2760–2766 (2010). [CrossRef] [PubMed]
7. J. Xia and L. V. Wang, “Small-animal whole-body photoacoustic tomography: A review,” IEEE Trans. Biomed. Eng. 61, 1380–1389 (2014). [CrossRef]
8. B. J. Vakoc, D. Fukumura, R. K. Jain, and B. E. Bouma, “Cancer imaging by optical coherence tomography: preclinical progress and clinical potential,” Nat. Rev. Cancer 12, 363–368 (2012). [CrossRef] [PubMed]
9. G. Ku and L. V. Wang, “Deeply penetrating photoacoustic tomography in biological tissues enhanced with an optical contrast agent,” Opt. Lett. 30, 507–509 (2005). [CrossRef] [PubMed]
10. M. Jeon, S. Jenkins, J. Oh, J. Kim, T. Peterson, J. Chen, and C. Kim, “Nonionizing photoacoustic cystography with near-infrared absorbing gold nanostructures as optical-opaque tracers,” Nanomed 9, 1–10 (2013).
11. R. a. Kruger, R. B. Lam, D. R. Reinecke, S. P. Del Rio, and R. P. Doyle, “Photoacoustic angiography of the breast,” Med. Phys. 37, 6096–6100 (2010). [CrossRef] [PubMed]
12. N. Wu, S. Ye, Q. Ren, and C. Li, “High-resolution dual-modality photoacoustic ocular imaging,” Opt. Lett. 39, 2451–2454 (2014). [CrossRef] [PubMed]
13. R. G. M. Kolkman, W. Steenbergen, and T. G. van Leeuwen, “In vivo photoacoustic imaging of blood vessels with a pulsed laser diode,” Lasers Med. Sci. 21, 134–139 (2006). [CrossRef] [PubMed]
14. Y. Yang, S. Wang, C. Tao, X. Wang, and X. Liu, “Photoacoustic tomography of tissue subwavelength microstructure with a narrowband and low frequency system,” Appl. Phys. Lett. 101, 034105 (2012). [CrossRef]
15. L. Liu, C. Tao, X. Liu, X. Li, and H. Zhang, “Pulse decomposition based analysis of PAT/TAT error caused by negative lobes in limited view conditions,” Chinese Phys. B 21, 024304 (2015). [CrossRef]
16. X. Gao, C. Tao, X. Wang, and X. Liu, “Quantitative imaging of microvasculature in deep tissue with a spectrum-based photo-acoustic microscopy,” Opt. Lett. 40, 970–973 (2015). [CrossRef] [PubMed]
17. P. Burgholzer, J. Bauer-Marschallinger, H. Grün, M. Haltmeier, and G. Paltauf, “Temporal back-projection algorithms for photoacoustic tomography with integrating line detectors,” Inverse Prob. 23, S65–S80 (2007). [CrossRef]
18. Y. Chao, “An implicit leastsquare method for the inverse problem of acoustic radiation,” J. Acoust. Soc. Am. 81, 1288–1292 (1987). [CrossRef]
19. G. T. Herman, A. V. Lakshminarayanan, and A. Naparstek, “Convolution reconstruction techniques for divergent beams,” Comput. Biol. Med. 6, 259–271 (1976). [CrossRef] [PubMed]
20. A. C. Kak, “Computerized tomography with X-ray, emission, and ultrasound sources,” Proc. IEEE 67, 1245–1272 (1979). [CrossRef]
