Open Access
Add to BibSonomyAbstract
The goal of this work is to resolve the limited-view problem of photoacoustic tomography (PAT). We report a two-loop iteration method to inverse the photoacoustic sources from the measured photoacoustic signals. PAT reconstruction with this method does not depend on the detection path. Therefore, the proposed method can provide recognizable image even when the detector only scans a small angle (about 20°~30°). The comparison with the delay-and-sum method shows the advantage of the proposed method in reconstructing image from incomplete data.
©2010 Optical Society of America
1. Introduction
Photoacoustic tomography (PAT) combines the merits both of pure ultrasound imaging and pure optical imaging [1–3]. Ultrasonic imaging can provide a better spatial resolution than optical imaging in opaque media with depths greater than ~1mm, because ultrasound scattering is two to three orders of magnitude weaker than optical scattering in biological tissues. On the other hand, optical imaging can provide high contrast functional imaging, referred to oxygen saturation and the concentration of hemoglobin in biological tissue. PAT yields images with simultaneously high contrast and high spatial resolution [2,3]. Moreover, the nonionizing wave used to excite photoacoustic effects is much safer for biological tissue than ionizing wave (e. g., x-ray radiation). These merits make PAT become a promising modality for the non-invasive imaging of animal or human organs [4–6], such as breast cancer detection [7], skin lesions examination [8], blood oxygenation study [9], brain imaging [10], molecular imaging [11,12], coronary artery stents imaging [13], cardiovascular dynamics and microcirculation study [14,15], moving target imaging [16] etc.
In many case of these studies, it was assumed the full view detection, i.e., the detector scans along a whole circle enclosing the objects. Incomplete scanning could blur the image or even make the reconstruction impossible [17,18]. However, inhomogenous complex structures in biological organs could always prevent the collection of the signals from all directions, e.g., the solid angle of detection is at most 180° for a breast. Therefore, ones have to face the limited-view problem in the practical application of PAT [17].
Jaeger et. al. proposed a pseudo-inverse method to inverse the linear projection between the initial energy distribution and the detected signal [19]. Based on this concept, they optimized the Fourier algorithm [20] and the iterative method [21]. The studies show that the pseudo-inverse method efficiently reduces the artifacts due to incomplete data [19–22]. However, the performance of this pseudo-inverse method for a narrow scanning angle has not been well examined.
In this study we report a two-loop iteration (TLI) reconstruction method to resolve the limited-view PAT problem. A numerical experiment based on finite-element analysis was used to examine the performance the TLI method in reconstructing PAT image from incomplete data. The delay-and-sum (DAS) method [23] was also employed to reconstruct the PAT image and compared with the TLI method.
2. Methods
PAT reconstruction can be considered to be an inverse problem [19]. It is known that there exists a linear projection between the initial optical energy distribution and the detected photoacoustic waves. Inversing this linear projection, the initial optical energy distribution can be obtained from the detected acoustic waves. The method proposed in this study follows this basic idea. Let us consider the acoustic wave p(r, t) prompted by the pulse Ie(t) [24,25]:
where c is the velocity of sound and p 0(r) describes the initial optical energy distribution. It has p(r, t) = (4πc 2)−1∫∫∫| r – r '|≤ct p 0(r')τ(t – |r – r'|/c)|r – r'|−1dV'. For a 2D tomography r = (x, y), the studied area is discretized into a grid of size Nx × Ny, as shown in Fig. 1 . The area of each element is ΔxΔy, where Δx and Δy are the element widths in the x and y directions, respectively. The initial optical energy distribution is described by the vector ξ = [ξ 1, …, ξj, …, ξJ]T, where J = Nx × Ny and ξj is the average optical energy p 0 in one element at position r j = (xnx, yny), j = (nx – 1)Ny + ny (nx = 1, …, Nx; nx = 1, …, Ny). The photoacoustic information detected at r k = (xk, yk) is stored in the vector Φ = [φ 1, …, φi, …, φI]T, where I = Nd × Nt and φi = ∫t ≤ l Δ t p(xk, yk, t)dt with i = (k – 1)Nt + l (k = 1, …, Nd; l = 1, …, Nt), Δt is the sampling interval in time domain, Nt and Nd are the numbers of the temporal and spatial sampling points. Then the linear projection between ξ and Φ is written as a matrix E and [18,19]:The element Ei j of E can be integrated by
where H = ψ -1/2 for ψ ≥ 0 and H = 0 for ψ < 0 [ψ = c2l2∆t2 – (x – xk)2 – (y – yk)2].With the measured Φ and the known matrix E, the spatial distribution of the optical energy p 0 can be obtained for imaging by solving the linear Eq. (2). For the PAT problem, p 0 is always physically nonnegative. In other words, the solution ξ of Eq. (2) should be nonnegative. Therefore, a two-loop iterative process is employed to estimate the nonnegative least-squares solution of Eq. (2) [26]. The details of this algorithm can be found in ref [26]. Here, we briefly introduce this algorithm as follows:
- (1) The iterative process begins with initializing ξ ← {0, …, 0}T and two sets P ← Ø and Q ← {1, 2,…, J}.
- (2) the J-vector w ← E T(Φ – Eξ) is calculated, where w = {w 1, …, wj, …, wJ}T.
- (3) Find an index q∈Q such that wq = max{wj: j∈Q}. And move the index q from set Q to P, i.e., Q ← Q–{q} and P←P∪{q}.
The following steps (4-7) are the sub-loop,
- (4) The Beginning of the sub-loop: Define a new matrix F, whose j-th column Fj is the j-th column of E if j∈P, otherwise, Fj ← {0,…,0}T,
- (5) Calculate the least-square solution ζ = (FTF)−1FTΦ of the linear equation F × ζ = Φ, where ζ = {ζ1, …, ζj, …, ζJ}T.
- (6) If there exists any ζj ≤ 0 for j∈P, ζ vector is iterated by ζ ←ζ–α(ζ–ξ), where α is the minimum of ξj/(ξj –ζj) with ζj ≤ 0 and j∈P. Otherwise (ζj > 0 for ∀j∈P), go to step 8 and the sub-loop iteration is completed.
- (7) The end of the sub-loop: Move from set P to set Q all indices j∈P for which ξj = 0. Repeat this sub-loop from step 4.
- (8) Set ξ←ζ. One iterative cycle of the main loop is completed and a new cycle begins again from the step (2).
- (9) Finally, the whole iterative process is finished until the set Q is empty or the vector w calculated in step 2 satisfies wj ≤ 0 for all j∈Q.
In the above description, the left arrow ‘A←B’ denotes that the value B is assigned to the variable A. The steps 2-9 are the main loop and the steps 4-7 are the sub-loop. Therefore, this is a two-loop algorithm.
The basic idea of our reconstruction is to solve the inverse problem described by the linear projection between the initial optical energy distribution and the detected photoacoustic signal. The effectiveness of this basic idea has been examined by many previous works [18–22]. Different from them, e.g. the pseudo-inverse method and the iterative method, we adopt a two-loop algorithm to directly get the nonnegative least-square solution of the inverse problem. It has been proved rigorously in mathematics that this two-loop algorithm is convergent and the solution ξ obtained by the two-loop iterative process minimizes ||Eξ – Φ|| subject to ξj ≥ 0 [26]. Therefore, the proposed reconstruction method is effective.
In addition, reviewing the above imaging method, it is noticeable that the relationship between the source p 0 and the detected signal p [Eq. (2)] is always true. The validity of Eq. (2) does not depend on the scanning orbit. If enough signals are detected (I > J), the inverse problem described by Eq. (2) could always be possible and the reconstruction by directly solving Eq. (2) could not depend on the scanning orbit too. Moreover, the nonnegative least-square solution obtained by the two-loop algorithm agrees with the physical picture of PAT, where the source is always nonnegative. It avoids the physically unreasonable solution due to the error of reconstruction with incomplete data. Therefore, the above new reconstruction method could provide a good PAT image with a limit-view scanning.
3. Numerical experiment and discussion
Numerical experiment was used to validate the proposed method. Figure 1 is the diagram of our numerical experiment. The finite-element method was applied to simulate the ultrasonic propagation in biological tissue with the density ρ = 1000 kg/m3 and c = 1500 m/s, where the ultrasound was excited by a rectangular laser pulse with an impulse width 5 × 10−7 s. There are two circular absorption sources in the studied area. Their radius are 2.0 mm and their centers are located at (x, y) = (0.0, 4.0) mm and (–4.0, –3.0) mm, respectively. A detector with a diameter 2.2 mm scanned the objects along the line from (10.0, –L/2) to (10.0, L/2) with a spatial sampling interval of 0.2 mm and a temporal sampling interval 10−7 s. The total imaging area has a size of 20 × 20 mm2 and was divided into a 50 × 50 grid. The range of the limited-view is quantified by the scanning angle θ = 2 × atan(L/2).
In the following discussion, for the sake of comparison, the weighted delay-and-sum (DAS) method [23] was also employed to reconstruct the PAT image. In this algorithm, the image value S(r) at the position r is determined as S(r) = ∑kw(k,r)p[t + δ(k, r)]/∑kw(k,r). Here, w(k,r) is the weighting factors, which takes into account the angular sensitivity of detector k for signals generated at the position r. p(t) is the photoacoustic signal, and δ(k,r) is the delay applied for detector element k and the source at r.
Figure 2 presents the reconstructed image with a limited-view scanning. When the scanning angle is relatively wide θ = 84°, the TLI method provides a high contrast image. Both the position and the shape of the sources can be clearly identified from the image [Fig. 2(a)]. However, in the same situation, the objects in the image reconstructed by the DAS method are dilated along the path around the center of the detection trace, which blurs the reconstructed image. This kind of image spread also occurs in the image reconstructed by the back projection method [17] and the iterative method [18]. The blurring due to the image spread decreases the contrast of the reconstructed image, moreover, makes it difficult to distinguish the two separated sources in the reconstructed image [Fig. 2(b)].
Fig. 2 The reconstructed image with a limited-view scanning. The images of the upper (a, c, e, g) and lower (b, d, f, h) rows are reconstructed by the TLI method and DAS method, respectively. The columns correspond to four scanning angle θ: (a, b) θ = 84°; (c, d) θ = 53°; (e, f) θ = 22°; (g, h) θ = 2.3°.
When the scanning angle is decreased to θ = 53° and 22°, the TLI method can still provide high contrast images. Although the shape of the sources is slightly distorted in these reconstructed images, the two sources can be undoubtedly recognized [Fig. 2(c) and 2(e)]. However, for the same scanning range, both the position and the shape of the sources cannot be clearly distinguished from the image reconstructed by the DAS method. The sources are distorted as concentric circles around the center of the detection trace [Fig. 2(d) and 2(f)]. When the detector only scanned an extreme narrow region θ = 2.3°, neither the TLI nor the DAS method can provide a recognizable image [Fig. 2(g), 2(h)].
Figure 3 quantifies the performance of the proposed method at different scanning angles. Figure 3(a) illustrates the full width at half-maximum (FWHM) of two sources. Because the incomplete data could cause the distortion of the reconstructed images, we use the FWHM along the x and y direction to quantify this characteristics. When θ > 60°, all FWHMs are close to the diameter 4mm of the sources. When θ is about 20°~60°, FWHMs along the x direction are still close to 4mm, FWHMs along the y direction are decreased. The two sources can still be distinguished, but their shape are compressed in the y-direction, as shown in Fig. 3(b) and Fig. 2(e). When θ < 20°, all FWHMs are much smaller than 4.0mm and the images of the sources are broken to many small pieces. Figure 3(c) gives the relative error inside the source 1 and 2. The relative error is defined as (∫∫|ps – pr|2 dxdy/∫∫ps 2 dxdy)1/2, where p s and pr are the normalized source distribution p 0 and the normalized recovered image. Figure 3(d) presents the average signal strengths (∫∫prdxdy/4π) of the sources 1 and 2 in the recovered image. Because the detector with a limited-view scanning can collect more photoacoustic signal from the close object (the source 1) than that from the distant object (the source 2), the signal strength of the source 1 is stronger than the source 2. This can also explain why the relative error of the source 2 became larger than that of the source 1 when θ is decreased to about 30°. Figure 3(e) proposes the correlation coefficient to assess the whole image quality. It is found that there is a significant linear correlation (the correlation coefficient > 0.75 for θ > 20°) between the recovered image and the initial optical energy distribution. All these measurements quantitatively show that the proposed TLI method could provide a distinguishable image, even when the detector only scanned narrow region (about 20°~30°) around the imaging region.
Fig. 3 Quantitative measurements of the reconstructed image quality as a function of the scanning angle: (a) the full width at half-maximum (FWHM), where xFWHM-n and yFWHW-n denote the FWHMs of the source n (n = 1, 2) along the x and y direction, respectively. (b) The 1D profile of source 2 with the scanning angle 33°. (c) the relative errors inside the sources 1 and 2, (d) the average signal strengths of the sources 1 and 2 in the recovered image (e) the correlation coefficients of the whole reconstructed image.
4. Conclusion
In summary, a new PAT reconstruction method is proposed in this study. In this method, a two-loop iterative method is used to directly estimate the initial optical energy distribution from its generated ultrasonic field. The reconstruction using this method is possible if only enough acoustic signals are collected, but not depended on the detection path. Moreover, the nonnegative least-square solution obtained by the two-loop algorithm agrees with the physical picture of PAT and avoids the physically unreasonable solution. Therefore, it can be a promising candidate for resolving the limited-view PAT problem. Our numerical experiments show that the proposed PAT method can provide acceptable images even though incomplete acoustic signals are collected only from a narrow view (about 20°~30°). Finally, it is worth to be mentioned that the reconstruction time for a 50 × 50 grid using the proposed method is about 18 minutes, which is much longer than the other method, such as the pseudo-inverse method [20]. This is one limitation of the proposed method.
Acknowledgments
This work was supported by NSF of China (No. 10874088 and No. 10904069) and the SRFDP200802840032.
References and links
1. T. Jetzfellner, D. Razansky, A. Rosenthal, R. Schulz, K.-H. Englmeier, and V. Ntziachristos, “Performance of iterative optoacoustic tomography with experimental data,” Appl. Phys. Lett. 95(1), 013703 (2009). [CrossRef]
2. F. Kong, Y. C. Chen, H. O. Lloyd, R. H. Silverman, H. H. Kim, J. M. Cannata, and K. K. Shung, “High-resolution photoacoustic imaging with focused laser and ultrasonic beams,” Appl. Phys. Lett. 94(3), 033902 (2009). [CrossRef]
3. T. D. Khokhlova, I. M. Pelivanov, and A. A. Karabutov, “Optoacoustic tomography utilizing focused transducers: The resolution study,” Appl. Phys. Lett. 92(2), 024105 (2008). [CrossRef]
4. 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] [PubMed]
5. J. R. Rajian, P. L. Carson, and X. Wang, “Quantitative photoacoustic measurement of tissue optical absorption spectrum aided by an optical contrast agent,” Opt. Express 17(6), 4879–4889 (2009). [CrossRef] [PubMed]
6. J. Gamelin, A. Maurudis, A. Aguirre, F. Huang, P. Guo, L. V. Wang, and Q. Zhu, “A real-time photoacoustic tomography system for small animals,” Opt. Express 17(13), 10489–10498 (2009). [CrossRef] [PubMed]
7. S. Manohar, A. Kharine, J. C. G. Hespen, W. Steenbergen, and T. G. Leeuwen, “The Twente Photoacoustic Mammoscope: system overview and performance,” Phys. Med. Biol. 50(11), 2543–2557 (2005). [CrossRef] [PubMed]
8. R. J. Talbert, S. H. Holan, and J. A. Viator, “Photoacoustic discrimination of viable and thermally coagulated blood using a two-wavelength method for burn injury monitoring,” Phys. Med. Biol. 52(7), 1815–1829 (2007). [CrossRef] [PubMed]
9. H. F. Zhang, K. Maslov, M. Sivaramakrishnan, J. Stoica, and L. V. Wang, “Imaging of hemoglobin oxygen saturation variations in single vessels in vivo using photoacoustic microscopy,” Appl. Phys. Lett. 90(5), 053901 (2007). [CrossRef]
10. 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]
11. P.-C. Li, C.-R. C. Wang, D.-B. Shieh, C. W. Wei, C. K. Liao, C. Poe, S. Jhan, A. A. Ding, and Y. N. Wu, “In vivo photoacoustic molecular imaging with simultaneous multiple selective targeting using antibody-conjugated gold nanorods,” Opt. Express 16(23), 18605–18615 (2008). [CrossRef]
12. P.-C. Li, C.-W. Wei, and Y.-L. Sheu, “Subband photoacoustic imaging for contrast improvement,” Opt. Express 16(25), 20215–20226 (2008). [CrossRef] [PubMed]
13. J. L.-S. Su, B. Wang, and S. Y. Emelianov, “Photoacoustic imaging of coronary artery stents,” Opt. Express 17(22), 19894–19901 (2009). [CrossRef] [PubMed]
14. R. J. Zemp, L. Song, R. Bitton, K. K. Shung, and L. V. Wang, “Realtime photoacoustic microscopy of murine cardiovascular dynamics,” Opt. Express 16(22), 18551–18556 (2008). [CrossRef] [PubMed]
15. L. Li, K. Maslov, G. Ku, and L. V. Wang, “Three-dimensional combined photoacoustic and optical coherence microscopy for in vivo microcirculation studies,” Opt. Express 17(19), 16450–16455 (2009). [CrossRef] [PubMed]
16. P. Ephrat, M. Roumeliotis, F. S. Prato, and J. J. L. Carson, “Four-dimensional photoacoustic imaging of moving targets,” Opt. Express 16(26), 21570–21581 (2008). [CrossRef] [PubMed]
17. Y. Xu, L. V. Wang, G. Ambartsoumian, and P. Kuchment, “Reconstructions in limited-view thermoacoustic tomography,” Med. Phys. 31(4), 724–733 (2004). [CrossRef] [PubMed]
18. G. Paltauf, J. A. Viator, S. A. Prahl, and S. L. Jacques, “Iterative reconstruction algorithm for optoacoustic imaging,” J. Acoust. Soc. Am. 112(4), 1536–1544 (2002). [CrossRef] [PubMed]
19. OAI optoacoustic Imaging Lab, “Image reconstruction as an inverse problem,” (University of Bern). http://fluor.unibe.ch/optoacoustics/reconstruction.htm
20. M. Jaeger, S. Schüpbach, A. Gertsch, M. Kitz, and M. Frenz, “Fourier reconstruction in optoacoustic imaging using truncated regularized inverse k-space interpolation,” Inverse Probl. 23(6), S51–S63 (2007). [CrossRef]
21. M. Jaeger, M. Frenz, and D. Schweizer, “Iterative Reconstruction Algorithm for Reduction of Echo Background in Optoacoustic Images,” Proc. SPIE 6856, 47 (2008).
22. M. Frenz and M. Jaeger, “Optimization of tissue irradiation in optoacoustic imaging using a linear transducer: theory and experiments,” Proc. SPIE 6856, 69 (2008).
23. R. I. Siphanto, K. K. Thumma, R. G. M. Kolkman, T. G. van Leeuwen, F. F. M. de Mul, J. W. van Neck, L. N. A. van Adrichem, and W. Steenbergen, “Serial noninvasive photoacoustic imaging of neovascularization in tumor angiogenesis,” Opt. Express 13(1), 89–95 (2005). [CrossRef] [PubMed]
24. M. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 71(1), 016706 (2005). [CrossRef] [PubMed]
25. G. J. Diebold, T. Sun, and M. I. Khan, “Photoacoustic monopole radiation in one, two, and three dimensions,” Phys. Rev. Lett. 67(24), 3384–3387 (1991). [CrossRef] [PubMed]
26. C. L. Lawson, and R. J. Hanson, Solving least squares problems (Prentice-Hall, Inc., Englewood Cliffs 1974), Chap. 23.
