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 ( 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 recommended1, 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 methods16 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) theory27 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.
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 at the detecting aperture, PAT tries to reconstruct the image of the initial pressure rise distribution . The forward and inverse problems in PAT express the reciprocal relationship between and . By solving the wave equation, the forward problem, which predicts by , can be derived as (assuming a delta pulse heating):is the speed of ultrasound, and is the position of the ultrasonic sensor.2 Sometimes, the velocity potential is employed to simplify Eq. 1:
The analytical inversion of Eq. 1 describes the inverse problem, which reconstructs with :, is the detecting aperture, and is the solid-angle weighting factor.2
To numerically model the forward and inverse problems, we need to properly discretize Eqs. 1, 2. We use a vector to represent , where each element of is the average value of initial pressure per unit volume. The size of depends on the field of view (FOV) and the desired spatial resolution of the reconstructed image. We use a vector to denote the velocity potential measured by all elements of the sensor array as a function of time. The size of is the number of detecting positions times the number of temporal samples at each position . Then, the forward problem can be described as , where matrix is the projection matrix. Similarly, the inverse problem can be written as , where represents the inverse process, and is the reconstructed image.
and are usually extremely large matrices (each containing data points), even for 2-D reconstruction problems. For example, when reconstructing a image with measurements from 512 detecting positions, where each position has 1024 time points, both and contain data points ( 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 and . For each detecting element , the forward and inverse operations are performed using a matrix of the same size as . Each element in this matrix stores an index of the temporal sample of measurement , and this index indicates where the corresponding element of should be projected. As a result, the storage space for both the forward and inverse operations for all elements is reduced to ( for the preceding example). To fully take advantage of parallel computing capability, the responses of all the elements can be calculated simultaneously.
Compressed Sensing for PAT
If the measurement is incomplete, matrix is ill-conditioned, and is obviously not an exact inversion of . 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.
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 to represent an -pixel image and to denote the wavelet basis set, then can be expanded as , where is the coefficient sequence of . 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 nature31 of PAT images further reduces the number of significant transform coefficients.
Since the object can be visually losslessly reconstructed with only a few large transform coefficients, the problem of sensing is equivalent to capturing these large coefficients in the representation domain . The forward problem of PAT can be seen as projecting the object to the sensing basis set , and the measurements are the resulting coefficients. The CS theory requires the two basis sets and to be incoherent, i.e., the sensing waveforms should have a dense representation in . In other words, the undersampled sensing basis should induce only incoherent artifacts that spread out and appear as random noise in .
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 , and we assume that a circular detecting aperture is uniformly sampled by multiple ultrasonic point sensors. The ’th transform coefficient in the domain [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 and transformed back to the image space [Fig. 1] with the inverse operator . Last, the reconstructed image is again transformed to the sparse domain [Fig. 1] with the forward discrete wavelet transform (FDWT). TPSF can be mathematically described as , and it measures the leakage of energy away from the ’th coefficient to other coefficients. The CS theory requires us to properly choose and so that these interferences can be minimized and spread out in . Readers are referred to Reference 29 for a quantitative comparison of the TPSF maps for various in PAT.
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 is obtained by solving the following constrained optimization problem:and are defined as earlier, is the measured data, and is the parameter that controls the fidelity of the reconstruction to . The parameter is usually set based on the expected noise level. The object function in Eq. 3 is the norm (defined as ). The norm is used here instead of the norm [defined as ], because the norm penalizes large coefficients heavily, and leads to nonsparsity. In the norm, many small coefficients tend to carry a much larger penalty than a few large coefficients; therefore, small coefficients are suppressed and solutions are often sparse. In Eq. 3, minimizing the norm of promotes sparsity, and the constraint enforces data consistency.
The algorithm is implemented with a nonlinear conjugate gradient descent method,32 as detailed in Sec. 5. On a laptop with a dual-core CPU and memory, the calculations usually take less than using MATLAB 2008a.
Results and Discussions
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 . Laser pulses with a repetition rate of 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 .
Figure 2 shows the reconstruction results with the BP [Figs. 2, 2, 2, 2], the CS [Figs. 2, 2, 4, 4], and the traditional iterative30 (IR) [Figs. 2, 2, 2, 2] methods, with 240, 120, 80, and 60 tomographic angles. The images are reconstructed with an FOV of . 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.
In Vivo Experiment
The first in vivo experiment was based on a custom-designed 512-element photoacoustic tomography array system.22 The 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 . 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 broadband linear transducer array has a total of 48 elements of dimensions with 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 ( 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.
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.
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.
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:35is a nonnegative regularization parameter, which determines the trade-off between the data consistency and the sparsity. In order for these two problems to be equivalent, and must satisfy a special relationship. However, it is difficult to find analytical solutions if the matrix is not orthogonal.35 Therefore, we solved a series of to find a suboptimal solution of the problem in Eq. 3. The process is described as follows:
1. Let Ref. 32, and solve problem (4) for .
2. Check the condition . If this condition holds, we increase to promote the sparsity ; otherwise, we decrease to enforce the data consistency .
3. Problem (4) is solved again with the new . Previously solved 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.
As mentioned in Sec. 2.1, the forward problem matrix is extremely large, and the direct matrix operation is computationally impractical. Therefore, the computations of both and its transpose were implemented as submodules.