OSA's Digital Library

Biomedical Optics Express

Biomedical Optics Express

  • Editor: Joseph A. Izatt
  • Vol. 5, Iss. 5 — May. 1, 2014
  • pp: 1363–1377
« Show journal navigation

Basis pursuit deconvolution for improving model-based reconstructed images in photoacoustic tomography

Jaya Prakash, Aditi Subramani Raju, Calvin B. Shaw, Manojit Pramanik, and Phaneendra K. Yalavarthy  »View Author Affiliations


Biomedical Optics Express, Vol. 5, Issue 5, pp. 1363-1377 (2014)
http://dx.doi.org/10.1364/BOE.5.001363


View Full Text Article

Acrobat PDF (2103 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

The model-based image reconstruction approaches in photoacoustic tomography have a distinct advantage compared to traditional analytical methods for cases where limited data is available. These methods typically deploy Tikhonov based regularization scheme to reconstruct the initial pressure from the boundary acoustic data. The model-resolution for these cases represents the blur induced by the regularization scheme. A method that utilizes this blurring model and performs the basis pursuit deconvolution to improve the quantitative accuracy of the reconstructed photoacoustic image is proposed and shown to be superior compared to other traditional methods via three numerical experiments. Moreover, this deconvolution including the building of an approximate blur matrix is achieved via the Lanczos bidagonalization (least-squares QR) making this approach attractive in real-time.

© 2014 Optical Society of America

1. Introduction

Photoacoustic (PA) imaging is an emerging, noninvasive, in vivo biomedical optical imaging modality that combines both optics and ultrasonic physics [1

1. V. E Gusev and A. A Karabutov, Laser Optoacoustics (AIP, 1993).

3

3. L. H. V. Wang and S. Hu, “Photoacoustic Tomography: In Vivo Imaging from Organelles to Organs,” Science 335, 1458–1462, (2012). [CrossRef] [PubMed]

]. Typically a nanosecond laser pulse irradiates biological tissue at near-infrared (NIR) wavelength for deep penetration. The rise in temperature (in the order of milli Kelvin) due to the light absorption by tissue causes emission of pressure waves via thermoelastic expansion. The initial pressure rise is proportional to the absorbed optical energy and the Grueneisen parameter (a dimensionless parameter of the tissue). This pressure wave travels in the soft biological tissues as an acoustic wave, also known as PA wave. A wideband ultrasonic transducer (UST) acquires the propagated PA waves outside the tissue boundary. The UST can be either an array of ultrasonic detectors [4

4. H. F. Zhang, K. Maslov, G. Stoica, and L. H. V. Wang, “Functional photoacoustic microscopy for high-resolution and noninvasive in vivo imaging,” Nat. Biotech. 24, 848–851 (2006). [CrossRef]

, 5

5. B. Yin, D. Xing, Y. Wang, Y. Zeng, Y. Tan, and Q. Chen, “Fast photoacoustic imaging system based on 320-element linear transducer array,” Phys. Med. Bio. 49, 1339–1346 (2004). [CrossRef]

] or a single-element ultrasonic transducer [6

6. G. Ku and L. H. V. Wang, “Scanning thermoacoustic tomography in biological tissue,” Med. Phys. 27, 1195–1202 (2000). [CrossRef] [PubMed]

, 7

7. M. Pramanik, G. Ku, C. H. Li, and L. H. V. Wang, “Design and evaluation of a novel breast cancer detection system combining both thermoacoustic (TA) and photoacoustic (PA) tomography,” Med. Phys. 35, 2218–2223 (2008). [CrossRef] [PubMed]

]. The PA wave that is collected outside the tissue boundary is used to map the initial pressure rise (or the absorbed optical energy density) within the tissue with the help of a reconstruction algorithm. Structural and functional PA imaging for both pre-clinical and clinical applications has been demonstrated in the literature [4

4. H. F. Zhang, K. Maslov, G. Stoica, and L. H. V. Wang, “Functional photoacoustic microscopy for high-resolution and noninvasive in vivo imaging,” Nat. Biotech. 24, 848–851 (2006). [CrossRef]

, 5

5. B. Yin, D. Xing, Y. Wang, Y. Zeng, Y. Tan, and Q. Chen, “Fast photoacoustic imaging system based on 320-element linear transducer array,” Phys. Med. Bio. 49, 1339–1346 (2004). [CrossRef]

, 8

8. M. Pramanik and L. H. V. Wang, “Thermoacoustic and photoacoustic sensing of temperature,” J. Biomed. Opt. 14, 054024 (2009). [CrossRef] [PubMed]

, 9

9. K. H. Song and L. H. V. Wang, “Noninvasive photoacoustic imaging of the thoracic cavity and the kidney in small and large animals,” Med. Phys. 35, 4524–4529 (2008). [CrossRef] [PubMed]

]. In addition, with the help of targeted contrast agents, the photoacoustic imaging has been shown to be a strong contender for the molecular imaging [10

10. Y. Wang, X. Y. Xie, X. D. Wang, G. Ku, K. L. Gill, D. P. O’Neal, G. Stoica, and L. H. V. Wang, “Photoacoustic Tomography of a Nanoshell Contrast Agent in the in Vivo Rat Brain,” Nano Lett. 4, 1689–1692 (2004). [CrossRef]

, 11

11. A. De la Zerda, C. Zavaleta, S. Keren, S. Vaithilingam, S. Bodapati, Z. Liu, J. Levi, B. R. Smith, T. J. Ma, O. Oralkan, Z. Cheng, X. Y. Chen, H. J. Dai, B. T. Khuri-Yakub, and S. S. Gambhir, “Carbon nanotubes as photoacoustic molecular imaging agents in living mice,” Nat. Nanotech. 3, 557–562 (2008). [CrossRef]

].

The critical step in PA imaging is the image reconstruction, where the recent emphasis has been on obtaining quantitatively accurate PA images [12

12. P. Kuchment and L. Kunyansky, “Mathematics of thermoacoustic and photoacoustic tomography,” European J. App. Math. 19, 191–224 (2008).

18

18. C. Huang, K. Wang, L. Nie, L. V. Wang, and M. A. Anastasio, “Full-Wave Iterative Image Reconstruction in Photoacoustic Tomography with Acoustically Inhomogeneous Media,” IEEE Trans. Med. Imag. 32, 1097–1110 (2013). [CrossRef]

]. Several photoacoustic image reconstruction algorithms exist in the literature, including analytic algorithms like filtered back projection (FBP), time-reversal, and Fourier transform based reconstruction [12

12. P. Kuchment and L. Kunyansky, “Mathematics of thermoacoustic and photoacoustic tomography,” European J. App. Math. 19, 191–224 (2008).

17

17. K. Wang, R. Su, A. A. Oraevsky, and M. A. Anastasio, “Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography,” Phys. Med. Bio. 57, 5399–5423 (2012). [CrossRef]

]. These algorithms are based on the spherical Radon transform model, which has an inherent limitation of requiring large number of data points around the target object for accurately estimating the initial pressure distribution. Large number of data points in turn require expensive transducer arrays or long data acquisition time (if single element transducers is used). Moreover, analytic reconstruction does not provide the desired quantitative accuracy with less number of data points. To overcome these limitations, several iterative image reconstruction algorithms [15

15. P. Ephrat, L. Keenliside, A. Seabrook, F. S. Prato, and J. J. L. Carson, “Three-dimensional photoacoustic imaging by sparse-array detection and iterative image reconstruction,” J. Biomed. Opt. 13, 054052 (2008). [CrossRef] [PubMed]

18

18. C. Huang, K. Wang, L. Nie, L. V. Wang, and M. A. Anastasio, “Full-Wave Iterative Image Reconstruction in Photoacoustic Tomography with Acoustically Inhomogeneous Media,” IEEE Trans. Med. Imag. 32, 1097–1110 (2013). [CrossRef]

] have been proposed to improve the reconstructed image quality, whose requirement for number of measurements are comparatively less.

As with any model-based image reconstruction methods, deployment of regularization parameter will inherently blur the reconstructed images [24

24. J. A. Fessler and W. L. Rogers, “Spatial resolution properties of penalized- likelihood image reconstruction methods: Space-invariant tomographs,” IEEE Trans. Image Process. 5, 1346–1358 (1996). [CrossRef]

, 25

25. J. Prakash, H. Dehghani, B. W. Pogue, and P. K. Yalavarthy, “Model-Resolution based Basis Pursuit Deconvolution Improves Diffuse Optical Tomographic Imaging,” IEEE Trans. Med. Imag. 2014 (in press). [CrossRef]

]. These reconstructed images can be deblurred for better recovery of the internal tissue structure as long as the source of the blurring is known or can be modeled. Recently, blind deconvolution method has shown to improve the resolution in optical-resolution photoacoustic microscopy (OR-PAM) [26

26. J. Chen, R. Lin, H. Wang, J. Meng, H. Zheng, and L. Song, “Blind-deconvolution optical-resolution photoacoustic microscopy in vivo,” Opt. Express 21, 7316–7327 (2013). [CrossRef] [PubMed]

], in which the point spread function (PSF) was estimated. If the PSF is known then a constrained deconvolution method or least square filtering methods could be used for performing the image deconvolution [27

27. B. R. Hunt, “The application of constrained least square estimation to image restoration by digital computer,” IEEE Trans. Comput. C-22, 805–812 (1973). [CrossRef]

]. Compressive sensing techniques like basis pursuit deconvolution (BPD) has been used in the field of image analysis, and found that this mechanism has the ability to retain the structural information [28

28. E. Candes and M. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag. 25, 21–30 (2008). [CrossRef]

, 29

29. J. Romberg, “Imaging via compressive sampling,” IEEE Signal Process. Mag. 25, 14–20 (2008). [CrossRef]

]. The BPD algorithms deploy l1-norm based minimization of the generalized least-square objective function. Basis pursuit denoising has been used previously in photoacoustic imaging to perform efficient reconstruction using fewer measurements [30

30. J. Provost and F. Lesage, “The application of compressed sensing for photo-acoustic tomography,” IEEE Trans. Med. Imag. 28, 585–594 (2009). [CrossRef]

, 31

31. Z. Guo, C. H. Li, L. Song, and L. H. V. Wang, “Compressed sensing in photoacoustic tomography in vivo,” J. Biomed. Opt. 15, 021311 (2010). [CrossRef] [PubMed]

].

The main aim of this work is to remove the blur introduced by regularization parameter, where the blurring model can be built by using the model-resolution characteristics, making the image reconstruction problem a two-step approach, with initial step being model-based traditional reconstruction step and second one being the deblurring step. The deblur-ring/deconvolution of the reconstructed image is achieved via BPD. More importantly, the computation of blur matrix is achieved through Least squares QR (LSQR) approach, which provides an computationally efficient dimensionality reduction procedure. Split augmented Lagrangian shrinkage algorithm (SALSA) is a well known algorithm for performing BPD and the same has been used here [32

32. M. Figueiredo, J. Bioucas-Dias, and M. Afonso, “Fast Frame-Based Image Deconvolution Using Variable Splitting and Constrained Optimization,” IEEE Worskhop on Statistical Signal Processing, Cardiff, Wales, 2009.

]. The original SALSA algorithm was based on sparsity signal processing toolbox [33

33. I. Selesnick, “Introduction to Sparsity in Signal Processing [Connexions Web site],” Available at: http://cnx.org/content/m43545/1.3/. (2012).

]. In this work, it has been shown that model-resolution based deblurring always results in qualitatively/quantitatively more accurate reconstruction compared to traditional model-based reconstruction, provided the regularization parameter was choosen within a reasonable bounds. The proposed method has the distinct advantage of removing the unwarranted bias in the reconstructed PA image introduced by the heuristic choice of regularization parameter.

2. Analytical and LSQR based reconstruction for photoacoustic tomography

2.1. k-wave based time reversal (analytical)

The k-Wave-based time-reversal is a one-step image reconstruction method, which is computationally efficient, for a planar measurement surface. The mathematical model governing photoacoustic tomography is [34

34. Y. Hristova, P. Kuchment, and L. V. Nguyen, “Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media,” Inv. Problems 24, 055006 (2008). [CrossRef]

]
p2t2=c2(r)Δrp,t0,rR2p(r,0)=f(r),pt(r,0)=0,p(y,t)=g(y,t),foryB,t0
(1)
Here g(y, t) is the time series pressure data measured at transducer’s location yB at time t. The pt and p2t2 are the first and second time derivatives of pressure (p) and c(r) represents the speed of the ultrasound propagation (kept constant in here). The Δr is the Laplacian function with respect to space r. The aim here is to estimate the initial value of f (r) = p(r, t) at t=0, given the measured ultrasound data g(y, t) and c(r).

There are many approaches to solve this problem, namely Filtered Back-projection, Fourier Reconstruction and Time reversal methods (all of these are analytical in nature) [17

17. K. Wang, R. Su, A. A. Oraevsky, and M. A. Anastasio, “Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography,” Phys. Med. Bio. 57, 5399–5423 (2012). [CrossRef]

]. The k-wave tool box, which is an open-source based on Matlab, provides analytical tools for performing the initial pressure reconstruction [35

35. B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt. 15, 021314 (2010). [CrossRef] [PubMed]

]. It is important to note that k-wave can perform full-wave and time-reversal based image reconstruction [18

18. C. Huang, K. Wang, L. Nie, L. V. Wang, and M. A. Anastasio, “Full-Wave Iterative Image Reconstruction in Photoacoustic Tomography with Acoustically Inhomogeneous Media,” IEEE Trans. Med. Imag. 32, 1097–1110 (2013). [CrossRef]

, 35

35. B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt. 15, 021314 (2010). [CrossRef] [PubMed]

], and in this work time-reversal based reconstruction was used for comparison. The time-reversal method states that for a time point L (the longest time for the PA wave to traverse the domain ω), the solution p(r, t) vanishes inside ω for any t > L [34

34. Y. Hristova, P. Kuchment, and L. V. Nguyen, “Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media,” Inv. Problems 24, 055006 (2008). [CrossRef]

]. A zero initial condition can now be imposed at t = L and boundary conditions equal to the measured data g and solve the above mathematical model in the time reversal direction, thus arriving at t = 0 to the function f (r) = p(r, 0) [34

34. Y. Hristova, P. Kuchment, and L. V. Nguyen, “Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media,” Inv. Problems 24, 055006 (2008). [CrossRef]

]. In this work, a interpolated measurement vector is used to estimate the initial pressure of the imaging domain under consideration [35

35. B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt. 15, 021314 (2010). [CrossRef] [PubMed]

].

2.2. System matrix approach

The process of collection of PA waves outside the boundary can be represented as a time-varying causal system [17

17. K. Wang, R. Su, A. A. Oraevsky, and M. A. Anastasio, “Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography,” Phys. Med. Bio. 57, 5399–5423 (2012). [CrossRef]

]. An image having a dimension of n × n is vectorized by stacking the columns into n2 × 1 vector, represented as x. Each column of the system matrix (A, having dimension of m × n2) is the system’s response to a corresponding entry in the vectorized image (x). It is important to note that the columns of the data (time-varying) are also stacked to result in a long vector having a dimension of m × 1. The system matrix is built as a block circulant matrix, as explained in Ref. [22

22. C. B. Shaw, J. Prakash, M. Pramanik, and P. K. Yalavarthy, “Least Squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography,” J. Biomed. Opt. 18, 080501:1–3 (2013).

], with limited transducer bandwidth.

The system response for a corner pixel is recorded using k-Wave MATLAB toolbox [35

35. B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt. 15, 021314 (2010). [CrossRef] [PubMed]

], similar to Ref. [22

22. C. B. Shaw, J. Prakash, M. Pramanik, and P. K. Yalavarthy, “Least Squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography,” J. Biomed. Opt. 18, 080501:1–3 (2013).

]. This system response was found using the k-space pseudo spectral solution method with two coupled first order equations [35

35. B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt. 15, 021314 (2010). [CrossRef] [PubMed]

]. The simulation geometry of data-collection is shown in Fig. 1. A computational grid indicated by a black square in Fig. 1, having a size of 501 × 501 pixels (0.1 mm/pixel) was used and the detectors (small circles) were placed on a circle of 22 mm radius. Sixty detectors were placed equidistantly around this full circle. The detectors were considered to be a point detector having a center frequency of 2.25 MHz and 70% bandwidth, which were given as set of input parameters into the k-wave simulation for the detection mechanism. Although in practice, large area detectors are used, for simplicity, we have assumed the detectors to be point detector. However, the center frequency and the bandwidth of the point detectors are kept same as the large transducers used in practical PAT systems. A perfectly matched layer (PML) was used to satisfy the boundary conditions. The imaging region (filled green square) was restricted to 201 × 201 pixels located at the center. For data collection, the time step was chosen to be 50 ns with recorded number of time steps being 500. For simplicity, the simulations assumed the speed of sound to be 1500 m/s and the medium to be homogeneous with no absorption or dispersion of sound.

Fig. 1 Schematic of Photoacoustic data acquisition set-up with shaded square region indicating the imaging domain.

The forward model of PA imaging can be summarized as
Ax=b
(2)
where x is a long column vector having a dimension of n2 × 1 (n = 201) and b is a measurement vector with a dimension of m × 1 (m = 60 × 500). Back projection (analytical) type image reconstruction scheme will be [36

36. G. Paltauf, J. A. Viator, S. A. Prahl, and S. L. Jacques, “Iterative reconstruction algorithm for optoacoustic imaging,” J. Acous. Soc. Am. 112, 1536–1544 (2002). [CrossRef]

]
xbp=ATb
(3)
where, T indicates the transpose of the matrix. Backprojection methods are non-iterative in nature making it computationally efficient, but known to provide only qualitative results, and hence may not be optimal for quantitative imaging [36

36. G. Paltauf, J. A. Viator, S. A. Prahl, and S. L. Jacques, “Iterative reconstruction algorithm for optoacoustic imaging,” J. Acous. Soc. Am. 112, 1536–1544 (2002). [CrossRef]

].

2.3. Model based method - Tikhonov regularization (ℓ2-norm based)

The model-based reconstruction relies on minimizing the residue of the linear equations along with a regularization functional having a smoothness constraint, known as Tikhonov Regularization, where the objective function can be written as
Ω=Axb22+λx22
(4)
where, λ is the regularization parameter, providing the balance between residue of the linear equations (first term on the right-hand side) and expected initial pressure distribution (x). The 2-norm is represented by .22. The function Ω is minimized with respect to x, resulting in,
xTikh=(ATA+λI)1ATb
(5)
The regularization parameter dictates the reconstructed initial pressure distribution characteristics. Higher regularization tends to over-smooth the image while lower λ values amplifies the noise in the images.

2.4. LSQR-based reconstruction method

Using Eqs. 6, 7, and 8 in the argument of Eq. 4, results in [23

23. J. Prakash and P. K. Yalavarthy, “A LSQR-type method provides a computationally efficient automated optimal choice of regularization parameter in diffuse optical tomography,” Med. Phys. 40, 033101 (2013). [CrossRef] [PubMed]

]
bAx=Uk+1(β0e1Bkx(k));x=Vkx(k),
(10)
with x(k) being the dimensionality reduced version of x with k < n2 [38

38. M. E. Kilmer and D. P. OLeary, “Choosing regularization parameters in iterative methods for ill-posed problems,” SIAM J. Matrix Anal. Appl. 22, 1204–1221 (2001). [CrossRef]

]. Substitution of Eq. 10 in Eq. 4 and using the property Uk+1T Uk+1 = I results in,
Ω˜=β0e1Bkx(k)22+λx(k)22
(11)
Considering the first order condition of the Eq. 11, the new update equation becomes [22

22. C. B. Shaw, J. Prakash, M. Pramanik, and P. K. Yalavarthy, “Least Squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography,” J. Biomed. Opt. 18, 080501:1–3 (2013).

, 23

23. J. Prakash and P. K. Yalavarthy, “A LSQR-type method provides a computationally efficient automated optimal choice of regularization parameter in diffuse optical tomography,” Med. Phys. 40, 033101 (2013). [CrossRef] [PubMed]

]
xest(k)=(BkTBk+λI)1β0BkTe1;xLSQR=Vkxest(k)
(12)
where xest(k) is the estimated version of x(k) for a given regularization parameter λ and fixed k as explained in Ref. [38

38. M. E. Kilmer and D. P. OLeary, “Choosing regularization parameters in iterative methods for ill-posed problems,” SIAM J. Matrix Anal. Appl. 22, 1204–1221 (2001). [CrossRef]

]. For the case, k = n2, xLSQR = xTikh for a given value of λ. For all other cases, xLSQR becomes an approximation to xTikh. It is well-known that there exists an optimal k with kn2, for which solutions xLSQR and xTikh are close within the precision limits of the computation [23

23. J. Prakash and P. K. Yalavarthy, “A LSQR-type method provides a computationally efficient automated optimal choice of regularization parameter in diffuse optical tomography,” Med. Phys. 40, 033101 (2013). [CrossRef] [PubMed]

, 37

37. C. C. Paige and M. A. Saunders, “LSQR: An algorithm for sparse linear equations and sparse least squares,” ACM Trans. Math. Software 8, 43–71 (1982). [CrossRef]

].

2.5. Estimation of optimal λ using a LSQR-type method

The LSQR type algorithm can be used for calculating the xest(k) (Eq. 12) and then multiplying it by the right Lanczos matrix (Vk) to obtain the approximate solution. This kind of evaluation of update turns out to be computationally more efficient, compared to traditional way of finding update using Eq. 5, as dimensionality reduction is achieved using Lanczos algorithm. Hence, a method of estimating the optimal regularization parameter (λ) was previously proposed in Ref. [22

22. C. B. Shaw, J. Prakash, M. Pramanik, and P. K. Yalavarthy, “Least Squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography,” J. Biomed. Opt. 18, 080501:1–3 (2013).

, 23

23. J. Prakash and P. K. Yalavarthy, “A LSQR-type method provides a computationally efficient automated optimal choice of regularization parameter in diffuse optical tomography,” Med. Phys. 40, 033101 (2013). [CrossRef] [PubMed]

] by minimizing the residue of the linear equations i.e., ||AxLSQRb||2. The estimation of the number of Lanczos iterations was also performed automatically in Ref. [22

22. C. B. Shaw, J. Prakash, M. Pramanik, and P. K. Yalavarthy, “Least Squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography,” J. Biomed. Opt. 18, 080501:1–3 (2013).

, 23

23. J. Prakash and P. K. Yalavarthy, “A LSQR-type method provides a computationally efficient automated optimal choice of regularization parameter in diffuse optical tomography,” Med. Phys. 40, 033101 (2013). [CrossRef] [PubMed]

]. In this work, the optimal number of Lanczos iterations was found to be 25. The Lanczos bidiagonalization was performed using the Matlab based regularization tools [39

39. P. C. Hansen, “Regularization tools version 4.0 for Matlab 7.3,” Numerical Algorithms 46, 189–194 (2007). [CrossRef]

].

3. Model-resolution based deconvolution

3.1. LSQR based model resolution matrix

The model resolution matrix can be estimated using the bidiagonal matrix. Substituting the Eq. 6 in Eq. 12 results in,
xest(k)=(BkTBk+λI)1BkTUk+1Tb
(13)
Substituting Eq. 2 in the above equation gives,
xest(k)=(BkTBk+λI)1BkTUk+1TAx
(14)
Using Eq. 10 converts Eq. 14 into,
xest(k)=(BkTBk+λI)1BkTUk+1TAVkx(k)
(15)
Rewriting Eq. 7 as,
Uk+1TAVk=Bk
(16)
Using the above (Eq. 16) in Eq. 15 results in,
xest(k)=(BkTBk+λI)1BkTBkx(k)
(17)
which relates the expected x(k) with its estimated version ( xest(k)). If the regularization parameter (λ) is equal to zero, then only xest(k)=x(k). The regularization parameter is always greater than 0, indicating that xest(k)x(k). Rewriting Eq. 17 as,
xest(k)=Mx(k)
(18)
with
M=(BkTBk+λI)1BkTBk
(19)
representing the model resolution matrix (also known as blur matrix) for x(k), having dimension of k × k. Note that the aim of our deconvolution/deblurring is to get an estimate of x(k) using M and xest(k). If the regularization parameter is equal to 0, then estimating the blur matrix (M) will not be possible (BkT Bk is a rank-deficient highly ill-conditioned matrix).

3.2. Basis pursuit deconvolution in photoacoustic tomography

Fig. 2 Flowchart showing the major steps used in the proposed method.

Algorithm 1:. Spilt augmented lagrangian shrinkage algorithm (SALSA) [41]

table-icon
View This Table
| View All Tables

4. Figures of merit

4.1. Pearson Correlation (PC)

The performance of the proposed method in terms of improving the reconstructed image quality was evaluated using a quantitative metric, namely Pearson Correlation (PC). The Pearson Correlation is defined as [42

42. J. Prakash, C. B. Shaw, R. Manjappa, R. Kanhirodan, and P. K. Yalavarthy, “Sparse Recovery Methods Hold Promise for Diffuse Optical Tomographic Image Reconstruction,” IEEE J. Sel. Top. Quantum Electron. 20, 6800609 (2014). [CrossRef]

]
PC(x,xrecon)=COV(x,xrecon)σ(x)σ(xrecon)
(22)
where x is the expected initial pressure distribution and xrecon is the reconstructed initial pressure distribution using one of the above explained methods. Here σ indicates the standard deviation and COV is the covariance. This measures the degree of correlation between the target (expected) and the reconstructed image, having values ranging from −1 to 1. The higher the value of PC, the better is the detectability of the targets in comparison to the expected image [42

42. J. Prakash, C. B. Shaw, R. Manjappa, R. Kanhirodan, and P. K. Yalavarthy, “Sparse Recovery Methods Hold Promise for Diffuse Optical Tomographic Image Reconstruction,” IEEE J. Sel. Top. Quantum Electron. 20, 6800609 (2014). [CrossRef]

].

4.2. Contrast to Noise Ratio (CNR)

The other figure of merit that is used is the contrast to noise ratio (CNR), which is defined as,
CNR=μroiμback(σroi2aroi+σback2aback)12
(23)
where μ and σ terms represents the mean and the variance corresponding to the region of interest (ROI, in here the objects) and the background in the reconstructed initial pressure (x). The aroi=AroiAtot and aback=AbackAtot represents the area ratio. Higher the value of the CNR better is the reconstruction performance.

5. Numerical experiments

It is important to note that measuring actual initial pressure rise in experimental phantom is challenging, leading to difficulties in comparing the quantitative accuracies of the reconstruction methods, making numerical experiments ideal for such comparisons. In order to show the effectiveness of the proposed method, a variation of the derenzo phantom was chosen. This phantom is shown in Fig. 3(a), consists of small and large targets distributions over the imaging region. The targets are grouped as 1–6 according to the sizes as shown in Fig. 3(a). The maximum initial pressure was kept at 1 kPa for the objects.

Fig. 3 (a) Derenzo phantom that is used for the study of resolution characteristics (dimensions are in millimeters). (b–h) Reconstructed photoacoustic images using k-wave interpolated, LSQR with optimal choice of λ, LSQR with heuristic choice of λ, Basis Pursuit Deconvolution (BPD) with optimal choice of λ in LSQR framework, BPD with heuristic choice of λ and 40 dB noise in LSQR framework, BPD with heuristic choice of λ and 30 dB noise in LSQR framework, BPD with heuristic choice of λ and 20 dB noise in LSQR framework, respectively.

Due to the high intrinsic optical contrast (e.g. blood), PA imaging is widely used for visualizing internal blood vessel structures, both in the brain as well as other areas. Another numerical blood vessel phantom as shown in Fig. 4(a), with initial pressure rise as 1 kPa was also used to demonstrate the performance of the algorithm. In terms of comparing the performance of the proposed method with others in terms of reconstruction of sharp edges, a numerical experiment consisting of targets as alphabets ‘PAT’ (shown in Fig. 5(a)) was performed. The k-wave toolbox was used to generate the PA signal with 60 detectors in both the experiments. The collected data was then added with 1% noise in all cases, resulting in signal-to-noise-ratio (SNR) of 40 dB. To test the effectiveness of the proposed scheme for increasing noise levels, the data collection for the derenzo phantom was tested with 20 and 30 dB SNR levels. This data was used for performing the reconstruction using traditional LSQR based reconstruction (without deconvolution) and BPD based two-step approaches (refer to Fig. 2). The k-wave based time reversal reconstruction is also performed for comparison, here this was based on interpolation of the detectors as explained in Ref. [35

35. B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt. 15, 021314 (2010). [CrossRef] [PubMed]

]. A Linux workstation with dual six-core Intel Xeon processor having a speed of 2.66 GHz with 64 GB RAM has been used for all computations performed in this work.

Fig. 4 (a) Numerical blood vessel phantom that is used for the study (dimensions are in millimeters). (b–f) Reconstructed photoacoustic images using k-wave interpolated, LSQR with heuristic choice of λ, LSQR with optimal choice of λ, Basis pursuit deconvolution (BPD) with heuristic choice of λ in LSQR framework, BPD with optimal choice of λ in LSQR framework, respectively. (g) One-dimensional cross-sectional plot for the results presented in (a),(c),(d),(e), and (f) along the dotted line shown in (a).
Fig. 5 (a) Numerical PAT phantom that is used for the study (dimensions are in millimeters). (b–f) Reconstructed photoacoustic images using k-wave interpolated, LSQR with heuristic choice of λ, LSQR with optimal choice of λ, Basis pursuit deconvolution (BPD) with heuristic choice of λ in LSQR framework, BPD with optimal choice of λ in LSQR framework, respectively. (g) One-dimensional cross-sectional plot for the results presented in (a),(c),(d),(e), and (f) along the dotted line shown in (a).

6. Results and discussion

The reconstructed initial pressure distribution using k-wave time-reversal reconstruction algorithm for derenzo phantom and blood vessel network is shown in Fig. 3(b) and Fig 4(b), respectively. The reconstruction was also performed using the model based reconstruction technique by using the regularization parameter as 0.01 and optimally calculating the regularization parameter using LSQR approach [22

22. C. B. Shaw, J. Prakash, M. Pramanik, and P. K. Yalavarthy, “Least Squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography,” J. Biomed. Opt. 18, 080501:1–3 (2013).

]. The reconstruction using the LSQR methods with heuristic choice of regularization for the derenzo and blood vessel phantom is shown in Fig. 3(d) and 4(c), respectively. Optimal choice of regularization parameter was estimated as explained in Ref. [22

22. C. B. Shaw, J. Prakash, M. Pramanik, and P. K. Yalavarthy, “Least Squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography,” J. Biomed. Opt. 18, 080501:1–3 (2013).

] and the pressure distribution for this case with derenzo and blood vessel phantom is indicated in Fig. 3(c) and 4(d), respectively. The optimal λ was found to be 0.004 and 0.0036 for the derenzo and the blood vessel phantom, respectively.

The LSQR based approximate model resolution matrix was computed with heuristic choice and optimal choice of regularization parameter. This model resolution matrix was then used to perform the deblurring operation. The initial pressure distribution after performing the deconvolution operation (heuristic choice of λ = 0.01) with the derenzo and the blood vessel phantom is shown in the Fig. 3(f) and 4(e), respectively. The pressure distribution with optimal choice of regularization and deconvolution procedure for the derenzo and the blood vessel phantom is shown in Fig. 3(e) and 4(f), respectively. The initial pressure distribution using the proposed method for varying data SNR levels i.e. 30 dB and 20 dB using the derenzo phantom is shown in Fig. 3(g) and 3(h) respectively. The quantitative metrics (figures of merit) of the reconstructions obtained using other techniques is compiled in Table-1 for comparison (actual reconstructed images are not shown). The results show that the proposed method provides the required noise-tolerance, making it deployable in an experimental scenario, and is able to obtain reasonably more accurate results compared to other methods, when the noise level is high (SNR = 20 dB).

Table 1. The Pearson Correlation (PC) and Contrast to Noise Ratio (CNR) of the reconstructed initial pressure distribution for the results presented in this work.

table-icon
View This Table
| View All Tables

From Fig. 3(b) it is clearly evident that for limited data points (60 transducer positions, whereas typical PAT system uses 120–250 or more transducer positions [5

5. B. Yin, D. Xing, Y. Wang, Y. Zeng, Y. Tan, and Q. Chen, “Fast photoacoustic imaging system based on 320-element linear transducer array,” Phys. Med. Bio. 49, 1339–1346 (2004). [CrossRef]

8

8. M. Pramanik and L. H. V. Wang, “Thermoacoustic and photoacoustic sensing of temperature,” J. Biomed. Opt. 14, 054024 (2009). [CrossRef] [PubMed]

]) k-wave based time reversal reconstruction is inadequate in recovering the original object, even with interpolated data. Also the resolution of the reconstructed objects are position dependent. For group 1 objects the resolution near the scanning centre is much better than those located near the boundary as shown in Fig. 3(b), this degradation in resolution is due to the low numbers of detector position used for data collection as well as the interpolation used. However, model based reconstruction was able to produce location independent object resolution as seen in Fig. 3(c–h).

Another challenge in photoacoustic image reconstruction is the object shape recovery. For object groups 3–6, which are large in size, the reconstructed object looks like donut/ring [Fig. 3(b–d)]. This can be attributed to the limited bandwidth of the detection system. As described earlier, the transducer used for the study have 2.25 MHz center frequency with 70% bandwidth. However with BPD based reconstruction, the object shape is recovered [Fig. 3(e–h)] very close to the original shape. It can be seen that the proposed BPD based reconstruction was able to provide better quantitative accuracy compared to analytical (k-wave time-reversal) based and model based methods (LSQR with heuristic λ and LSQR with optimal λ). The initial pressure rise values show the quantitative accuracy of the BPD based reconstruction.

The line profile for the target and the reconstructed pressure distribution using blood vessel phantom is shown in Fig. 4(g). The line profile was taken along the dotted line as shown in Fig. 4(a). The line profile for the k-wave based time reversal reconstructed image is not shown as visually it is evident that k-wave based time reversal reconstructed image was very poor in quality without delineating the blood vessel structures. Arrows in Fig. 4(b) indicates the blood vessels which were not reconstructed by k-wave. This is due to the limited data points (60) used in k-wave reconstruction. From the line profile, it is clearly evident that BPD based reconstruction, was able to recover the initial pressure rise upto 0.8 kPa compared to the other methods which gave a quantitative value of only 0.3 kPa. The reconstruction results pertaining to ‘PAT’ phantom are as shown in Fig. 5(b–f). It can be seen that the reconstruction using the proposed method yields better contrast compared to the established methods. The optimal choice of regularization parameter (λ) in this case was found to be 1.205e-4. A line profile is shown in Fig. 5(g) that was taken along the dotted line in Fig. 5(a).

Figure 6 shows a plot of residual error ( ω=AxLSQRb22) versus the number of Lanczos iterations for the results pertaining to Fig. 4 in terms of choosing optimal k, in here k = 25. Similar behavior was observed for other cases discussed in this work. Note that the relative difference in the residual error (ωkωk−1) beyond k = 25 was less than 10−6 (single precision limits). The optimal k (25 for all cases discussed here) is chosen to be the one, when the change in residual error is in single precision limits. Note that the typical computational time taken for each of the method that was discussed in this work is reported in Table-2. This indicates that performing the additional step of BPD has negligible computational burden (2 seconds).

Fig. 6 Plot of the variation of the Residual error and the number of Lanczos iterations for numerical blood vessel phantom (Fig. 4).

Table 2. Typical computational time (in seconds) for reconstructing the initial pressure distribution for the various methods presented in this work. Note that the time taken for building the system matrix (375 seconds) was excluded from this.

table-icon
View This Table
| View All Tables

The aim of this work is to introduce a basis pursuit deconvolution approach for improving the reconstructed pressure distribution in PAT, where an approximate blur matrix is built using the efficient LSQR approach. The line profile in Fig. 4 indicates that the BPD method results in similar reconstruction with heuristic and optimal choice of regularization parameter (Figs. 4(e) and (f)). The reconstructed pressure distribution in Figs. 3 and 4 indicate that the inclusion of additional step of deconvolution provides better reconstruction than just performing the image reconstruction using the LSQR type approach. The deconvolution process does not require the optimal choice of regularization parameter (λ) as long as it is chosen within the reasonable bounds. These bounds are the values of λ for which
x(k)xd(k)2<x(k)xest(k)2
(24)
where xest(k) is a function of λ. In simple words, unless xest(k) is close to the expected x(k), this relation will not hold good, which essentially means that the BPD methods works well for all cases where xest(k) is a meaningful result.

There were earlier attempts of performing deconvolution of photoacoustic images [26

26. J. Chen, R. Lin, H. Wang, J. Meng, H. Zheng, and L. Song, “Blind-deconvolution optical-resolution photoacoustic microscopy in vivo,” Opt. Express 21, 7316–7327 (2013). [CrossRef] [PubMed]

] using blind deconvolution techniques assuming the blur kernel was not known. This work introduced a framework to build the blur matrix based on model-resolution and applied basis pursuit deconvolution to improve the reconstructed pressure distribution. The building of the model resolution matrix was done using a LSQR approach (using bidiagonal matrix) hence providing the computational efficiency than building the model resolution matrix using the original system matrix as explained in Ref. [25

25. J. Prakash, H. Dehghani, B. W. Pogue, and P. K. Yalavarthy, “Model-Resolution based Basis Pursuit Deconvolution Improves Diffuse Optical Tomographic Imaging,” IEEE Trans. Med. Imag. 2014 (in press). [CrossRef]

].

It is also important to note that the proposed method was able to provide better quantitation compared all other techniques with limited number of detectors (60). Thus, deblurring holds a promise in providing better reconstructed images with less number of detectors, additionally reducing the acquisition time and system cost. The blur that arises due to the laser excitation pulse width or the limited detector bandwidth can also be corrected in the signal domain using a different approach [43

43. N. A. Rejesh, H. Pullagurla, and M. Pramanik, “Deconvolution based deblurring of reconstructed images in photoacoustic/thermoacoustic tomography,” J. Opt. Soc. Am. A 30, 1994–2001 (2013). [CrossRef]

]. There, the PA signal itself has been corrected for blurring using Tikhonov regularization. Such methods can also be used on the recorded PA signals and then can be fed to the reconstruction algorithm, where deblurring due to regularization can be corrected using our proposed technique.

7. Conclusions

Photoacoustic imaging has been shown to be an invaluable imaging modality both in pre-clinical and clinical settings. The quantitative accuracy of the photoacoustic images can be considerably improved with the usage of model-based reconstruction techniques, especially in cases of limited data. These model-based reconstruction techniques utilizes regularization, typically 2-norm based, to obtain more meaningful reconstruction results, with a caveat of smoothening (blurring) the solutions. This work utilized the least-squares QR-based framework to build a approximate model-resolution (blur) matrix, which in turn got utilized in basis pursuit deconvolution approach to restore the sharp features of the reconstructed photoacoustic image. More importantly, it was shown that the proposed methodology works well when the regularization parameter is chosen within reasonable bounds (rather than at only optimal value), reducing the computational overhead considerably. It was also shown that using the proposed framework, the quantitative accuracy of the reconstructed photoacoustic image has improved by more than 50%.

Acknowledgments

This work is supported by Department of Biotechnology (DBT) Rapid Grant for Young investigator (RGYI) (No: BT/PR6494/GDB/27/415/2012) and DBT Bioengineering Grant (No: BT/PR7994/MED/32/284/2013). JP acknowledges the support by Microsoft Corporation and Microsoft Research India under the Microsoft Research India PhD Fellowship Award and SPIE Optics and Photonics Education Scholarship.

References and links

1.

V. E Gusev and A. A Karabutov, Laser Optoacoustics (AIP, 1993).

2.

A. A. Karabutov, E. V. Savateeva, and A. A. Oraevsky, “Optoacoustic tomography: New modality of laser diagnostic systems,” Laser Phys. 13, 711–723 (2003).

3.

L. H. V. Wang and S. Hu, “Photoacoustic Tomography: In Vivo Imaging from Organelles to Organs,” Science 335, 1458–1462, (2012). [CrossRef] [PubMed]

4.

H. F. Zhang, K. Maslov, G. Stoica, and L. H. V. Wang, “Functional photoacoustic microscopy for high-resolution and noninvasive in vivo imaging,” Nat. Biotech. 24, 848–851 (2006). [CrossRef]

5.

B. Yin, D. Xing, Y. Wang, Y. Zeng, Y. Tan, and Q. Chen, “Fast photoacoustic imaging system based on 320-element linear transducer array,” Phys. Med. Bio. 49, 1339–1346 (2004). [CrossRef]

6.

G. Ku and L. H. V. Wang, “Scanning thermoacoustic tomography in biological tissue,” Med. Phys. 27, 1195–1202 (2000). [CrossRef] [PubMed]

7.

M. Pramanik, G. Ku, C. H. Li, and L. H. V. Wang, “Design and evaluation of a novel breast cancer detection system combining both thermoacoustic (TA) and photoacoustic (PA) tomography,” Med. Phys. 35, 2218–2223 (2008). [CrossRef] [PubMed]

8.

M. Pramanik and L. H. V. Wang, “Thermoacoustic and photoacoustic sensing of temperature,” J. Biomed. Opt. 14, 054024 (2009). [CrossRef] [PubMed]

9.

K. H. Song and L. H. V. Wang, “Noninvasive photoacoustic imaging of the thoracic cavity and the kidney in small and large animals,” Med. Phys. 35, 4524–4529 (2008). [CrossRef] [PubMed]

10.

Y. Wang, X. Y. Xie, X. D. Wang, G. Ku, K. L. Gill, D. P. O’Neal, G. Stoica, and L. H. V. Wang, “Photoacoustic Tomography of a Nanoshell Contrast Agent in the in Vivo Rat Brain,” Nano Lett. 4, 1689–1692 (2004). [CrossRef]

11.

A. De la Zerda, C. Zavaleta, S. Keren, S. Vaithilingam, S. Bodapati, Z. Liu, J. Levi, B. R. Smith, T. J. Ma, O. Oralkan, Z. Cheng, X. Y. Chen, H. J. Dai, B. T. Khuri-Yakub, and S. S. Gambhir, “Carbon nanotubes as photoacoustic molecular imaging agents in living mice,” Nat. Nanotech. 3, 557–562 (2008). [CrossRef]

12.

P. Kuchment and L. Kunyansky, “Mathematics of thermoacoustic and photoacoustic tomography,” European J. App. Math. 19, 191–224 (2008).

13.

M. Xu and L. H. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Review E 71, 016706 (2005). [CrossRef]

14.

K. P. Kostli and P. C. Beard, “Two-dimensional photoacoustic imaging by use of Fourier transform image reconstruction and a detector with an anisotropic response,” App. Opt. 42, 1899–1908 (2003). [CrossRef]

15.

P. Ephrat, L. Keenliside, A. Seabrook, F. S. Prato, and J. J. L. Carson, “Three-dimensional photoacoustic imaging by sparse-array detection and iterative image reconstruction,” J. Biomed. Opt. 13, 054052 (2008). [CrossRef] [PubMed]

16.

M. A. Anastasio, J. Zhang, X. Pan, Y. Zou, G. Ku, and L. H. V. Wang, “Half-time image reconstruction in thermoacoustic tomography,” IEEE Trans. Med. Imag. 24, 199–210 (2005). [CrossRef]

17.

K. Wang, R. Su, A. A. Oraevsky, and M. A. Anastasio, “Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography,” Phys. Med. Bio. 57, 5399–5423 (2012). [CrossRef]

18.

C. Huang, K. Wang, L. Nie, L. V. Wang, and M. A. Anastasio, “Full-Wave Iterative Image Reconstruction in Photoacoustic Tomography with Acoustically Inhomogeneous Media,” IEEE Trans. Med. Imag. 32, 1097–1110 (2013). [CrossRef]

19.

X. L. Dean-Ben, V. Ntziachristos, and D. Razansky, “Acceleration of Optoacoustic Model-Based Reconstruction Using Angular Image Discretization,” IEEE Trans. Med. Imag. 31, 1154–1162 (2012). [CrossRef]

20.

X. L. Dean-Ben, A. Buehler, V. Ntziachristos, and D. Razansky, “Accurate Model-Based Reconstruction Algorithm for Three-Dimensional Optoacoustic Tomography,” IEEE Trans. Med. Imag. 31, 1922–1928 (2012). [CrossRef]

21.

A. Buehler, A. Rosenthal, T. Jetzfellner, A. Dima, D. Razansky, and V. Ntziachristos, “Model-based optoacoustic inversions with incomplete projection data,” Med. Phys. 38, 1694–1704 (2011). [CrossRef] [PubMed]

22.

C. B. Shaw, J. Prakash, M. Pramanik, and P. K. Yalavarthy, “Least Squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography,” J. Biomed. Opt. 18, 080501:1–3 (2013).

23.

J. Prakash and P. K. Yalavarthy, “A LSQR-type method provides a computationally efficient automated optimal choice of regularization parameter in diffuse optical tomography,” Med. Phys. 40, 033101 (2013). [CrossRef] [PubMed]

24.

J. A. Fessler and W. L. Rogers, “Spatial resolution properties of penalized- likelihood image reconstruction methods: Space-invariant tomographs,” IEEE Trans. Image Process. 5, 1346–1358 (1996). [CrossRef]

25.

J. Prakash, H. Dehghani, B. W. Pogue, and P. K. Yalavarthy, “Model-Resolution based Basis Pursuit Deconvolution Improves Diffuse Optical Tomographic Imaging,” IEEE Trans. Med. Imag. 2014 (in press). [CrossRef]

26.

J. Chen, R. Lin, H. Wang, J. Meng, H. Zheng, and L. Song, “Blind-deconvolution optical-resolution photoacoustic microscopy in vivo,” Opt. Express 21, 7316–7327 (2013). [CrossRef] [PubMed]

27.

B. R. Hunt, “The application of constrained least square estimation to image restoration by digital computer,” IEEE Trans. Comput. C-22, 805–812 (1973). [CrossRef]

28.

E. Candes and M. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag. 25, 21–30 (2008). [CrossRef]

29.

J. Romberg, “Imaging via compressive sampling,” IEEE Signal Process. Mag. 25, 14–20 (2008). [CrossRef]

30.

J. Provost and F. Lesage, “The application of compressed sensing for photo-acoustic tomography,” IEEE Trans. Med. Imag. 28, 585–594 (2009). [CrossRef]

31.

Z. Guo, C. H. Li, L. Song, and L. H. V. Wang, “Compressed sensing in photoacoustic tomography in vivo,” J. Biomed. Opt. 15, 021311 (2010). [CrossRef] [PubMed]

32.

M. Figueiredo, J. Bioucas-Dias, and M. Afonso, “Fast Frame-Based Image Deconvolution Using Variable Splitting and Constrained Optimization,” IEEE Worskhop on Statistical Signal Processing, Cardiff, Wales, 2009.

33.

I. Selesnick, “Introduction to Sparsity in Signal Processing [Connexions Web site],” Available at: http://cnx.org/content/m43545/1.3/. (2012).

34.

Y. Hristova, P. Kuchment, and L. V. Nguyen, “Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media,” Inv. Problems 24, 055006 (2008). [CrossRef]

35.

B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt. 15, 021314 (2010). [CrossRef] [PubMed]

36.

G. Paltauf, J. A. Viator, S. A. Prahl, and S. L. Jacques, “Iterative reconstruction algorithm for optoacoustic imaging,” J. Acous. Soc. Am. 112, 1536–1544 (2002). [CrossRef]

37.

C. C. Paige and M. A. Saunders, “LSQR: An algorithm for sparse linear equations and sparse least squares,” ACM Trans. Math. Software 8, 43–71 (1982). [CrossRef]

38.

M. E. Kilmer and D. P. OLeary, “Choosing regularization parameters in iterative methods for ill-posed problems,” SIAM J. Matrix Anal. Appl. 22, 1204–1221 (2001). [CrossRef]

39.

P. C. Hansen, “Regularization tools version 4.0 for Matlab 7.3,” Numerical Algorithms 46, 189–194 (2007). [CrossRef]

40.

P. C. Hansen, J. G. Nagy, and D. P. O. Leary, Deblurring Images: Matrices, Spectra, and Filtering, 1 (SIAM, Philadelphia, 2006). [CrossRef]

41.

M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo, “Fast Image Recovery Using Variable Splitting and Constrained Optimization,” IEEE Trans. Image Process. 19, 2345–2356 (2010). [CrossRef] [PubMed]

42.

J. Prakash, C. B. Shaw, R. Manjappa, R. Kanhirodan, and P. K. Yalavarthy, “Sparse Recovery Methods Hold Promise for Diffuse Optical Tomographic Image Reconstruction,” IEEE J. Sel. Top. Quantum Electron. 20, 6800609 (2014). [CrossRef]

43.

N. A. Rejesh, H. Pullagurla, and M. Pramanik, “Deconvolution based deblurring of reconstructed images in photoacoustic/thermoacoustic tomography,” J. Opt. Soc. Am. A 30, 1994–2001 (2013). [CrossRef]

OCIS Codes
(170.0170) Medical optics and biotechnology : Medical optics and biotechnology
(170.3010) Medical optics and biotechnology : Image reconstruction techniques
(170.5120) Medical optics and biotechnology : Photoacoustic imaging

ToC Category:
Image Reconstruction and Inverse Problems

History
Original Manuscript: January 6, 2014
Revised Manuscript: February 18, 2014
Manuscript Accepted: March 17, 2014
Published: April 2, 2014

Citation
Jaya Prakash, Aditi Subramani Raju, Calvin B. Shaw, Manojit Pramanik, and Phaneendra K. Yalavarthy, "Basis pursuit deconvolution for improving model-based reconstructed images in photoacoustic tomography," Biomed. Opt. Express 5, 1363-1377 (2014)
http://www.opticsinfobase.org/boe/abstract.cfm?URI=boe-5-5-1363


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. V. E Gusev and A. A Karabutov, Laser Optoacoustics (AIP, 1993).
  2. A. A. Karabutov, E. V. Savateeva, and A. A. Oraevsky, “Optoacoustic tomography: New modality of laser diagnostic systems,” Laser Phys.13, 711–723 (2003).
  3. L. H. V. Wang and S. Hu, “Photoacoustic Tomography: In Vivo Imaging from Organelles to Organs,” Science335, 1458–1462, (2012). [CrossRef] [PubMed]
  4. H. F. Zhang, K. Maslov, G. Stoica, and L. H. V. Wang, “Functional photoacoustic microscopy for high-resolution and noninvasive in vivo imaging,” Nat. Biotech.24, 848–851 (2006). [CrossRef]
  5. B. Yin, D. Xing, Y. Wang, Y. Zeng, Y. Tan, and Q. Chen, “Fast photoacoustic imaging system based on 320-element linear transducer array,” Phys. Med. Bio.49, 1339–1346 (2004). [CrossRef]
  6. G. Ku and L. H. V. Wang, “Scanning thermoacoustic tomography in biological tissue,” Med. Phys.27, 1195–1202 (2000). [CrossRef] [PubMed]
  7. M. Pramanik, G. Ku, C. H. Li, and L. H. V. Wang, “Design and evaluation of a novel breast cancer detection system combining both thermoacoustic (TA) and photoacoustic (PA) tomography,” Med. Phys.35, 2218–2223 (2008). [CrossRef] [PubMed]
  8. M. Pramanik and L. H. V. Wang, “Thermoacoustic and photoacoustic sensing of temperature,” J. Biomed. Opt.14, 054024 (2009). [CrossRef] [PubMed]
  9. K. H. Song and L. H. V. Wang, “Noninvasive photoacoustic imaging of the thoracic cavity and the kidney in small and large animals,” Med. Phys.35, 4524–4529 (2008). [CrossRef] [PubMed]
  10. Y. Wang, X. Y. Xie, X. D. Wang, G. Ku, K. L. Gill, D. P. O’Neal, G. Stoica, and L. H. V. Wang, “Photoacoustic Tomography of a Nanoshell Contrast Agent in the in Vivo Rat Brain,” Nano Lett.4, 1689–1692 (2004). [CrossRef]
  11. A. De la Zerda, C. Zavaleta, S. Keren, S. Vaithilingam, S. Bodapati, Z. Liu, J. Levi, B. R. Smith, T. J. Ma, O. Oralkan, Z. Cheng, X. Y. Chen, H. J. Dai, B. T. Khuri-Yakub, and S. S. Gambhir, “Carbon nanotubes as photoacoustic molecular imaging agents in living mice,” Nat. Nanotech.3, 557–562 (2008). [CrossRef]
  12. P. Kuchment and L. Kunyansky, “Mathematics of thermoacoustic and photoacoustic tomography,” European J. App. Math.19, 191–224 (2008).
  13. M. Xu and L. H. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Review E71, 016706 (2005). [CrossRef]
  14. K. P. Kostli and P. C. Beard, “Two-dimensional photoacoustic imaging by use of Fourier transform image reconstruction and a detector with an anisotropic response,” App. Opt.42, 1899–1908 (2003). [CrossRef]
  15. P. Ephrat, L. Keenliside, A. Seabrook, F. S. Prato, and J. J. L. Carson, “Three-dimensional photoacoustic imaging by sparse-array detection and iterative image reconstruction,” J. Biomed. Opt.13, 054052 (2008). [CrossRef] [PubMed]
  16. M. A. Anastasio, J. Zhang, X. Pan, Y. Zou, G. Ku, and L. H. V. Wang, “Half-time image reconstruction in thermoacoustic tomography,” IEEE Trans. Med. Imag.24, 199–210 (2005). [CrossRef]
  17. K. Wang, R. Su, A. A. Oraevsky, and M. A. Anastasio, “Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography,” Phys. Med. Bio.57, 5399–5423 (2012). [CrossRef]
  18. C. Huang, K. Wang, L. Nie, L. V. Wang, and M. A. Anastasio, “Full-Wave Iterative Image Reconstruction in Photoacoustic Tomography with Acoustically Inhomogeneous Media,” IEEE Trans. Med. Imag.32, 1097–1110 (2013). [CrossRef]
  19. X. L. Dean-Ben, V. Ntziachristos, and D. Razansky, “Acceleration of Optoacoustic Model-Based Reconstruction Using Angular Image Discretization,” IEEE Trans. Med. Imag.31, 1154–1162 (2012). [CrossRef]
  20. X. L. Dean-Ben, A. Buehler, V. Ntziachristos, and D. Razansky, “Accurate Model-Based Reconstruction Algorithm for Three-Dimensional Optoacoustic Tomography,” IEEE Trans. Med. Imag.31, 1922–1928 (2012). [CrossRef]
  21. A. Buehler, A. Rosenthal, T. Jetzfellner, A. Dima, D. Razansky, and V. Ntziachristos, “Model-based optoacoustic inversions with incomplete projection data,” Med. Phys.38, 1694–1704 (2011). [CrossRef] [PubMed]
  22. C. B. Shaw, J. Prakash, M. Pramanik, and P. K. Yalavarthy, “Least Squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography,” J. Biomed. Opt.18, 080501:1–3 (2013).
  23. J. Prakash and P. K. Yalavarthy, “A LSQR-type method provides a computationally efficient automated optimal choice of regularization parameter in diffuse optical tomography,” Med. Phys.40, 033101 (2013). [CrossRef] [PubMed]
  24. J. A. Fessler and W. L. Rogers, “Spatial resolution properties of penalized- likelihood image reconstruction methods: Space-invariant tomographs,” IEEE Trans. Image Process.5, 1346–1358 (1996). [CrossRef]
  25. J. Prakash, H. Dehghani, B. W. Pogue, and P. K. Yalavarthy, “Model-Resolution based Basis Pursuit Deconvolution Improves Diffuse Optical Tomographic Imaging,” IEEE Trans. Med. Imag. 2014 (in press). [CrossRef]
  26. J. Chen, R. Lin, H. Wang, J. Meng, H. Zheng, and L. Song, “Blind-deconvolution optical-resolution photoacoustic microscopy in vivo,” Opt. Express21, 7316–7327 (2013). [CrossRef] [PubMed]
  27. B. R. Hunt, “The application of constrained least square estimation to image restoration by digital computer,” IEEE Trans. Comput.C-22, 805–812 (1973). [CrossRef]
  28. E. Candes and M. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag.25, 21–30 (2008). [CrossRef]
  29. J. Romberg, “Imaging via compressive sampling,” IEEE Signal Process. Mag.25, 14–20 (2008). [CrossRef]
  30. J. Provost and F. Lesage, “The application of compressed sensing for photo-acoustic tomography,” IEEE Trans. Med. Imag.28, 585–594 (2009). [CrossRef]
  31. Z. Guo, C. H. Li, L. Song, and L. H. V. Wang, “Compressed sensing in photoacoustic tomography in vivo,” J. Biomed. Opt.15, 021311 (2010). [CrossRef] [PubMed]
  32. M. Figueiredo, J. Bioucas-Dias, and M. Afonso, “Fast Frame-Based Image Deconvolution Using Variable Splitting and Constrained Optimization,” IEEE Worskhop on Statistical Signal Processing, Cardiff, Wales, 2009.
  33. I. Selesnick, “Introduction to Sparsity in Signal Processing [Connexions Web site],” Available at: http://cnx.org/content/m43545/1.3/ . (2012).
  34. Y. Hristova, P. Kuchment, and L. V. Nguyen, “Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media,” Inv. Problems24, 055006 (2008). [CrossRef]
  35. B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt.15, 021314 (2010). [CrossRef] [PubMed]
  36. G. Paltauf, J. A. Viator, S. A. Prahl, and S. L. Jacques, “Iterative reconstruction algorithm for optoacoustic imaging,” J. Acous. Soc. Am.112, 1536–1544 (2002). [CrossRef]
  37. C. C. Paige and M. A. Saunders, “LSQR: An algorithm for sparse linear equations and sparse least squares,” ACM Trans. Math. Software8, 43–71 (1982). [CrossRef]
  38. M. E. Kilmer and D. P. OLeary, “Choosing regularization parameters in iterative methods for ill-posed problems,” SIAM J. Matrix Anal. Appl.22, 1204–1221 (2001). [CrossRef]
  39. P. C. Hansen, “Regularization tools version 4.0 for Matlab 7.3,” Numerical Algorithms46, 189–194 (2007). [CrossRef]
  40. P. C. Hansen, J. G. Nagy, and D. P. O. Leary, Deblurring Images: Matrices, Spectra, and Filtering, 1 (SIAM, Philadelphia, 2006). [CrossRef]
  41. M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo, “Fast Image Recovery Using Variable Splitting and Constrained Optimization,” IEEE Trans. Image Process.19, 2345–2356 (2010). [CrossRef] [PubMed]
  42. J. Prakash, C. B. Shaw, R. Manjappa, R. Kanhirodan, and P. K. Yalavarthy, “Sparse Recovery Methods Hold Promise for Diffuse Optical Tomographic Image Reconstruction,” IEEE J. Sel. Top. Quantum Electron.20, 6800609 (2014). [CrossRef]
  43. N. A. Rejesh, H. Pullagurla, and M. Pramanik, “Deconvolution based deblurring of reconstructed images in photoacoustic/thermoacoustic tomography,” J. Opt. Soc. Am. A30, 1994–2001 (2013). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited