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 artifact-free imaging with fewer measurements. In this paper, an effective acceleration framework using the alternating direction method (ADM) was proposed for recovering images from limited-view 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 high-performance and cost-effective 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 [57] in both preclinical and clinical studies. Photoacoustic tomography (PAT) provides speckle-free 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 single-element 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 back-projection (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 limited-view PAT yet. To resolve such limiting factors, based on the compressive sensing (CS) theory [11] can be used to achieve artifact-free imaging from limited-view 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 limited-view 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 arc-direction compressed-sensing 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 CS-based reconstruction techniques can reduce the number of ultrasonic transducers of the PAT system significantly and obtain high-resolution results with limited-view 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 L1-norm 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 temporal-frequency 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 temporal-frequency 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 first-order primal-dual 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 one-dimensional 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 Limited-View Data

In this section, based on both the Symmlet wavelet with an order of 4 and L1-regularization, 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 single-element 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.

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 90-degree 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 90-degree 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 signal-to-noise (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

3.3. Discussion

It has been demonstrated that 90 degrees is sufficient for high-quality 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.

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.

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 state-of-the-art CS algorithms. Moreover, the proposed algorithm has been shown to be robust to noise in limited-view imaging. Ongoing work includes a more thorough experimental evaluation of ADM.


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. KGCX2-YW-907, 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.