## Abstract

Compressed sensing (CS) can recover sparse signals from under-sampled measurements. In this work, we have developed an advanced CS framework for photoacoustic computed tomography (PACT). During the reconstruction, a small part of the nonzero signals’ locations in the transformed sparse domain is used as partially known support (PKS). PACT reconstructions have been performed with under-sampled *in vivo* image data of human hands and a rat. Compared to PACT with basic CS, PACT with CS-PKS can recover signals using fewer ultrasonic transducer elements and can improve convergence speed, which may ultimately enable high-speed, low-cost PACT for various biomedical applications.

©2012 Optical Society of America

## 1. Introduction

Photoacoustic imaging is a novel noninvasive imaging modality combining the advantages of both optical imaging and ultrasonic imaging [1]. Compared with ultrasonic imaging, it provides high optical contrast sensitive to important physiological parameters (such as the oxygen saturation of hemoglobin and the metabolic rate of oxygen); compared with purely optical imaging, it offers better spatial resolution for deep imaging. Although young, It has attracted much attention with its wide applications in various preclinical imaging settings [2,3]. Moreover, clinical photoacoustic imaging technologies are also emerging [4]. By now, it has been explored for the early diagnosis of cancer, the imaging of tumor angiogenesis, and identification of cardiovascular vulnerable plaques, in both preclinical and clinical studies [4–7]. Photoacoustic computed tomography (PACT) is one form of photoacoustic imaging [1], based on the reconstruction of the internal photoacoustic source distributions from the ultrasonic measurements over a plane or a volume. Generally, it uses either an unfocused ultrasonic transducer with mechanical scanning or an ultrasonic transducer array to acquire the image data [8–10]; then, the optical properties of the imaged objects can be recovered using a reconstruction algorithm.

For many biomedical applications, PACT with a high-frequency ultrasonic array is of particular interest, owing to its capability of both high-speed and high-resolution imaging [11]. However, the cost of a high-frequency ultrasonic array can be very high. In addition, due to challenges in fabrication, a fair amount of today’s high-frequency ultrasonic arrays are not dense enough to spatially fully sample the signals, which requires a pitch that is half of the center ultrasonic wavelength or smaller. For example, in one high-frequency linear-array PACT system, the pitch of the ultrasonic array is as large as twice of the center ultrasonic wavelength [12], resulting in reconstruction artifacts. Generally, to achieve high-quality PACT with conventional back-projection reconstruction, a number of densely packed ultrasonic transducers (or array elements) are needed, which can give rise to a high cost, and a long data acquisition and processing time.

Compressed sensing (CS) is a novel information theory capable of recovering signals that are sparse from under-sampled measurements. To date, CS has been widely studied and used in MRI [13], diffuse optical tomography [14] and CT [15]. Recently, explorations aiming to leverage the advantages of CS for PACT have also been reported. For example, J. Provost *et al.* have conducted numerical and phantom studies on CS-based PACT [16]; D. Liang *et al.* have numerically demonstrated CS-based photoacoustic imaging under random optical illumination [17]; Z. Guo *et al.* have reported CS-based photoacoustic reconstruction in time domain with images acquired *in vivo* [18]; M. Sun *et al.* have developed an arc-direction compressed-sensing PACT algorithm with numerical phantoms [19]. All these studies have shown that CS-based reconstruction techniques can recover photoacoustic signals sampled at a (temporal and/or spatial) rate lower than that required by the Nyquist sampling theory. Therefore, both the number of ultrasonic transducers required for the PACT system and the amount of acquired data needed for reconstruction can be significantly reduced.

However, although CS has been used for PACT, this exploration is in its infancy. To our knowledge, the use of CS-based PACT in frequency domain has not been studied and validated with *in vivo* image data yet. More importantly, many advanced CS theories that may of great potential benefits for PACT still remain unexplored. For instance, recently, N. Vaswani *et al.* developed a modified-CS method, in which partially known support (PKS), i.e., some known nonzero locations in the transformed sparse domain, was used as an additional prior besides sparsity [20]. Theoretical studies and applications in MRI have shown that PKS can further reduce the number of measurements needed for high-fidelity reconstruction [21], a significant benefit of great interest to PACT as well—including the potential to reduce the number of ultrasonic transducers (or array elements) and the image data acquisition time.

In this work, first, using *in vivo* photoacoustic image data acquired with a linear ultrasonic array, we have validated our basic CS reconstruction methods for PACT in both time and frequency domains. The CS-based PACT has produced significantly better images than those obtained with the conventional back-projection (BP) PACT. In addition, compared to CS-based PACT in time domain, its counterpart in frequency domain has produced images of a better contrast-to-noise ratio (CNR), with a significantly reduced computation time. Most importantly, for the first time to our knowledge, we have developed an advanced CS reconstruction model that incorporated PKS for PACT reconstruction. To solve the complex inverse problem, an iteratively reweighed conjugate gradient descent method (IR-CGD) has been developed. Results with *in vivo* image data have shown that, even in the case of heavily under-sampling, the developed CS-PKS strategy can recover the photoacoustic signals accurately, suppress the under-sampling artifacts effectively, and converge rapidly with a relatively small number of iterations, all of which are not attainable with previously developed basic CS approaches. These advantages of CS-PKS may ultimately enable high-speed, low-cost PACT to be constructed with less densely packed ultrasonic arrays or less number of ultrasonic transducers for various biomedical applications.

## 2. Theory

#### 2.1 Photoacoustic computed tomography

According to the photoacoustic generation theory, the acoustic pressure $p(r,t)$at position $r$and time *t* satisfies the wave equation [22]

*c*is the speed of sound. By solving the wave equation, the forward problem, which predicts $p(r,t)$ by ${p}_{0}(r)$, can be written as

The above time dependent wave equation can be transformed into the temporal-frequency domain. Denoting the Fourier transform of $p$ as $\overline{p}$, we have [23]

#### 2.2 Compressed sensing

One important prerequisite of compressed sensing is that either the signal or its transform in a certain domain is sparse or compressible. A signal $x\in {R}^{N}$ is called *s*-sparse when *s* of its coefficients are nonzero (usually *s*<<*N*). A signal is *s*-compressible if its ordered set of coefficients decays rapidly and $x$ can be well approximated by the first *s* large coefficients. Fortunately, most medical images are sparse in a certain domain by finding an appropriate sparse transform $\psi :x=\psi \theta $, where $\theta $ is the original image, and $x$ is the transformed one. It has been proven that photoacoustic images can be transformed into sparse domains by a variety of transforms, such as the numerical derivative (ND), the wavelet transform, etc [16–19].

In PACT, the goal of reconstruction is to recover the photoacoustic image signal $\theta $ through the measurement data $y$ from ultrasonic transducers. Given that $y$ is obtained with a measurement matrix $K$, we have $y=K\theta $. Photoacoustic images with CS can usually be reconstructed by solving the following constrained optimization problem for $x$:

According to Eq. (2) and the back projection principle, the measurement matrix in time domain can be designed as [18]

From Eq. (4), the measurement matrix in temporal-frequency domain can be written as [16]

In Eq. (6) and Eq. (7), ${r}_{i,j}$ represents the Cartesian coordinates of image pixels, ${r}_{m}$ represents the positions of transducers, $p$ represents the number of transducers, ${q}_{s}$ and ${q}_{n}$ represent the number of sampling points in time and frequency domains, respectively.

#### 2.3 Compressed sensing with partially known support

Recently, a new type of CS reconstruction framework—CS with partially known support (CS-PKS)—has been developed to incorporate certain prior support of a sparse signal into the CS reconstruction process, where the support is defined as the locations of non-zero signals in the sparse domain [20]. Compared with basic CS, CS-PKS has been demonstrated in MRI to recover signals with fewer measurements or obtain higher quality images with the same number of measurements [21,24].

Let $x\in {R}^{N}$ be the sparse coefficients of the desired image and denote the support $T=supp(x)={T}_{0}\cup \Delta $, where ${T}_{0}\subset \left\{1,2,\cdots ,N\right\}$is the known part of $T$with a size of $\left|{T}_{0}\right|$ and $\Delta \subset \left\{1,2,\cdots ,N\right\}$ is the unknown part of $T$with a size of $\left|\Delta \right|$. With the partially known support ${T}_{0}$, the candidates for the *s*-sparse signal set $x$ are restricted in a signal space smaller than that in basic CS, and thus requires fewer samples to yield an accurate reconstruction [25]. Since the signal $x$ is known to be nonzero at some locations, CS-PKS allows us to minimize the number of nonzeros out of the known part during searching for a sparse solution to $\Phi x=y$. Using convex relaxation, this procedure can be formulated as [21]

## 3. Methods

#### 3.1 Reconstruction model

Based on the above theory of CS and PKS, the reconstruction models of basic CS and CS-PKS were developed for PACT, which are expressed by the Eq. (9) and Eq. (10), respectively.

*F*is the objective function consisting of three parts—the first part represents the square error between the estimated measurements from the reconstructed signal and the experimentally acquired measurements, the second part represents the ${l}_{1}$ norm of signal in a sparse domain, and the third part represents the total variation (TV) penalty of the signal, and $\alpha $ and $\beta $ are the regularization parameters determining the trade-off between the data consistency and the sparsity. The difference between the CS reconstruction model and the CS-PKS model is that the PKS prior information is incorporated into the second part of

*F*in Eq. (10). Please note that it is often useful to include a TV item in CS even when other sparse transforms are used in the objective function [13,17,26]. The definition of ${x}_{\Delta}$ in Eq. (10) is given by

To solve Eq. (10), an iterative reweighed method is employed, and Eq. (10) is rewritten as

**W**is a diagonal matrix whose diagonal element is set to 1 when the corresponding signal element belongs to the unknown support $\Delta $ in the ${(i-1)}^{th}$ iteration, and set to 0 otherwise. Note that

**W**is determined by the partial support ${T}_{0}$, which can be obtained by thresholding the reconstructed sparse coefficients [21]where ${x}_{z}^{(i)}$is the ${z}^{th}$element of reconstructed signal $x$ in the ${i}^{th}$ iteration and ${\tau}^{(i)}$ is a threshold used to determine the set of partially known support ${T}_{0}^{(i)}$. It has been found that, the iterative reconstruction process can work well for fast-decaying signal ${x}^{(i)}$by setting ${\tau}^{(i)}={\Vert {x}^{(i)}\Vert}_{\infty}/{\delta}^{(i)}$, where $\delta $ is a hyper-parameter.

#### 3.2 Optimization algorithm

In order to solve Eq. (12), an iteratively reweighed conjugate gradient descent method (IR-CGD) was developed. The result from the first few basic CS iterations, denoted as ${x}^{(0)}$, is used as the initial value for subsequent CS-PKS reconstruction. The major steps of the algorithm are specified below.

IR-CGD algorithm for CS-PKS based photoacoustic computed tomography |

Inputs: y, K, $\Psi $, $\alpha $, $\beta $. |

Step 1: Initialize ${\tau}^{(0)}$, ${T}_{0}^{(0)}$, ${W}^{(0)}$, and ${\Delta}^{(0)}$ from ${x}^{(0)}$ based on Eq. (11) and Eq. (13); set the outer loop variable $i=1$ and the total outer loop number = $I$. |

Step 2: Do the inner loop |

2.1 Initialize the inner loop variable $l=1$, set ${x}_{l}^{(i)}={x}^{(i-1)}$(${x}_{l}^{(i)}$represents the reconstructed signal of the ${l}^{th}$ inner loop in the ${i}^{th}$ outer loop). |

2.2 Compute the gradient $\nabla F$of the objective function $F$ with respect to signal ${x}_{l}^{(i)}$ in Eq. (11). |

2.3 Obtain ${x}_{l+1}^{(i)}$ with the CGD method. |

2.4 If ${\Vert {x}_{l+1}^{(i)}-{x}_{l}^{(i)}\Vert}_{2}<\epsilon $, go to step3; otherwise, let $l=l+1$, and go to step 2.2. |

Step 3: Set $i=i+1$, ${x}^{(i-1)}={x}_{l+1}^{(i-1)}$, update ${\tau}^{(i-1)}$, ${T}_{0}^{(i-1)}$, ${W}^{(i-1)}$, and ${\Delta}^{(i-1)}$. |

Step 4: If $i>I$, finish; otherwise, go to step2. |

#### 3.3 In vivo imaging

Noninvasive photoacoustic imaging of a rat and human hands was performed *in vivo* using a linear-array PACT system. The major components of the system include: (1) a tunable dye laser (Cobra, Sirah Laser-und Plasmatechnik GmbH, Germany) pumped by a Q-switched Nd:YLF laser (INNOSLAB, Edgewave GmbH, Germany), (2) a custom-built linear ultrasonic array with a center frequency of 30 MHZ, and (3) an 8-core PC (Dell Precision 490) equipped with an 8-channel high-speed PCI data acquisition (DAQ) card (Octopus CompuScope 8389, GaGe, USA). Specifically, the high-frequency ultrasonic array was fabricated from a 2-2-piezo-composite by the NIH Resource Center for Medical Ultrasonic Transducer Technology at the University of Southern California. The array had 48 elements (82 $\mu m$ × 2 mm) with 100 $\mu m$ spacing. The elements were elevationally focused at 8.2 mm, providing an elevational resolution of 200 – 300 $\mu m$within a ~3.5-mm focal zone. The mean fractional bandwidth was ~70% for receiving-only operation, as used in our PACT system. The axial and lateral resolutions of the system at 8 mm normal depth from the transducer surface were ~25 and ~80 $\mu m$, respectively. Further details of the imaging system can be found in our previous publications [10,12]. Note that, (1) the name “”microscopy” was used in references [10,12] only because high microscopic spatial resolution was achieved, the imaging system (which was also used for this study) is essentially a linear-detection geometry PACT; (2) because a linear array was used, our imaging system is considered as a limited-view PACT; (3) as the pitch of the ultrasonic array was 100 $\mu m$, about twice of the center ultrasonic wavelength, the array is considered as not spatially fully sampled [27].

For *in vivo* imaging experiments, the dye laser output was tuned to 584 nm—an isosbestic point at which the oxy- and deoxy-hemoglobin absorb equally. For small animal imaging, a Sprague Dawley rat (Harlan Laboratories, Inc., USA) was used. Intra-dermal injection of a mixture of ketamine (85 mg/kg) and xylazine (15 mg/kg) was implemented for anesthesia. Before photoacoustic imaging, the hair in the imaged region was removed with commercial depilatory lotion. For human hand imaging, the optical fluence on the skin surface was set to ~0.5 mJ/cm^{2} per pulse, well below the ANSI [28] recommended Maximum Permissible Exposure (MPE) of 20 mJ/cm^{2} for a single pulse in the visible spectral range. The time-averaged light intensity during image acquisition was ~150 mW/cm^{2}, also below the ANSI recommended MPE of 196 mW/cm^{2}, calculated by 1.1*t*^{0.25} W/cm^{2}. As the ANSI safety limit for this pulse width region is dominantly based on the thermal mechanism, our compliance with the ANSI standards guarantees no thermal damage to the tissue. All animal and human experiments described here were carried out in compliance with Washington University approved protocols.

## 4. Results

#### 4.1 Results in time and frequency domains

In this section, Eq. (9) was adopted as the basic CS reconstruction model, and a four-level Daubechies wavelet was used as the sparse transform. A photograph of the imaged hand (human hand-1) is shown in Fig. 1(A) , in which a dashed box indicates the imaged area (of ~6.4 mm × 9.6 mm). For each 3D image, 160 B-scan frames were acquired, with each reconstructed B-scan image consisting of 128 × 128 pixels (corresponding to a cross-section of ~6.4 mm × 1.6 mm). Note that, in all B-scans, although Hilbert transform was used for envelop-detection to display the images, it was applied only after reconstructions. The rest of Fig. 1 shows the reconstructed images with data from all 48 ultrasonic transducer elements. Figure 1(B) is the maximum amplitude projection image (MAP) —the maximum photoacoustic amplitudes projected along the depth direction to the surface of the hand—reconstructed by the conventional photoacoustic BP method. Figure 1(C) is the MAP image reconstructed by CS in time domain (TD), and Fig. 1(D) is the MAP image reconstructed by CS in frequency domain (FD), in which a band-pass filter of 0 – 40 MHz is applied. Figures 1(a) – (f) are the B-scan images along the positions represented by the vertical dash lines in Figs. 1(B) – (D). With data from all 48 transducer elements, it can be seen that, although the linear-array is not considered as spatially fully sampled, the photoacoustic images can be reconstructed well using any of the three reconstruction methods (note that the imaging system using BP for reconstruction—with data from all 48 transducer elements—has already been validated with optical microscopy of excised tissue in small animals [12]). In addition, compared with the TD-CS reconstruction, the FD-CS reconstruction has produced images of a slightly higher contrast-to-noise ratio (CNR), yet has consumed ~50% less computation time (for each B-scan reconstruction, it takes about 40 s using FD-CS vs. 75 s using TD-CS on a PC with a Intel Core i5-2400 CPU of 3.1 GHz).

Figure 2 shows the reconstructed images in the case of more heavily undersampling—with data from only 16 uniformly distributed transducer elements of the array. Similar to the approaches above, while all 128 time points per transducer element are used for BP and TD-CS reconstructions, only signals representing frequencies between 0 – 40 MHz (corresponding to 40 points in frequency domain) are used for the FD-CS reconstruction. In this case, the images reconstructed by BP (Figs. 2(A), (a) and (b)) are becoming significantly worse, while the images reconstructed by TD-CS and FD-CS (Figs. 2(B) – (C) and (c) – (f)) still retain the major features seen in the images reconstructed with data from all 48 transducer elements (Fig. 1). Similarly, compared with TD-CS, the FD-CS approach has shown a higher CNR with ~50% less reconstruction time.

The results suggest that, in the case of under-sampling, the FD-CS reconstruction can produce images of a higher CNR, compared to those obtained with BP and TD-CS. However, when using data from only 16 transducer elements (i.e., heavily undersampling), significant under-sampling artifacts emerge in all reconstructed images (Fig. 2).

#### 4.2 Results with CS-PKS method

To further improve the image quality, PKS was incorporated into our FD-CS photoacoustic reconstruction process (section 3). Specifically, the result after a few iterations of traditional CS was used as the initial image. To start the iterative reconstruction process, it was transformed into a sparse domain using a four-level Daubechies wavelet. According to Eq. (13), in each iteration, the locations (in the sparse domain) of which the coefficients are larger than a threshold were used as the partially known support. *In vivo* imaging of a human hand (human hand-2) and a rat were performed to validate and demonstrate the advantages of the developed CS-PKS strategy. All CS-PKS reconstructions were performed based on the optimization model defined by Eq. (12). The imaged volume is about 6.4 mm × 10 mm × 1.6 mm (in *x*, *y*, and *z* directions, respectively). For each of the following experiments, 166 B-scan frames were acquired, with each reconstructed B-scan image consisting of 128 × 128 pixels.

Figure 3 shows the reconstructed photoacoustic images of human hand-2. Figures 3(A) and (B) are the MAP images reconstructed by BP with data from 48 and 16 transducer elements, respectively. Figures 3(C) and (D) are the MAP images reconstructed by the basic FD-CS and the FD-CS-PKS methods, respectively, with data from only 16 transducer elements. Figures 3(a) – (h) are the corresponding B-scan images along the vertical dashed lines shown in Figs. 3(A) – (D). Compared with the reconstructed images using data from all 48 transducer elements (Figs. 3(A), (a), and (b)), the images reconstructed by BP with data from only 16 transducer elements are much worse (Figs. 3(B), (c) and (d))—the features are almost submerged in noise. With the use of regular CS, the image quality improves significantly (Figs. 3(C)), (e) and (f)), although some under-sampling artifacts are still notable. Surprisingly, when the PKS prior information was exploited, the under-sampling artifacts are greatly suppressed (Figs. 3(D), (g) and (h)), and the image can be recovered almost to a degree obtainable only with data from all 48 transducer elements using the conventional BP reconstruction. In addition, our experiments have also shown that the CS-PKS method is able to improve the convergence speed. For each B-scan reconstruction, it takes 7 s using FD-CS-PKS vs. 12 s using FD-CS, on a PC with an Intel Core i5-2400 CPU of 3.1 GHz.

*In vivo* imaging of a rat was also performed to further demonstrate the advantages of the developed FD-CS-PKS strategy. Figure 4
shows the reconstructed images of the subcutaneous vasculature of the back of the rat. Figures 4(A) and (B) are the MAP images reconstructed using BP with data from 48 and 16 transducer elements, respectively. Figures 4(C) and (D) are the MAP images reconstructed using the basic CS and the FD-CS-PKS methods, respectively, with data from 16 transducer elements. Figures 4(a) – (h) are the corresponding B-scan images along the dashed lines shown in Figs. 4(A) – (D). In this case, similar to the results obtained with human hand-2, BP recovers only a small part of useful features. Furthermore, while FD-CS clearly improves the image quality, FD-CS-PKS produces the best images of a high CNR and significantly reduced artifacts.

#### 4.3 Quantitative comparisons

From the above results, the advantages of our CS-PKS reconstruction strategy can be seen visually. In this section, to further compare the results, three quantitative analyses are performed (in which all images were normalized by dividing their corresponding maximum values).

The first quantitative approach for comparison is to compute the histograms of the difference images—images obtained by subtracting the control from the reconstructed images and then taking the absolute value (Fig. 5 ). In our calculation, the FD-CS image reconstructed with data from all 48 transducer elements was used as the control (using the image reconstructed by BP with data from all 48 transducer elements as the control leads to the same conclusions). Figures 5(A) – (C) are the difference histograms between the control and the images reconstructed by BP, TD-CS, and FD-CS, respectively, with data from 16 transducer elements. From Fig. 5(A) to Fig. 5(C), the number of counts lying in the lower ranges of the difference amplitudes increases significantly, suggesting that the images reconstructed using the CS-based methods (in particular, FD-CS) possessing higher fidelity. From Figs. 5(D) – (F), it is obvious that, for the images reconstructed with CS, the majority of the difference amplitudes lies in the lower ranges. In particular, in the case of FD-CS-PKS, the number of counts in [0, 0.1] represents >90% of all pixel numbers. Figures 5(G) – (I), obtained with the rat model, show very similar results.

Further, localized comparisons were performed by plotting the photoacoustic amplitudes along chosen lines in the MAP images (Fig. 6
). Figures 6(A) and (B) are the plots for human hand-2 and the rat along the pink lines shown in Fig. 3(A) and Fig. 4(A), respectively. In order to have a better quantitative comparison, the CNRs of selected signal peaks were also computed (Fig. 6, insets) (in our computation, the range from 0.0 mm to 1.5 mm along the *x*-axis is used as an estimation of the background). It can be seen that, (1) images from CS-based reconstructions generally have higher CNRs compared to those by the BP method, and (2) the images reconstructed by CS-PKS with only 16 transducers have shown equivalent CNRs to that of BP images with data from all 48 ultrasonic elements.

The third quantitative approach for comparison is based on the computation of the cross-correlation coefficient (*corr*), as defined by Eq. (14):

*i*and

*j*represent the row and column of the images. ${\overline{I}}^{r}$and ${\overline{I}}^{o}$are the mean values of the reconstructed images and the control, respectively. For all models (“Hand-1”, “Hand-2”, and “Rat”), using the images reconstructed by FD-CS48 as the control, the computed cross-correlation coefficients are shown in Table 1 . It can be seen clearly that the images reconstructed using CS have larger cross-correlation coefficients than those reconstructed by BP. Furthermore, the cross-correlation coefficients of the images reconstructed using our CS-PKS strategy are among the largest of all approaches.

## 5. Discussion and conclusions

The above results have clearly demonstrated the superiority of our advanced CS-PKS strategy over the basic CS approaches. Here, we discuss a few steps during the reconstruction process that are critical for achieving optimal results. First, in the optimization process with the IR-CGD method, two regularization parameters $\alpha $ and $\beta $—representing the weight of PKS and the ND regularization item, respectively—need to be determined appropriately. Overweighed values will result in the distortion of the reconstruction. On the contrary, if their proportions to the objective function are too small, their support would become ineffective. What’s more, the weight distribution between PKS and ND, i.e., the relative magnitude of $\alpha $ and $\beta $, also matters. Unfortunately, at present, there is no general method to effectively compute the optimal values of these regularization parameters, which vary with different reconstruction algorithms and raw data sets. In our experiments, they are empirically determined by trying different combinations and choosing the best one. Second, the parameter $\delta $, which is used to determine the set of the detection support, is another key factor affecting the final reconstruction results. According to the relationship between $\delta $ and $\tau $(section 3.1), a large $\delta $ will result in too many false locations in the detection support, and a small $\delta $ will result in too few locations in the detection support. Finally, several other factors may also influence the reconstructed image quality, including the step for updating signals in each iteration, the sparsity of the signals, and the choice of the sparse transform basis. The best values of the reconstruction parameters used in our experiments are listed in Table 2 .

In summary, we have developed CS-based reconstruction methods for PACT in time and frequency domains. Supported by *in vivo* image data, we have found that, in the case of under-sampling, the CS-based photoacoustic reconstructions are superior to the conventional photoacoustic BP reconstruction. Furthermore, the CS-PKS reconstruction framework developed by us has demonstrated accurate reconstruction results even with heavily under-sampled data from only 16 ultrasonic transducer elements. Compared to PACT with basic CS, PACT with the developed CS-PKS strategy has demonstrated significant advantages, including a higher CNR, a faster convergence, and a more effective suppression of under-sampling artifacts. The promising results demonstrated in this work suggest that the development of novel advanced CS reconstruction strategies has the potential to enable high-quality PACT with a significantly reduced number of ultrasonic transducers (or array elements). In return, this will reduce the cost and the amount of image data required for accurate reconstruction, which may ultimately have a profound impact in promoting various preclinical and clinical applications of PACT.

## Acknowledgments

This work was supported in part by the Shenzhen Development and Reform Commission Fund [2010] 1599 (Key Technologies for Medical Imaging), the Guangdong Innovation Research Team Fund for Low-cost Healthcare Technologies (GIRTF-LCHT), the Low-cost Health Engineering Research Program of the Chinese Academy of Sciences [2011] to L.S., as well as the National Natural Science Foundation of China 61102042, 81120108012, and the Basic Research Program of Shenzhen JC201104220219A to D.L.. J. Meng would like to thank H. F. Shi and J. H. Chen for beneficial discussions on the analysis of frequency signals.

## References and links

**1. **L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science **335**(6075), 1458–1462 (2012). [CrossRef] [PubMed]

**2. **G. Ku, X. D. Wang, X. Y. Xie, G. Stoica, and L. V. Wang, “Imaging of tumor angiogenesis in rat brains in vivo by photoacoustic tomography,” Appl. Opt. **44**(5), 770–775 (2005). [CrossRef] [PubMed]

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

**4. **S. A. Ermilov, T. Khamapirad, A. Conjusteau, M. H. Leonard, R. Lacewell, K. Mehta, T. Miller, and A. A. Oraevsky, “Laser optoacoustic imaging system for detection of breast cancer,” J. Biomed. Opt. **14**(2), 024007 (2009). [CrossRef] [PubMed]

**5. **C. Kim, C. Favazza, and L. V. Wang, “In vivo photoacoustic tomography of chemicals: high-resolution functional and molecular optical imaging at new depths,” Chem. Rev. **110**(5), 2756–2782 (2010). [CrossRef] [PubMed]

**6. **Z. X. Xie, W. Roberts, P. Carson, X. Liu, C. Tao, and X. Wang, “Evaluation of bladder microvasculature with high-resolution photoacoustic imaging,” Opt. Lett. **36**(24), 4815–4817 (2011). [CrossRef] [PubMed]

**7. **W. Wei, X. Li, Q. F. Zhou, K. K. Shung, and Z. Chen, “Integrated ultrasound and photoacoustic probe for co-registered intravascular imaging,” J. Biomed. Opt. **16**(10), 106001 (2011). [CrossRef] [PubMed]

**8. **X. D. Wang, Y. J. Pang, G. Ku, X. Y. Xie, G. Stoica, and L. V. Wang, “Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain,” Nat. Biotechnol. **21**(7), 803–806 (2003). [CrossRef] [PubMed]

**9. **C. H. Li, A. Aguirre, J. Gamelin, A. Maurudis, Q. Zhu, and L. V. Wang, “Real-time photoacoustic tomography of cortical hemodynamics in small animals,” J. Biomed. Opt. **15**(1), 010509 (2010). [CrossRef] [PubMed]

**10. **L. Song, K. Maslov, K. K. Shung, and L. V. Wang, “Ultrasound-array-based real-time photoacoustic microscopy of human pulsatile dynamics in vivo,” J. Biomed. Opt. **15**(2), 021303 (2010). [CrossRef] [PubMed]

**11. **L. Song, K. Maslov, and L. V. Wang, “Multifocal optical-resolution photoacoustic microscopy *in vivo*,” Opt. Lett. **36**(7), 1236–1238 (2011). [CrossRef] [PubMed]

**12. **L. Song, K. Maslov, R. Bitton, K. K. Shung, and L. V. Wang, “Fast 3-D dark-field reflection-mode photoacoustic microscopy in vivo with a 30-MHz ultrasound linear array,” J. Biomed. Opt. **13**(5), 054028 (2008). [CrossRef] [PubMed]

**13. **M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The application of compressed sensing for rapid MR imaging,” Magn. Reson. Med. **58**(6), 1182–1195 (2007). [CrossRef] [PubMed]

**14. **D. Han, J. Tian, K. Liu, J. C. Feng, B. Zhang, X. Ma, and C. H. Qin, “Sparsity-promoting tomographic fluorescence imaging with simplified spherical harmonics approximation,” IEEE Trans. Biomed. Eng. **57**(10), 2564–2567 (2010). [CrossRef] [PubMed]

**15. **G. H. Chen, J. Tang, and S. Leng, “Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets,” Med. Phys. **35**(2), 660–663 (2008). [CrossRef] [PubMed]

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

**17. **D. Liang, H. F. Zhang, and L. Ying, “Compressed-sensing photoacoustic imaging based on random optical illumination,” IJFIPM **2**(4), 394–406 (2009). [CrossRef]

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

**19. **M. J. Sun, N. Z. Feng, Y. Shen, X. L. Shen, L. Y. Ma, J. G. Li, and Z. H. Wu, “Photoacoustic imaging method based on arc-direction compressed sensing and multi-angle observation,” Opt. Express **19**(16), 14801–14806 (2011). [CrossRef] [PubMed]

**20. **N. Vaswani and W. Lu, “Modified-CS: Modifying compressive sensing for problems with partially known support,” IEEE T Signal Process. **58**, 4595–4607 (2010).

**21. **D. Liang, E. V. R. DiBella, R. R. Chen, and L. Ying, “k-t ISD: Dynamic cardiac MR imaging using compressed sensing with iterative support detection,” Magn. Reson. Med. , doi:. [CrossRef]

**22. **M. H. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **71**(1), 016706 (2005). [CrossRef] [PubMed]

**23. **Y. Xu, D. Z. Feng, and L. V. Wang, “Exact frequency-domain reconstruction for thermoacoustic tomography--I: Planar Geometry,” IEEE Trans. Med. Imaging **21**(7), 823–828 (2002). [CrossRef] [PubMed]

**24. **Y. Wang and W. Yin, “Sparse signal reconstruction via iterative support detection,” SIAM J. Imaging Sciences **3**(3), 462–491 (2010). [CrossRef]

**25. **L. Jacques, “A short note on compressed sensing with partially known signal support,” Signal Process. **90**(12), 3308–3312 (2010). [CrossRef]

**26. **Y. Tsaig and D. Donoho, “Extensions of compressed sensing,” Signal Process. **86**, 549–571 (2006). [CrossRef]

**27. **K. K. Shung and M. Zippuro, “Ultrasonic transducers and arrays,” IEEE Eng. Med. Biol. **15**(6), 20–30 (1996). [CrossRef]

**28. **The Laser Institute of America, American National Standard for Safe Use of Lasers, (ANSI Z136.1–2000), The Laser Institute of America (2000).