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 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 [5–7] 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]: ∇²−1𝑐²𝜕²𝜕𝑡²𝑝()𝛽𝑟,𝑡=−𝐶_𝑝𝜕𝐻(),𝜕𝑡𝑟,𝑡 (1) 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: 𝑝()=𝛽𝑟,𝑡4𝜋𝐶_𝑝𝐴𝑟^′𝛿𝑡−𝑟−𝑟^′𝑐𝑑³𝑟^′𝑟−𝑟^′. (2)

Taking the Fourier transform on variable 𝑡 of (2) and denoting 𝑘=𝜔/𝑐, the forward problem in the temporal-frequency domain is expressed as 𝑝(𝑟,𝑘)=−𝑖𝑐𝑘𝛽4𝜋𝐶_𝑝𝐴𝑟^′exp𝑖𝑘𝑟−𝑟^′𝑟−𝑟^′𝑑³𝑟^′. (3) 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 Φ_{(𝑚,𝑛)(𝑖,𝑗)}=−𝑖𝑐𝑘_𝑛𝑒^{𝑖𝑘𝑛|𝑟𝑚−𝑟𝑖𝑗|}𝑟_𝑚−𝑟_𝑖𝑗, (4) 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]: min‖𝜃‖₁s.t.𝑦=ΦΨ⁻¹𝜃. (5) 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 min_𝜃‖𝜃‖₁s.t.𝑦−ΦΨ⁻¹𝜃²₂<𝜀, (6) where 𝜀 is the noise level. From the optimization theory, problem (6) is equivalent to the following problem with a suitable parameter: min_𝜃‖𝜃‖₁+12𝜇𝑦−ΦΨ⁻¹𝜃²₂. (7)

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 min_𝜃‖𝜃‖₁+1‖𝑟‖2𝜇²²s.t.𝑦=ΦΨ⁻¹𝜃+𝑟. (8) Equation (8) has an augmented Lagrange subproblem of the form min_𝜃,𝑟‖𝜃‖₁+1‖𝑟‖2𝜇²²−𝜆^𝑇ΦΨ⁻¹+𝛽𝜃+𝑟−𝑦2ΦΨ⁻¹𝜃+𝑟−𝑦², (9) where 𝜆 is a Lagrange multiplier and 𝛽>0 is a penalty parameter. Given (𝜃^𝑘,𝜆^𝑘), the minimization of (9) with respect to 𝑟 is given by 𝑟^𝑘+1=𝜇𝛽𝜆1+𝜇𝛽𝛽^𝑘−ΦΨ⁻¹𝜃−𝑦. (10) For 𝑟=𝑟^𝑘+1 and 𝜆=𝜆^𝑘 are fixed, the minimization of (9) with respect to 𝜃 is equivalent to min_𝜃‖𝜃‖₁+𝛽2ΦΨ⁻¹𝜃+𝑟^𝑘+1𝜆−𝑦−^𝑘𝛽²₂. (11) The solution of (11) can be given explicitly by one-dimensional shrinkage 𝜃^𝑘+1=Shink𝜃^𝑘−𝜏𝑔^𝑘,𝜏𝛽=max𝜃^𝑘−𝜏𝑔^𝑘−𝜏𝛽𝜃,0^𝑘−𝜏𝑔^𝑘𝜃^𝑘−𝜏𝑔^𝑘. (12) Finally, with a constant 𝛾>0 we updated the multiplier 𝜆 by 𝜆^𝑘+1=𝜆^𝑘−𝛾𝛽ΦΨ⁻¹𝜃^𝑘+1+𝑟^𝑘+1−𝑦. (13) In short, ADM applied to (7) produces the iteration: 𝑟^𝑘+1=𝜇𝛽𝜆1+𝜇𝛽𝛽^𝑘−ΦΨ⁻¹,𝜃𝜃−𝑦_𝑘+1=Shink𝜃^𝑘−𝜏𝑔^𝑘,𝜏𝛽,𝜆^𝑘+1=𝜆^𝑘−𝛾𝛽ΦΨ⁻¹𝜃^𝑘+1+𝑟^𝑘+1−𝑦, (14) 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 𝑘_𝑛/2𝜋𝑐∈[0.2,2.5] and 𝑛∈[1,128]. At every detection angle, 64 randomly chosen 𝑘_𝑛/2𝜋𝑐’s inside the [0.2,2.5] MHz window were used to completely define Φ_{(𝑚,𝑛)(𝑖,𝑗)}. By rescaling the intensity values of the phantom to [0,1], we generated measurements 𝑦 using Φ_{(𝑚,𝑛)(𝑖,𝑗)} in the frequency domain.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1 

Image reconstructions using the FBP method and ADM algorithm. (a) Original phantom. (b) and (c) image reconstructions using the FBP method with 200 and 80 transducers evenly covering the circle. (d) and (e) image reconstructions using the ADM method with 80 and 40 transducers uniformly covering the 90-degree view. (f) Center lines extracted from (a) to (e).

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 𝛿=‖𝑥^𝑘+1−𝑥^𝑘‖/‖𝑥^𝑘‖<0.005 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.

Table 1 

Numerical results for L1magic, SPGL1, and ADM methods on PAI images with different sampling angles and 64 Fourier samples per angle uniformly covering the 90-degree view.

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 (𝑥)=𝑥RelErr^′−𝑥‖𝑥‖×100%. (15)

(a)

(b)

(c)

Figure 2 

Image reconstructions in (a) with L1magic, (b) with SPGL1, and (c) with ADM, using 56 detection angles and 64 Fourier samples per angle with a uniform distribution at 90-degree curve noisy observation with SNR=40 dB.

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.

(a)

(b)

(c)

(d)

Figure 3 

Image reconstructions using the ADM algorithm with 40 detection angles that were uniformly distributed at a 90-degree curve (a) noiseless observation; noisy observation with (b) SNR=30 dB; (c) SNR=20 dB; (d) SNR=10 dB.

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 SNR=20dB, the ADM algorithm sometimes reduced relative errors slower than the other two algorithms.

(a)

(b)

(c)

(d)

Figure 4 

Relative errors for different CS reconstruction algorithms at different positions for a uniform distribution at a 90-degree curve. (a) Noiseless observation; noisy observation with (b) SNR=40 dB; (c) SNR=30 dB; (d) SNR=20 dB.

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.

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. 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.