## 1.

## Introduction

The field of photoacoustic tomography (PAT) has been expanding rapidly in the past few years.^{1} By combining strong optical absorption contrast and high ultrasonic resolution in a single modality, PAT can achieve much better spatial resolution at depths beyond the optical ballistic regime (
$\sim 1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
in the skin) than the traditional optical modalities.^{2, 3} In PAT, biological tissues are usually irradiated by a pulsed laser. Absorbed energy is converted heat, which is further converted to a pressure rise via thermoelastic expansion. The initial pressure rise then propagates as ultrasonic waves, which are detected by ultrasound sensors, and the received ultrasonic signals are used to form an image. When the excitation laser is replaced by microwave or RF sources, the technique is called thermoacoustic tomography (TAT).^{4, 5} Both PAT and TAT have been used successfully in a variety of applications, including high-quality *in vivo* vascular structural imaging, hemodynamic functional imaging,^{6, 7} visualization of breast tumors,^{8, 9} and molecular imaging of biomarkers with exogenous contrast agents.^{10, 11, 12}

PAT has been implemented in various forms, and each form has its own advantages and applications.^{1} In this paper, we focus on photoacoustic computed tomography (PACT, or simply PAT), in which an array of unfocused ultrasonic transducers is placed outside the object, and an inverse algorithm is used to reconstruct the image. Closed-form reconstruction formulas have been reported in both the frequency and time domains for spherical, planar, and cylindrical detecting geometries.^{13, 14, 15, 16, 17, 18, 19, 20} However, a fundamental assumption of all these algorithms is that the spatial sampling of the detecting aperture is sufficient; otherwise, undersampling artifacts, such as streaking artifacts or grating lobes, appear.

Reliable image reconstruction with sparse sampling of the detecting aperture is desirable. In practical PAT systems, it is recommended^{1, 21} to set the discrete spatial sampling period to be two to five times smaller than the sensing aperture of the detector. For a scanning PAT system, it may require hundreds or even thousands of scanning steps to acquire an image, depending on the sizes of both the detector and the detecting aperture. Such scanning usually takes several minutes to complete. To reach real-time imaging, PAT is implemented with an array of ultrasonic transducers, where all or groups of the array elements can detect photoacoustic signals simultaneously. However, the data acquisition speed is still limited by the number of parallel data acquisition (DAQ) channels, and employing a large number of DAQ channels greatly increases the system cost. For example, for a fast 512-element ring array PAT system with a 64-channel data acquisition module,^{22} it takes 8 laser shots to collect data from all 512 elements. For direct 3-D reconstruction PAT applications,^{23, 24} the data from a 2-D ultrasonic array is usually an extremely sparse sampling of the detecting aperture. Moreover, channel cross talk is also related to the space between neighboring elements (kerf), and an extensive spatial sampling may increase the cross talk.

Imaging an object in PAT can be understood as sensing the object in a certain domain. For example, with the Fourier-shell identity,^{25} PAT can also be seen as detecting the spatial frequencies of the object (sensing in the Fourier domain). Sparse spatial sampling of the detection aperture implies that the spatial frequencies cannot be exactly determined. Traditional backprojection (BP) reconstruction methods^{16} reconstruct the image of “minimal energy” under the observation constraints. An improved reconstruction algorithm should be able to “guess” these undetermined frequency components. However, interpolation in the Fourier domain is a critical issue and usually creates artifacts in reconstructed images.^{26} The recently developed compressed sensing (CS) theory^{27} enables us to recover these unobserved components under certain conditions. The theory has been successfully applied in MRI,^{28} where MRI images were able to be reconstructed from significantly undersampled K-space measurements. Reference 29 introduced the CS theory into the field of PAT, and the idea was tested in phantoms using a circular scanning PAT system. In this paper, we improve the speed of the reconstruction algorithm by adopting a nonlinear conjugate gradient descent method and demonstrate the algorithm with both phantom and animal data, using various detecting geometries.

## 2.

## Method

## 2.1.

### Forward and Inverse Problems in PAT and Their Numerical Implementations

In PAT, pulsed laser irradiation creates pressure rises as a result of thermoelastic expansion. These initial pressure rises propagate as photoacoustic waves, which can be detected by ultrasonic sensors. Based on the pressure measurements $p(\stackrel{\u20d1}{r},t)$ at the detecting aperture, PAT tries to reconstruct the image of the initial pressure rise distribution ${p}_{0}\left(\stackrel{\u20d1}{r}\right)$ . The forward and inverse problems in PAT express the reciprocal relationship between ${p}_{0}\left(\stackrel{\u20d1}{r}\right)$ and $p(\stackrel{\u20d1}{r},t)$ . By solving the wave equation, the forward problem, which predicts $p(\stackrel{\u20d1}{r},t)$ by ${p}_{0}\left(\stackrel{\u20d1}{r}\right)$ , can be derived as (assuming a delta pulse heating):

## Eq. 1

$$p(\stackrel{\u20d1}{r},t)=\frac{\partial}{\partial t}[\frac{1}{4\pi {c}^{3}t}\int \phantom{\rule{0.2em}{0ex}}\mathrm{d}{\stackrel{\u20d1}{r}}^{\prime}{p}_{0}\left({\stackrel{\u20d1}{r}}^{\prime}\right)\delta (t-\frac{|\stackrel{\u20d1}{r}-{\stackrel{\u20d1}{r}}^{\prime}|}{c})],$$^{2}Sometimes, the velocity potential $\phi (\stackrel{\u20d1}{r},t)={\int}_{0}^{t}p(\stackrel{\u20d1}{r},{t}^{\prime})\mathrm{d}{t}^{\prime}$ is employed to simplify Eq. 1:

The analytical inversion of Eq. 1 describes the inverse problem, which reconstructs ${p}_{0}\left(\stackrel{\u20d1}{r}\right)$ with $p(\stackrel{\u20d1}{r},t)$ :

## Eq. 2

$${p}_{0}\left(\stackrel{\u20d1}{r}\right)={\int}_{{S}_{0}}{\phantom{|}[2p({\stackrel{\u20d1}{r}}_{0},\overline{t})-2\overline{t}\partial p({\stackrel{\u20d1}{r}}_{0},\overline{t})\u2215\partial \overline{t}]|}_{\overline{t}=|r-{r}_{0}|}\mathrm{d}{\Omega}_{0}\u2215{\Omega}_{0},$$^{2}

To numerically model the forward and inverse problems, we need to properly discretize Eqs. 1, 2. We use a vector $x$ to represent ${p}_{0}\left(\stackrel{\u20d1}{r}\right)$ , where each element of $x$ is the average value of initial pressure per unit volume. The size of $x$ $({N}_{x}\times {N}_{y}\times {N}_{z})$ depends on the field of view (FOV) and the desired spatial resolution of the reconstructed image. We use a vector $y$ to denote the velocity potential $\phi (\stackrel{\u20d1}{r},t)$ measured by all elements of the sensor array as a function of time. The size of $y$ is the number of detecting positions $\left(L\right)$ times the number of temporal samples at each position $\left(M\right)$ . Then, the forward problem can be described as $y=\Phi x$ , where matrix $\Phi $ is the projection matrix. Similarly, the inverse problem can be written as $\overline{x}={\Phi}^{-1}y$ , where ${\Phi}^{-1}$ represents the inverse process, and $\overline{x}$ is the reconstructed image.

$\Phi $ and ${\Phi}^{-1}$ are usually extremely large matrices (each containing ${N}_{x}\times {N}_{y}\times {N}_{z}\times L\times M$ data points), even for 2-D reconstruction problems. For example, when reconstructing a $256\times 256$ image with measurements from 512 detecting positions, where each position has 1024 time points, both $\Phi $ and ${\Phi}^{-1}$ contain $3.436\times {10}^{10}$ data points ( $\sim 256\phantom{\rule{0.3em}{0ex}}\mathrm{GB}$ if each point is expressed in double-precision), which makes direct matrix operations computationally impractical. Reference 30 tried to solve this problem by saving only nonzero elements of metrics $\Phi $ and ${\Phi}^{-1}$ . For each detecting element $i$ $(i=\mathrm{1,2},\dots ,L)$ , the forward and inverse operations are performed using a matrix of the same size as $x$ . Each element in this matrix stores an index of the temporal sample $(k=\mathrm{1,2},\dots ,M)$ of measurement $i$ , and this index indicates where the corresponding element of $x$ should be projected. As a result, the storage space for both the forward and inverse operations for all elements is reduced to ${N}_{x}\times {N}_{y}\times {N}_{z}\times L$ ( $\sim 256\phantom{\rule{0.3em}{0ex}}\mathrm{MB}$ for the preceding example). To fully take advantage of parallel computing capability, the responses of all the elements can be calculated simultaneously.

## 2.2.

### Compressed Sensing for PAT

If the measurement is incomplete, matrix $\Phi $ is ill-conditioned, and ${\Phi}^{-1}$ is obviously not an exact inversion of $\Phi $ . Intuitively, an incomplete dataset usually leads to uncertainties in the recovery of the signals. In the case of PAT reconstruction with insufficient measurements, the BP method usually generates streaking artifacts or grating lobes. However, these uncertainties can be eliminated by incorporating prior information, such as sparsity constraints. The CS theory was rigorously formulated to reconstruct images from incomplete datasets. To make this possible, the CS theory relies on two principles: sparsity, which pertains to the object of interest, and incoherence, which pertains to the sensing modality. Moreover, a nonlinear reconstruction is used to enforce both sparsity of the image representation and consistency with the acquired data.

## 2.2.1.

#### Sparsity of PAT images

Unlike ultrasound imaging and all other coherent imaging technologies, PAT is devoid of speckle artifacts and sensitive to boundaries because of its optical absorption contrast.^{31} Therefore, computing the finite difference (FD) of PAT images in the spatial domain sometimes directly results in a sparse representation. When imaging objects with rich absorbing structures such as blood vessels in the mouse brain cortex, however, PAT images may not be sparse in the spatial domain. In these cases, we need to project the images onto an appropriate basis set, such as the wavelet basis. Mathematically speaking, if we use a vector
$x\in {R}^{n}$
to represent an
$n$
-pixel image and
$\Psi $
to denote the wavelet basis set, then
$x$
can be expanded as
$x={\sum}_{i=1}^{n}{a}_{i}{\psi}_{i}$
, where
${a}_{i}=\u27e8x,{\psi}_{i}\u27e9$
is the coefficient sequence of
$x$
. Even when most of the image pixels have nonzero values, the wavelet coefficients may provide a concise representation of the original image: most coefficients are small, and the relatively few large coefficients capture most of the information. The speckle-free nature^{31} of PAT images further reduces the number of significant transform coefficients.

## 2.2.2.

#### Incoherence

Since the object $x$ can be visually losslessly reconstructed with only a few large transform coefficients, the problem of sensing $x$ is equivalent to capturing these large coefficients in the representation domain $\Psi $ . The forward problem of PAT can be seen as projecting the object $x$ to the sensing basis set $\Phi $ , and the measurements are the resulting coefficients. The CS theory requires the two basis sets $\Psi $ and $\Phi $ to be incoherent, i.e., the sensing waveforms should have a dense representation in $\Psi $ . In other words, the undersampled sensing basis $\Phi $ should induce only incoherent artifacts that spread out and appear as random noise in $\Psi $ .

It is difficult to mathematically demonstrate that a physical system satisfies the incoherence condition. The transform point-spread function (TPSF)^{28, 29} was introduced to measure incoherence. Figure 1
illustrates the definition of TPSF in PAT. A wavelet transform is adopted as the sparsifying transform
$\Psi $
, and we assume that a circular detecting aperture is uniformly sampled by multiple ultrasonic point sensors. The
$i$
’th transform coefficient
${e}_{i}$
in the domain
$\Psi $
[Fig. 1] is transformed to the image space [Fig. 1] by the inverse discrete wavelet transform (IDWT). Then, the measurements are generated with the forward operator
$\Phi $
and transformed back to the image space [Fig. 1] with the inverse operator
${\Phi}^{-1}$
. Last, the reconstructed image is again transformed to the sparse domain
$\Psi $
[Fig. 1] with the forward discrete wavelet transform (FDWT). TPSF can be mathematically described as
$\mathit{TPSF}(i,j)={e}_{j}^{*}\Psi {\Phi}^{-1}\Phi {\Psi}^{*}{e}_{i}$
, and it measures the leakage of energy away from the
$i$
’th coefficient to other coefficients. The CS theory requires us to properly choose
$\Phi $
and
$\Psi $
so that these interferences can be minimized and spread out in
$\Psi $
. Readers are referred to Reference 29 for a quantitative comparison of the TPSF maps for various
$\Psi $
in PAT.

## 2.2.3.

#### Reconstruction method

If it satisfies the preceding two conditions, a sparse signal can be accurately recovered from highly incomplete datasets by solving a nonlinear convex optimization problem. We now describe in more detail the CS reconstruction method for PAT. In the CS theory, the reconstruction of image $x$ is obtained by solving the following constrained optimization problem:

## Eq. 3

$$\underset{x}{\mathrm{min}}{\left|\Psi x\right|}_{1}\phantom{\rule{0.3em}{0ex}}\mathrm{s.t.}\phantom{\rule{0.3em}{0ex}}{|\Phi x-y|}_{2}<\epsilon .$$The algorithm is implemented with a nonlinear conjugate gradient descent method,^{32} as detailed in Sec. 5. On a laptop with a dual-core
$2\text{-}\mathrm{GHz}$
CPU and
$3\text{-}\mathrm{GB}$
memory, the calculations usually take less than
$10\phantom{\rule{0.3em}{0ex}}\mathrm{min}$
using MATLAB 2008a.

## 3.

## Results and Discussions

## 3.1.

### Tissue-Mimicking Phantom Experiment

We first demonstrate the CS method using a tissue-mimicking phantom experiment. Tissue phantoms were imaged by scanning a virtual point detector in a setup similar to that of Ref. 33 The PA source contained three black human hair crosses glued on top of optical fibers, with an interval between the hair crosses of about $10\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ . Laser pulses with a repetition rate of $10\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$ were diverged by a ground glass to achieve a relatively uniform illumination. The virtual point detectors evenly scanned the object along a horizontal circle, stopping at 240 points, and the signals were averaged over 20 times at each stop. The total data acquisition time was $8\phantom{\rule{0.3em}{0ex}}\mathrm{min}$ .

Figure 2
shows the reconstruction results with the BP [Figs.
2, 2, 2, 2], the CS [Figs.
2, 2, 4, 4], and the traditional iterative^{30} (IR) [Figs.
2, 2, 2, 2] methods, with 240, 120, 80, and 60 tomographic angles. The images are reconstructed with an FOV of
$30\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\times 15\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
. We can observe that the CS method is clearly superior to the BP and the IR methods. This can be shown by extracting and comparing lines from the reconstructed images [Fig. 2]. The interference level has been reduced significantly with the CS reconstruction. Moreover, as predicted by the theory, the CS scheme is robust to inaccurate measurements, so the noise level has also been suppressed. We took Fig. 2 as the gold standard, and calculated the mean-square errors (MSEs) of all other images from the standard, as shown in Fig. 2. Using the CS reconstruction method, we improved the data acquisition time in the circular scanning geometry by fourfold.

## 3.2.

###
*In Vivo* Experiment

The first *in vivo* experiment was based on a custom-designed 512-element photoacoustic tomography array system.^{22} The
$5\text{-}\mathrm{MHz}$
piezocomposite transducer array was formed into a complete circular aperture. With a 64-channel data acquisition module, the system could provide full tomographic imaging at up to
$8\phantom{\rule{0.3em}{0ex}}\text{frames}\u2215\mathrm{s}$
. We used this system to image mouse cortical blood vessels. The images were reconstructed by the BP [Figs.
3, 3, 3, 3
] and the CS [Figs.
3, 3, 3, 3] algorithms, with 512, 256, 171, and 128 detecting elements, respectively. To achieve the optimal reconstruction results, we simultaneously used both the total variance (TV) and the wavelet as sparsifying transforms in the CS method. The undersampling artifacts appear in the outer region in Fig. 3, which is a natural result of the spatial variant PSF in PAT. Figure 3 shows the images reconstructed with 128 tomographic angles using only the TV regularization, which promotes sharp boundary features and suppresses small variances. Figure 3 shows the images reconstructed using only the wavelet regularization, and the image is “blurred.” Since 128 tomographic angles do not contain enough information to capture all the important transform coefficients, the reconstruction artifacts started to appear and some object features started to disappear. The image reconstruction result with the CS method can be further improved if a more accurate system matrix is employed. For example, unlike the virtual point detectors used in the phantom study, the detecting elements of the circular ultrasonic array system are not perfect point detectors. Therefore, further improvements can be achieved by taking into consideration of the limited detecting angles of the detecting elements.^{25}

The second *in vivo* experiment demonstrates the capability of the CS method with linear array detecting geometry. The
$30\text{-}\mathrm{MHz}$
broadband linear transducer array has a total of 48 elements of dimensions
$82\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}\times 2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
with
$100\text{-}\mu \mathrm{m}$
pitch.^{34} The linear array is focused in the elevation direction to perform cross-sectional (B-scan) imaging, and 3-D volume imaging can be achieved by scanning the probe in the third dimension. We scanned the upper dorsal region of a rat to image the subcutaneous vasculature and acquired a total of 166 B-scan slices. Each B-scan image was reconstructed with both the BP and the CS methods, and the Hilbert transform was taken after the reconstruction. After processing all the B-scans, the maximum amplitude projection (MAP) images were acquired through projecting the B-scans along the axial direction. Figures 4 and 4
show MAP images reconstructed with the BP and the CS methods, respectively, and one typical B-scan was extracted as shown in Figs. 4 and 4. We observed a significantly reduced noise level with the CS reconstruction. To further demonstrate the ability of the CS method in reducing the undersampling artifacts, we aggressively reconstruct the image with only 16 elements (
$1\u22153$
of total 48 elements) with both the BP and the CS method. The results are shown in Figs.
4, 4, 4, 4. With the BP reconstructions, the extremely sparse linear array generates significant grating lobe artifacts. By comparison, these undersampling artifacts were effectively reduced with the CS reconstruction.

## 4.

## Conclusions

Both the phantom and the *in vivo* results show that the CS method can effectively reduce the undersampling artifacts. By incorporating the CS theory in the PAT reconstruction, we can effectively reduce the system cost, or cover a larger FOV with the same number of measurements. Although the CS method is demonstrated here only with 2-D problems, the generalization to 3-D reconstructions is straightforward.

## Acknowledgements

This work was sponsored in part by National Institutes of Health Grant Nos. R01 EB008085, U54 CA136398, R01 CA113453901, R01 NS46214 (BRP), and R01 EB000712. L.W. has a financial interest in Microphotoacoustic, Inc., and Endra, Inc., which, however, did not support this work.

## References

**,” Rev. Sci. Instrum., 77 (4), 041011 (2006). 0034-6748 Google Scholar**

*Photoacoustic imaging in biomedicine***,” IEEE J. Sel. Top. Quantum Electron., 14 (1), 171 –179 (2008). https://doi.org/10.1109/JSTQE.2007.913398 1077-260X Google Scholar**

*Tutorial on photoacoustic microscopy and computed tomography***,” Med. Phys., 26 (9), 1832 –1837 (1999). https://doi.org/10.1118/1.598688 0094-2405 Google Scholar**

*Thermoacoustic computed tomography—technical considerations***,” Rev. Sci. Instrum., 70 (9), 3744 –3748 (1999). https://doi.org/10.1063/1.1149986 0034-6748 Google Scholar**

*Microwave-induced acoustic imaging of biological tissues***,” Nat. Biotechnol., 21 (7), 803 –806 (2003). https://doi.org/10.1038/nbt839 1087-0156 Google Scholar**

*Noninvasive laser-induced photoacoustic tomography for structural and functional**in vivo*imaging of the brain**,” Nat. Biotechnol., 24 (7), 848 –851 (2006). https://doi.org/10.1038/nbt1220 1087-0156 Google Scholar**

*Functional photoacoustic microscopy for high-resolution and noninvasive**in vivo*imaging**,” Radiology, 210 (2), 587 –587 (1999). 0033-8419 Google Scholar**

*Thermoacoustic computed tomography for breast imaging***,” Proc. SPIE, 4618 81 –94 (2002). https://doi.org/10.1117/12.469851 0277-786X Google Scholar**

*Optoacoustic imaging of blood for visualization and diagnostics of breast cancer***,” J. Biomed. Opt., 12 (2), 020504 (2007). https://doi.org/10.1117/1.2717531 1083-3668 Google Scholar**

*Photoacoustic imaging of lacZ gene expression**in vivo***,” Proc. IEEE, 96 (3), 481 –489 (2008). https://doi.org/10.1109/JPROC.2007.913515 0018-9219 Google Scholar**

*Simultaneous molecular and hypoxia imaging of brain tumors**in vivo*using spectroscopic photoacoustic tomography**,” Nat. Nanotechnol., 3 (9), 557 –562 (2008). https://doi.org/10.1038/nnano.2008.231 1748-3387 Google Scholar**

*Carbon nanotubes as photoacoustic molecular imaging agents in living mice***,” IEEE Trans. Med. Imaging, 21 (7), 823 –828 (2002). https://doi.org/10.1109/TMI.2002.801172 0278-0062 Google Scholar**

*Exact frequency-domain reconstruction for thermoacoustic tomography—I: Planar geometry***,” IEEE Trans. Med. Imaging, 21 (7), 829 –833 (2002). https://doi.org/10.1109/TMI.2002.801171 0278-0062 Google Scholar**

*Exact frequency-domain reconstruction for thermoacoustic tomography—II: Cylindrical geometry***,” IEEE Trans. Biomed. Eng., 50 (9), 1086 –1099 (2003). https://doi.org/10.1109/TBME.2003.816081 0018-9294 Google Scholar**

*Time-domain reconstruction algorithms and numerical simulations for thermoacoustic tomography in various geometries***,” Phys. Rev. E, 75 (5), 016706 (2005). https://doi.org/10.1103/PhysRevE.71.016706 1063-651X Google Scholar**

*Universal back-projection algorithm for photoacoustic computed tomography***,” SIAM J. Math. Anal., 35 (5), 1213 –1240 (2004). https://doi.org/10.1137/S0036141002417814 0036-1410 Google Scholar**

*Determining a function from its mean values over a family of spheres***,” Med. Phys., 22 (10), 1605 –1609 (1995). https://doi.org/10.1118/1.597429 0094-2405 Google Scholar**

*Photoacoustic ultrasound (Paus)—reconstruction tomography***,” Phys. Med. Biol., 46 (7), 1863 –1872 (2001). https://doi.org/10.1088/0031-9155/46/7/309 0031-9155 Google Scholar**

*Temporal backward projection of optoacoustic pressure transients using Fourier transform methods***,” Opt. Lett., 23 (8), 648 –650 (1998). https://doi.org/10.1364/OL.23.000648 0146-9592 Google Scholar**

*Three-dimensional photoacoustic imaging of blood vessels in tissue***,” Inverse Probl., 23 (6), S95 –S112 (2007). https://doi.org/10.1088/0266-5611/23/6/S08 0266-5611 Google Scholar**

*Photoacoustic tomography with a limited-aperture planar sensor and a reverberant cavity***,” Opt. Express, 17 (13), 10489 –10498 (2009). https://doi.org/10.1364/OE.17.010489 1094-4087 Google Scholar**

*A real-time photoacoustic tomography system for small animals***,” Opt. Express, 16 (26), 21570 –21581 (2008). https://doi.org/10.1364/OE.16.021570 1094-4087 Google Scholar**

*Four-dimensional photoacoustic imaging of moving targets***,” J. Biomed. Opt., 13 (5), 054052 (2008). https://doi.org/10.1117/1.2992131 1083-3668 Google Scholar**

*Three-dimensional photoacoustic imaging by sparsearray detection and iterative image reconstruction***,” Inverse Probl., 23 (6), S21 –S35 (2007). https://doi.org/10.1088/0266-5611/23/6/S03 0266-5611 Google Scholar**

*Application of inverse source concepts to photoacoustic tomography***,” IEEE Trans. Med. Imaging, 28 (11), 1727 –1735 (2009). https://doi.org/10.1109/TMI.2009.2022623 0278-0062 Google Scholar**

*A reconstruction algorithm for photoacoustic imaging based on the nonuniform FFT***,” IEEE Trans. Inf. Theory, 52 (2), 489 –509 (2006). https://doi.org/10.1109/TIT.2005.862083 0018-9448 Google Scholar**

*Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information***,” Magn. Reson. Med., 58 (6), 1182 –1195 (2007). https://doi.org/10.1002/mrm.21391 0740-3194 Google Scholar**

*Sparse MRI: the application of compressed sensing for rapid MR imaging***,” IEEE Trans. Med. Imaging, 28 (4), 585 –594 (2009). https://doi.org/10.1109/TMI.2008.2007825 0278-0062 Google Scholar**

*The Application of compressed sensing for photo-acoustic tomography***,” J. Acoust. Soc. Am., 112 (4), 1536 –1544 (2002). https://doi.org/10.1121/1.1501898 0001-4966 Google Scholar**

*Iterative reconstruction algorithm for optoacoustic imaging***,” Med. Phys., 36 (9), 4084 (2009). https://doi.org/10.1118/1.3187231 0094-2405 Google Scholar**

*On the speckle-free nature of photoacoustic tomography***,” IEEE J. Sel. Topics Signal Process., 1 (4), 586 –597 (2007). https://doi.org/10.1109/JSTSP.2007.910281 Google Scholar**

*Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems***,” Appl. Phys. Lett., 93 (3), 033902 (2008). https://doi.org/10.1063/1.2963365 0003-6951 Google Scholar**

*High-numerical-aperture-based virtual point detectors for photoacoustic tomography***,” J. Biomed. Opt., 13 (5), 054028 (2008). https://doi.org/10.1117/1.2976141 1083-3668 Google Scholar**

*Fast 3-D dark-field reflection-mode photoacoustic microscopy**in vivo*with a $30\text{-}\mathrm{MHz}$ ultrasound linear array**,” SIAM J. Sci. Comput. (USA), 31 (2), 890 –912 (2008). https://doi.org/10.1137/080714488 1064-8275 Google Scholar**

*Probing the Pareto frontier for basis pursuit solutions*## Appendices

### Appendix

In this section, we describe the reconstruction algorithm for solving the constrained optimization problem in Eq. 3, which has been proved to be closely related to solving the following convex unconstrained optimization problem:^{35}

## Eq. 4

$$\underset{x}{\mathrm{min}}\phantom{\rule{0.2em}{0ex}}f\left(x\right)=\frac{1}{2}{|\Phi x-y|}_{2}+\lambda {\left|\Psi x\right|}_{1},$$^{35}Therefore, we solved a series of $\lambda $ to find a suboptimal solution of the problem in Eq. 3. The process is described as follows:

1. Let $\lambda =0.05{\left|{A}^{T}y\right|}_{\infty}$ Ref. 32, and solve problem (4) for $x$ .

2. Check the condition ${|\Phi x-y|}_{2}<\epsilon $ . If this condition holds, we increase $\lambda $ to promote the sparsity ${\left|\Psi x\right|}_{1}$ ; otherwise, we decrease $\lambda $ to enforce the data consistency ${|\Phi x-y|}_{2}$ .

3. Problem (4) is solved again with the new $\lambda $ . Previously solved $x$ is used as the initial guess. By using this warm starting technique,

^{32}the current optimization process takes for fewer numbers of iterations than the previous one.4. Steps 2 and 3 are repeated multiple times.

^{28}which is computationally efficient to implement.

As mentioned in Sec. 2.1, the forward problem matrix $\Phi $ is extremely large, and the direct matrix operation is computationally impractical. Therefore, the computations of both $\Phi $ and its transpose ${\Phi}^{T}$ were implemented as submodules.