Research Article  Open Access
Compressed Sensing Photoacoustic Imaging Based on Fast Alternating Direction Algorithm
Abstract
Photoacoustic imaging (PAI) has been employed to reconstruct endogenous optical contrast present in tissues. At the cost of longer calculations, a compressive sensing reconstruction scheme can achieve artifactfree imaging with fewer measurements. In this paper, an effective acceleration framework using the alternating direction method (ADM) was proposed for recovering images from limitedview and noisy observations. Results of the simulation demonstrated that the proposed algorithm could perform favorably in comparison to two recently introduced algorithms in computational efficiency and data fidelity. In particular, it ran considerably faster than these two methods. PAI with ADM can improve convergence speed with fewer ultrasonic transducers, enabling a highperformance and costeffective PAI system for biomedical applications.
1. Introduction
As an emerging biomedical imaging technique, photoacoustic imaging (PAI) has experienced considerable growth in the past decade [1]. It has been explored for molecular imaging of biomarkers [2], functional imaging of physiological parameters [3, 4] and the imaging of tumor angiogenesis [5–7] in both preclinical and clinical studies. Photoacoustic tomography (PAT) provides specklefree imaging with high contrast and high resolution which is one form of PAI. When biological tissues are irradiated by short laser pulses, some optical energy is absorbed and converted into heat. The resultant thermoplastic expansion leads to the emission of ultrasonic waves which are acquired by a singleelement focused ultrasonic transducer with mechanical scanning or an ultrasonic transducer array from a full view [8, 9]. Then the information of the tissue’s optical absorption properties can be recovered using a reconstruction algorithm.
Many algorithms have been developed to exactly or approximately reconstruct the image with a full view of data [10]. A limiting factor for the traditional filtered backprojection (FBP) algorithm is the great number of measurements made with transducers, implying long acquisition times. In addition, it is almost impossible to cover the entire surface of the tissues in many practices. To the best of our knowledge, there is no exact formula reported for limitedview PAT yet. To resolve such limiting factors, based on the compressive sensing (CS) theory [11] can be used to achieve artifactfree imaging from limitedview acquisition.
An image can be reconstructed from far fewer measurements than what the Shannon sampling theory requires if the image is sparse or can be compressed [12]. By using data from a small number of angles and an L1magic convex optimization algorithm, Provost et al. introduced the CS theory into the field of PAT [13, 14]. The issue of artifacts and loss of resolution in limitedview imaging can be addressed by using random optical illuminations for fast data acquisition via the SPGL1 algorithm [15, 16]. Sun et al. have developed an arcdirection compressedsensing PAT algorithm with numerical phantoms [17]. Both phantom and in vivo results showed that the CS method can effectively reduce undersampling artifacts via the nonlinear conjugate gradient descent algorithm [18, 19]. All of these studies have shown that CSbased reconstruction techniques can reduce the number of ultrasonic transducers of the PAT system significantly and obtain highresolution results with limitedview photoacoustic data. However, one of the critical issues that used to hinder the application of CS in PAT is the computational cost of the underlying image reconstruction process.
In this paper, we proposed a fast CS reconstruction algorithm to overcome this difficulty, leading to acceptable computational times. We studied the use of the alternating direction method (ADM) for L1norm minimization compressive sensing problems arising from sparse PAT reconstruction [20]. The proposed algorithm was used to improve the speed of the reconstruction from highly incomplete data [21, 22]. The numerical simulation results showed that the ADM algorithm was efficient and robust. In particular, the ADM can generally reduce relative errors faster than all of the other tested algorithms.
2. Method
2.1. Photoacoustic Imaging
According to the photoacoustic signal generation theory, the acoustic pressure at location and time in an acoustically homogeneous medium obeys the following wave equation [23]: where is the sound speed, is pressure, is the isobaric volume expansion coefficient, is the specific heat, and is the heating function that can be written as the product of the initial absorbed optical energy density and a temporal function of illumination . If the pulse pumping can be regarded as a Dirac delta function , the following problem can be solved using the Green’s functions to obtain the pressure:
Taking the Fourier transform on variable of (2) and denoting , the forward problem in the temporalfrequency domain is expressed as To numerically model the previously mentioned problem we used a vector to represent and a vector to represent the detected acoustic pressure . Then (3) can be expressed as , and the forward projection matrix in the temporalfrequency domain can be written as where indicates the position of the transducer, represents the sampling point in the frequency domain, and indicates the Cartesian coordinates of the image pixels.
2.2. CS Application in PAI
Mathematically, the projection matrix is ill conditioned if the measurement is insufficient. This will lead to uncertainties during the reconstruction. Fortunately, the CS theory tells us that a sparse signal can be exactly reconstructed from incomplete datasets if satisfying some requirements. It has been proven that photoacoustic image is sparse or compressed enough in a certain domain [13, 15]. By finding an appropriate sparse transform , the photoacoustic image can be reconstructed by solving a convex optimization problem in the following form [12]: When contains noise, or is not exactly sparse but only compressible, as in most practical applications, the constraint in must be relaxed, resulting in the constrained basis pursuit denoising problem where is the noise level. From the optimization theory, problem (6) is equivalent to the following problem with a suitable parameter:
As and approach zero, both problem (6) and (7) converge to problem (5).
2.3. Reconstruction Method
Based on the classical ADM technique, the firstorder primaldual algorithm that updates both primal and dual variables at each iteration was used. With an auxiliary variable problem (7) is clearly equivalent to Equation (8) has an augmented Lagrange subproblem of the form where is a Lagrange multiplier and is a penalty parameter. Given , the minimization of (9) with respect to is given by For and are fixed, the minimization of (9) with respect to is equivalent to The solution of (11) can be given explicitly by onedimensional shrinkage Finally, with a constant we updated the multiplier by In short, ADM applied to (7) produces the iteration: where both the primal and the dual variables are updated at each and every iteration.
3. Results and Discussion
To demonstrate the efficiency and superiority of the ADM algorithm in PAI reconstruction with limited angle observations, computer simulations were conducted in 2D where the imaged sources are approximately located within the transducer focal plane. All of the experiments were performed using MATLAB (MathWorks, Natick, MA, USA).
Although the existing CS algorithms provided accuracy results, the computational cost of the optimization process was significantly higher, hindering practical application. By using the ADM algorithm, in the following section, the problem with the computational cost of the image reconstruction could be overcome, leading to acceptable reconstruction computational times.
3.1. Reconstruction from Simulated LimitedView Data
In this section, based on both the Symmlet wavelet with an order of 4 and L1regularization, the numerical experiments have been conducted on a sparse 30 mm × 30 mm phantom (shown in Figure 1(a)) with a 128 × 128 resolution. In each experiment, the singleelement focused ultrasonic transducer was used to record the photoacoustic signals. To model the transducer response, the domains of and were restricted to certain values and . At every detection angle, 64 randomly chosen ’s inside the MHz window were used to completely define . By rescaling the intensity values of the phantom to , we generated measurements using in the frequency domain.
(a)
(b)
(c)
(d)
(e)
(f)
The results reconstructed by the FBP and ADM algorithms based on the multiangles observation with a different detection position are shown in Figures 1(b)1(c) and Figures 1(d)1(e), respectively. Figures 1(b) and 1(c) show image reconstruction using the FBP algorithm with measurements along a horizontal circle, stopping at the 200 and 80 positions. Figures 1(d) and 1(e) show the reconstruction results using the ADM algorithm with 80 and 40 transducers uniformly covering the 90degree view. It is shown that the results of the CS method are clearly superior to those of the FBP method. This can be seen by extracting and comparing lines from the reconstructed images in Figure 1(f).
3.2. Comparison of CS Reconstruction Algorithms
We compared ADM with L1magic and SPGL1 for the solution model (6). For all of the experiments, the same number of transducers and Fourier samples per angle uniformly covering the 90degree view was used. In order to compare these three algorithms in a way that is as independent as possible, the same iteration stopping criteria were used. The quality of the reconstructed image including the number of iterations, the CPU times, and the signaltonoise (SNR) achieved by each of the algorithms is presented in Table 1; all of which are the average values over 10 runs for each experiment.

We can conclude from Table 1 that there are large differences between the algorithm execution times: ADM can be roughly 6 times faster than SPGL1, which itself is about 20 times faster than L1magic. The proposed algorithm is not only faster, but also maintains the best SNR ratio.
In order to demonstrate the data fidelity of the ADM algorithm, Figure 2 shows the reconstruction results obtained with these three minimization schemes with measurements from 56 detection angles polluted by average SNRs of 40. In terms of reconstruction quality, all of the algorithms produced similar results overall. The denoising effect of the CS can be observed on the residual images. For a computed solution , its relative error to is defined as
(a)
(b)
(c)
3.3. Discussion
It has been demonstrated that 90 degrees is sufficient for highquality reconstruction. To study the effects of white Gaussian noise polluted measurements on the reconstruction performance for the ADM algorithm, acquisitions were simulated of different SNRs. The reconstructions are shown in Figure 3. As predicted in theory [12, 21], the proposed algorithm based on the CS scheme is robust to inaccurate measurements. By contrasting the influence of different noise levels on the reconstruction results, we can see that the increased measurement noise only increases the reconstruction noise. No major artifacts can be observed.
(a)
(b)
(c)
(d)
Figure 4 shows the tendency chart of relative errors in reconstructed photoacoustic images obtained with these three CS algorithms. It is clear that the ADM algorithm is efficient and robust. In particular, the proposed algorithm cannot only use fewer measurements to obtain better performance, but also reduce relative errors faster than other tested algorithms. However, while keeping the noise level at , the ADM algorithm sometimes reduced relative errors slower than the other two algorithms.
(a)
(b)
(c)
(d)
4. Conclusion
We presented a novel fast algorithm for PAI using a small number of angles. The proposed algorithm is based on the CS theory. Numerical simulations showed that our algorithm produced better images than FBP and other stateoftheart CS algorithms. Moreover, the proposed algorithm has been shown to be robust to noise in limitedview imaging. Ongoing work includes a more thorough experimental evaluation of ADM.
Acknowledgments
This study was supported by the National Basic Research Program of China (973 Program) under Grant no. 2011CB707700, the National Natural Science Foundation of China under Grant nos. 81227901, 61231004, 81027002, and 81261120414, the Chinese Academy of Sciences Visiting Professorship for Senior International Scientists under Grant no. 2012T1G0036, the Knowledge Innovation Project of the Chinese Academy of Sciences under Grant no. KGCX2YW907, and the Beijing Natural Science Foundation no. 4111004. X. Liu would like to thank J. Provost and Z. Yuan for beneficial discussions on the analysis of frequency signals.
References
 L. V. Wang, “Prospects of photoacoustic tomography,” Medical Physics, vol. 35, no. 12, pp. 5758–5767, 2008. View at: Publisher Site  Google Scholar
 M. Navas, A. Jimenez, and A. Asuero, “Human biomarkers in breath by photoacoustic spectroscopy,” Clinica Chimica Acta, vol. 413, no. 1516, pp. 1171–1178, 2012. View at: Publisher Site  Google Scholar
 J. Yang, C. Favazza, R. Chen et al., “Simultaneous functional photoacoustic and ultrasonic endoscopy of internal organs in vivo,” Nature Medicine, vol. 18, no. 8, pp. 1297–1302, 2012. View at: Google Scholar
 K. Homan, S. Kim, Y. S. Chen, B. Wang, S. Mallidi, and S. Emelianov, “Prospects of molecular photoacoustic imaging at 1064 nm wavelength,” Optics Letters, vol. 35, no. 15, pp. 2663–2665, 2010. View at: Publisher Site  Google Scholar
 G. Ku, X. Wang, X. Xie, G. Stoica, and L. V. Wang, “Imaging of tumor angiogenesis in rat brains in vivo by photoacoustic tomography,” Applied Optics, vol. 44, no. 5, pp. 770–775, 2005. View at: Publisher Site  Google Scholar
 R. I. Siphanto, K. K. Thumma, R. G. M. Kolkman et al., “Serial noninvasive photoacoustic imaging of neovascularization in tumor angiogenesis,” Optics Express, vol. 13, no. 1, pp. 89–95, 2005. View at: Publisher Site  Google Scholar
 A. A. Oraevsky, E. V. Savateeva, S. V. Solomatin et al., “Optoacoustic imaging of blood for visualization and diagnostics of breast cancer,” in Biomedical Optoacoustics III, Proceedings of SPIE, pp. 81–94, San Jose, Calif, USA, January 2002. View at: Publisher Site  Google Scholar
 X. Wang, Y. Pang, G. Ku, X. Xie, G. Stoica, and L. V. Wang, “Noninvasive laserinduced photoacoustic tomography for structural and functional in vivo imaging of the brain,” Nature Biotechnology, vol. 21, no. 7, pp. 803–806, 2003. View at: Publisher Site  Google Scholar
 C. Li, A. Aguirre, J. Gamelin, A. Maurudis, Q. Zhu, and L. V. Wang, “Realtime photoacoustic tomography of cortical hemodynamics in small animals,” Journal of Biomedical Optics, vol. 15, no. 1, Article ID 010509, pp. 1–3, 2010. View at: Google Scholar
 M. Xu and L. V. Wang, “Universal backprojection algorithm for photoacoustic computed tomography,” Physical Review E, vol. 71, no. 1, Article ID 016706, pp. 1–7, 2005. View at: Publisher Site  Google Scholar
 D. L. Donoho, “Compressed sensing,” IEEE Transactions on Information Theory, vol. 52, no. 4, pp. 1289–1306, 2006. View at: Publisher Site  Google Scholar
 E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Transactions on Information Theory, vol. 52, no. 2, pp. 489–509, 2006. View at: Publisher Site  Google Scholar
 J. Provost and F. Lesage, “The application of compressed sensing for photoacoustic tomography,” IEEE Transactions on Medical Imaging, vol. 28, no. 4, pp. 585–594, 2009. View at: Publisher Site  Google Scholar
 E. van den Berg and M. P. Friedlander, “Probing the pareto frontier for basis pursuit solutions,” SIAM Journal on Scientific Computing, vol. 31, no. 2, pp. 890–912, 2008. View at: Publisher Site  Google Scholar
 D. Liang, H. Zhang, and L. Ying, “Compressedsensing photoacoustic imaging based on random optical illumination,” International Journal of Functional Informatics and Personalized Medicine, vol. 2, no. 4, pp. 394–406, 2009. View at: Google Scholar
 E. Candes and J. Romberg, “L1Magic: Recovery of Sparse Signals via Convex Programming,” http://www.cms.caltech.edu/. View at: Google Scholar
 M. Sun, N. Feng, Y. Shen et al., “Photoacoustic imaging method based on arcdirection compressed sensing and multiangle observation,” Optics Express, vol. 19, no. 16, pp. 14801–14806, 2011. View at: Google Scholar
 Z. Guo, C. Li, L. Song, and L. V. Wang, “Compressed sensing in photoacoustic tomography in vivo,” Journal of Biomedical Optics, vol. 15, no. 2, Article ID 021311, pp. 1–6, 2010. View at: Google Scholar
 J. Meng, L. V. Wang, L. Ying, D. Liang, and L. Song, “Compressedsensing photoacoustic computed tomography in vivo with partially known support,” Optics Express, vol. 20, no. 15, pp. 16510–16523, 2012. View at: Google Scholar
 Y. Zhang, “YALL1: Your algorithms for L1,” http://yall1.blogs.rice.edu/. View at: Google Scholar
 J. Yang and Y. Zhang, “Alternating direction algorithms for ℓ1problems in compressive sensing,” SIAM Journal on Scientific Computing, vol. 33, no. 1, pp. 250–278, 2011. View at: Publisher Site  Google Scholar
 E. J. Candès, J. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate information,” Communications on Pure and Applied Mathematics, vol. 59, pp. 1207–1233, 2005. View at: Google Scholar
 M. Xu and L. V. Wang, “Timedomain reconstruction for thermoacoustic tomography in a spherical geometry,” IEEE Transactions on Medical Imaging, vol. 21, no. 7, pp. 814–822, 2002. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2012 Xueyan Liu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.