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

Biomedical Optics Express, Vol. 5, Issue 5, pp. 1363-1377 (2014)

http://dx.doi.org/10.1364/BOE.5.001363

Acrobat PDF (2103 KB)

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

*in vivo*biomedical optical imaging modality that combines both optics and ultrasonic physics [1–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]

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]

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]

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]

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]

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]

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

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]

*l*

_{1}-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. 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]

*Split augmented Lagrangian shrinkage algorithm*(SALSA) is a well known algorithm for performing BPD and the same has been used here [32]. 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).

## 2. Analytical and LSQR based reconstruction for photoacoustic tomography

### 2.1. k*-wave based time reversal (analytical)*

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]

*g*(

*y*,

*t*) is the time series pressure data measured at transducer’s location

*y*∈

*B*at time t. The

*p*and

_{t}*p*) and

*c*(

*r*) represents the speed of the ultrasound propagation (kept constant in here). The Δ

*is the Laplacian function with respect to space*

_{r}*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*).

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]

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]

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

*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]

*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]

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

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]

*n*×

*n*is vectorized by stacking the columns into

*n*

^{2}× 1 vector, represented as

*x*. Each column of the system matrix (

**A**, having dimension of

*m*×

*n*

^{2}) 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], with limited transducer bandwidth.

*k*-Wave MATLAB toolbox [35

**15**, 021314 (2010). [CrossRef] [PubMed]

*k*-space pseudo spectral solution method with two coupled first order equations [35

**15**, 021314 (2010). [CrossRef] [PubMed]

*x*is a long column vector having a dimension of

*n*

^{2}× 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]

*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)

*λ*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

*x*, resulting in, 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

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]

*x*

^{(}

^{k}^{)}being the dimensionality reduced version of

*x*with

*k*<

*n*

^{2}[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]

**U**

_{k+1}

^{T}**U**=

_{k+1}**I**results in, Considering the first order condition of the Eq. 11, the new update equation becomes [22, 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]

*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]

*k*=

*n*

^{2},

*x*=

_{LSQR}*x*for a given value of

_{Tikh}*λ*. For all other cases,

*x*becomes an approximation to

_{LSQR}*x*. It is well-known that there exists an optimal

_{Tikh}*k*with

*k*≪

*n*

^{2}, for which solutions

*x*and

_{LSQR}*x*are close within the precision limits of the computation [23

_{Tikh}**40**, 033101 (2013). [CrossRef] [PubMed]

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

**V**) 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 (

_{k}*λ*) was previously proposed in Ref. [22, 23

**40**, 033101 (2013). [CrossRef] [PubMed]

*Ax*−

_{LSQR}*b*||

_{2}. The estimation of the number of Lanczos iterations was also performed automatically in Ref. [22, 23

**40**, 033101 (2013). [CrossRef] [PubMed]

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

*x*

^{(}

^{k}^{)}with its estimated version (

*λ*) is equal to zero, then only

*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

**M**) will not be possible (

**B**

_{k}

^{T}**B**is a rank-deficient highly ill-conditioned matrix).

_{k}### 3.2. Basis pursuit deconvolution in photoacoustic tomography

40. P. C. Hansen, J. G. Nagy, and D. P. O. Leary, *Deblurring Images: Matrices, Spectra, and Filtering*, 1 (SIAM, Philadelphia, 2006). [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]

*ℓ*

_{1}-norm, which promotes sparseness and sharp features [28

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]

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]

**M**is the model resolution matrix derived using the least squares QR (LSQR) algorithm (Eq. 19), and

*λ*

_{l}_{1}) is the deblurred estimate of the

*x̃*is the final output (estimate of deblurred

_{LSQR}*x*). It is possible to get exact model-resolution matrix using the Tikhonov formulation (which is computationally expensive) [25

_{LSQR}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]

*x*

^{(}

^{k}^{)}compared to

*Spilt augmented lagrangian shrinkage algorithm*(SALSA) [32, 41

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]

**M**) dimensionality is small, among all existing

*ℓ*

_{1}-norm based algorithms achieved via the usage of variable splitting in the minimization problem. This method uses an alternating direction method of multipliers (ADMM), which is based on the augmented Lagrangian method (ALM) [32, 41

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]

*λ*

_{l}_{1}) is chosen as 0.00001 for all cases. The ADMM parameter (

*d*), which has similar size of

*α*in the SALSA algorithm weighs the

*ℓ*

_{2}-norm of the unknown parameter, which was chosen to be 0.1 for SNR of the data being 40 dB and 10 for all other SNR values. The number of iterations

*N*is set to 5000. In here,

_{it}*N*was chosen to be a large value beyond which the solution did not improve. Note that the computational time for each of these sub-iterations is in the order of milli-seconds. The

_{it}**M**and

*ℓ*

_{1}-norm of

*ℓ*

_{1}-norm penalty from being null, when

*λ*

_{l}_{1})/

*α*. The step-3 in here is a direct translation of Tikhonov-kind of update, which utilizes the normal equations. The step-4 updates the ADMM parameter, ideally when converged to solution there is no update in the

*d*[32, 41

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]

## 4. Figures of merit

### 4.1. Pearson Correlation (PC)

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]

*x*is the expected initial pressure distribution and

*x*is the reconstructed initial pressure distribution using one of the above explained methods. Here

^{recon}*σ*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)

*μ*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

## 5. Numerical experiments

*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

**15**, 021314 (2010). [CrossRef] [PubMed]

## 6. Results and discussion

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

*λ*= 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).

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

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

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

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

*λ*) as long as it is chosen within the reasonable bounds. These bounds are the values of

*λ*for which where

*λ*. In simple words, unless

*x*

^{(k)}, this relation will not hold good, which essentially means that the BPD methods works well for all cases where

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]

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]

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]

## 7. Conclusions

*ℓ*

_{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

## References and links

1. | V. E Gusev and A. A Karabutov, |

2. | A. A. Karabutov, E. V. Savateeva, and A. A. Oraevsky, “Optoacoustic tomography: New modality of laser diagnostic systems,” Laser Phys. |

3. | L. H. V. Wang and S. Hu, “Photoacoustic Tomography: In Vivo Imaging from Organelles to Organs,” Science |

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

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

6. | G. Ku and L. H. V. Wang, “Scanning thermoacoustic tomography in biological tissue,” Med. Phys. |

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

8. | M. Pramanik and L. H. V. Wang, “Thermoacoustic and photoacoustic sensing of temperature,” J. Biomed. Opt. |

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

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

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

12. | P. Kuchment and L. Kunyansky, “Mathematics of thermoacoustic and photoacoustic tomography,” European J. App. Math. |

13. | M. Xu and L. H. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Review E |

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

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

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

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

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

19. | X. L. Dean-Ben, V. Ntziachristos, and D. Razansky, “Acceleration of Optoacoustic Model-Based Reconstruction Using Angular Image Discretization,” IEEE Trans. Med. Imag. |

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

21. | A. Buehler, A. Rosenthal, T. Jetzfellner, A. Dima, D. Razansky, and V. Ntziachristos, “Model-based optoacoustic inversions with incomplete projection data,” Med. Phys. |

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

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

24. | J. A. Fessler and W. L. Rogers, “Spatial resolution properties of penalized- likelihood image reconstruction methods: Space-invariant tomographs,” IEEE Trans. Image Process. |

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 |

27. | B. R. Hunt, “The application of constrained least square estimation to image restoration by digital computer,” IEEE Trans. Comput. |

28. | E. Candes and M. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag. |

29. | J. Romberg, “Imaging via compressive sampling,” IEEE Signal Process. Mag. |

30. | J. Provost and F. Lesage, “The application of compressed sensing for photo-acoustic tomography,” IEEE Trans. Med. Imag. |

31. | Z. Guo, C. H. Li, L. Song, and L. H. V. Wang, “Compressed sensing in photoacoustic tomography in vivo,” J. Biomed. Opt. |

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 |

35. | B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt. |

36. | G. Paltauf, J. A. Viator, S. A. Prahl, and S. L. Jacques, “Iterative reconstruction algorithm for optoacoustic imaging,” J. Acous. Soc. Am. |

37. | C. C. Paige and M. A. Saunders, “LSQR: An algorithm for sparse linear equations and sparse least squares,” ACM Trans. Math. Software |

38. | M. E. Kilmer and D. P. OLeary, “Choosing regularization parameters in iterative methods for ill-posed problems,” SIAM J. Matrix Anal. Appl. |

39. | P. C. Hansen, “Regularization tools version 4.0 for Matlab 7.3,” Numerical Algorithms |

40. | P. C. Hansen, J. G. Nagy, and D. P. O. Leary, |

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

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

43. | N. A. Rejesh, H. Pullagurla, and M. Pramanik, “Deconvolution based deblurring of reconstructed images in photoacoustic/thermoacoustic tomography,” J. Opt. Soc. Am. A |

**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: Year | Journal | Reset

### References

- V. E Gusev and A. A Karabutov, Laser Optoacoustics (AIP, 1993).
- A. A. Karabutov, E. V. Savateeva, and A. A. Oraevsky, “Optoacoustic tomography: New modality of laser diagnostic systems,” Laser Phys.13, 711–723 (2003).
- L. H. V. Wang and S. Hu, “Photoacoustic Tomography: In Vivo Imaging from Organelles to Organs,” Science335, 1458–1462, (2012). [CrossRef] [PubMed]
- 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]
- 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]
- G. Ku and L. H. V. Wang, “Scanning thermoacoustic tomography in biological tissue,” Med. Phys.27, 1195–1202 (2000). [CrossRef] [PubMed]
- 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]
- M. Pramanik and L. H. V. Wang, “Thermoacoustic and photoacoustic sensing of temperature,” J. Biomed. Opt.14, 054024 (2009). [CrossRef] [PubMed]
- 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]
- 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]
- 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]
- P. Kuchment and L. Kunyansky, “Mathematics of thermoacoustic and photoacoustic tomography,” European J. App. Math.19, 191–224 (2008).
- M. Xu and L. H. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Review E71, 016706 (2005). [CrossRef]
- 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]
- 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]
- 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]
- 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]
- 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]
- 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]
- 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]
- 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]
- 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).
- 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]
- 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]
- 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]
- 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]
- 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]
- E. Candes and M. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag.25, 21–30 (2008). [CrossRef]
- J. Romberg, “Imaging via compressive sampling,” IEEE Signal Process. Mag.25, 14–20 (2008). [CrossRef]
- J. Provost and F. Lesage, “The application of compressed sensing for photo-acoustic tomography,” IEEE Trans. Med. Imag.28, 585–594 (2009). [CrossRef]
- 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]
- 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.
- I. Selesnick, “Introduction to Sparsity in Signal Processing [Connexions Web site],” Available at: http://cnx.org/content/m43545/1.3/ . (2012).
- 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]
- 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]
- 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]
- 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]
- 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]
- P. C. Hansen, “Regularization tools version 4.0 for Matlab 7.3,” Numerical Algorithms46, 189–194 (2007). [CrossRef]
- P. C. Hansen, J. G. Nagy, and D. P. O. Leary, Deblurring Images: Matrices, Spectra, and Filtering, 1 (SIAM, Philadelphia, 2006). [CrossRef]
- 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]
- 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]
- 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.