High-resolution imaging with optical contrast in deep scattering tissue has been a sought-after goal for a long time, challenged by multiple scattering. Photoacoustic imaging is a rapidly developing modality providing ultrasonic acoustic resolution in deep tissue with optical absorption contrast [1,2]. A standard beamforming technique used to form photoacoustic images is delay-and-sum beamforming (DSB). DSB is effective in a broad range of contexts, but does not naturally support the incorporation of a priori information about the imaging target. In this Letter, we demonstrate photoacoustic imaging resolution improvement at depth by incorporating a priori spatial sparsity in the beamforming process, referring to the approach as “sparsity-based reconstruction (SBR).”

Some recent work has sought to improve photoacoustic resolution, such as in the context of optically resolved photoacoustic microscopy (OR-PAM) [3–6]. OR-PAM can achieve high resolution by using a scanned focused beam but, consequently, is limited by scattering to imaging depths less than the transport length. However, since we focus on improving resolution at greater depths than the transport length, we restrict our attention to the photoacoustic computed tomography (PACT) setting. PACT uses an unfocused excitation and reconstructs an image using a beamforming algorithm on data recorded from an array of transducer elements. Consequently, PACT can image deeper than OR-PAM, but it is limited to acoustic resolution due to the presence of significant scattering.

It should be noted that some work has focused on using spatial sparsity to improve resolution in an ultrasound imaging context [7,8]. We aim to apply similar ideas to a PACT context. In this context, some recent work has sought to incorporate a priori sparsity, but primarily in non-spatial domains or for the purpose of reducing the number of measurements [9–13]. The work in Ref. [14] uses spatial sparsity to improve PACT resolution on post-beamformed data (with mixed results), but primarily focuses on using randomly structured illumination together with spatial sparsity to increase resolution. In addition, Ref. [14] limits testing to the imaging of point-like beads and uses DSB output as input for the reconstruction algorithm. We extend this work by imaging wires (allowing for resolution characterization), and by using raw channel data (avoiding information loss) as input to the reconstruction algorithm. Beyond these conceptual differences, we use a high-frequency (21 MHz) transducer, which is important for high-resolution imaging in small animals. This is significant because our long-term objective is to image to centimeter-scale depths with cellular-scale resolution, which is important for many small animal imaging studies.

The key idea of SBR is to incorporate known spatial sparsity during beamforming. To do this, we first write a linear model for the channel data generated by photoacoustic excitation, illustrated in Fig. 1:

 

Fig. 1. Received pressure-over-time channel data are modeled as the linear superpositions of the responses due to individual point targets. Estimating $f$ produces an image and can be done using SBR or DSB.

Download Full Size | PDF

Here $g$ is the received pressure-over-time channel data generated by all the optical absorbers when a single laser pulse is fired. The $j_{\mathrm{th}}$ column $H_j$ of the “dictionary” matrix $H$ is the receive channel data obtained from a unit strength point absorber located at location $r_j$, reshaped into a column vector. We define $f$ as a column vector of absorber strengths so its $j_{\mathrm{th}}$ component $f_j$ is the strength of the source at location $r_j$. Finally, $n$ is a noise vector. The dictionary matrix relates input information (absorber locations) to output information (photoacoustic channel data). Thus, we model the receive channel data due to a collection of optically absorbing point targets as a linear combination of the receive channel data due to each source individually, weighted by source strength.

DSB can be used to solve Eq. (1), but it does not incorporate sparsity information. This is because DSB simply consists of summation across channel data aligned in time so that pressure data that originated from the receive focus point of interest interfere constructively, but so that other information interferes incoherently. To incorporate sparsity information during the image formation process, we use the following optimization program [15]:

Here the ${\ell }_1$-norm term (${\Vert f\Vert }_1=\sum _i|f_i|$) encourages a sparse solution, and the ${\ell }_2$-norm term (${\Vert f\Vert }_2=(\sum _i{|f_i|}^2{)}^{1/2}$) encourages the selection of ultrasound sources $\widehat{f}$ with an expected response similar to the measured response, so $g\approx H\widehat{f}$. The weighting parameter $\tau \in \mathbb{R}$ sets the relative importance of each term. Reshaping $\widehat{f}$ results in an image optimized for sparsity. This optimization program forms an image while incorporating sparsity and working with pre-beamformed raw channel data. For an investigation of the theoretical performance of this beamforming system in an ocean acoustic data context; see Ref. [16].

To solve Eq. (2), we used the L1-homotopy solver [17], as we found it to quickly generate reasonable estimates $\widehat{f}$. Its out-of-the-box ability to make effective use of multiple cores was especially helpful. To use the solver, we had to provide the inputs $\tau$, $H$, and $g$; and now we discuss how these were generated.

To set $\tau$, we manually selected a value that simultaneously generated qualitatively acceptable output across a range of cross-sectional sub-experiments. This took a handful of iterations.

To construct the dictionary elements of the $H$ matrix, we first tried to generate expected point responses using a simulation package. However, we found that experimental data worked better in practice. To generate dictionary elements experimentally, we extracted channel data from a single wire at a fixed location. The SBR reconstruction (using an initial, imprecise dictionary) for this single wire yielded a small cluster of points (when the sparsity weighting was de-emphasized), rather than an individual point. We then used the location of these cluster points to generate a superposition of channel data point responses to represent a single experimental dictionary element for a fixed wire location. Applying this technique in a two-wire context allowed us to extract an experimental estimate for the response due to a single point in a two-wire setting. Dictionary elements for other wire locations were then generated using appropriate delays and amplitude scaling of the channel data.

We used a dictionary $H$ with pressure-over-time responses from a $60\times 60$ grid, balancing the need for precision in the formed image against the computational time required. The grid stretched from $-0.25$ to 0.35 mm laterally and stretched from 8.5 to 8.75 mm axially (where the transducer surface was defined to be at 0 mm axially).

We generated the channel data $g$ using the experimental setup we now describe. We performed photoacoustic imaging using a 10 Hz repetition rate Nd:YAG laser (Surelite OPO Plus, Continuum) and a high-frequency ultrasound transducer (LZ-250, FUJIFILM Visualsonics) controlled by a programmable ultrasound system (Vantage 256, Verasonics, US). The experimental configuration is illustrated in Fig. 2. The ultrasound transducer used had a center frequency of 21 MHz, a bandwidth spanning 13–24 MHz, 256 channels, a kerf of 18 μm, and an element width of 72 μm. Channel data were sampled at 62.5 MHz and upsampled to 200 MHz, and the channel data vector $g$ was formed as an average across 16 measurements to reduce noise. The excitation laser had a wavelength of 710 nm, and imaging was conducted through $\approx 8.5\mathrm{mm}$ of 1% intralipid solution, a common choice for a tissue mimicking medium. The imaging target was chosen to test the resolving ability of the imaging technique; therefore, we imaged successive cross sections of two non-parallel wires. We used a programmable stage to achieve a 58.6 μm step size between cross sections, striking a compromise between image fineness and imaging time. The wires used were 17.8 μm in diameter and made of aluminum (ALW-29S, Heraeus).

 

Fig. 2. Experimental apparatus to measure the photoacoustic response of the cross section of two converging wires.

Download Full Size | PDF

To conduct the optimization routine on the experimental data, both $g$ and the columns of $H$ were first normalized to the same maximum. In addition, we used a sparsity weighting value of $\tau =115$ for all transverse slices. This value of $\tau$ was selected manually after a few iterations, made possible by the fact that the results were not very sensitive to $\tau$. However, it should be noted that as $\tau$ is increased, both low-intensity artifacts and low-intensity portions of the desired image tend to disappear.

We also performed DSB on the collected channel data, allowing for a comparison with SBR. DSB was performed using unity apodization, and the receive beamforming aperture size was selected to set $f_{\#}\approx 1$ at the depth of the target.

The raw experimental SBR results are shown in Fig. 3. The C-scan images were formed by taking maximum amplitude projections of B-scans. Each B-scan shows the intensity estimated for each element in the $60\times 60$ dictionary. We also show the results of applying a Gaussian blurring operation to the raw image and upsampling by a factor of 10 in Fig. 4. Finally, the 10-times upsampled DSB image is shown in Fig. 5.

 

Fig. 3. (a) and (b) are experimental B-scans formed using SBR; (c) is a maximum amplitude projection C-scan formed using SBR.

Download Full Size | PDF

 

Fig. 4. (a) and (b) are experimental B-scans formed using SBR; (c) is a maximum amplitude projection C-scan formed using SBR. Gaussian smoothing followed by 10-times upsampling was applied to generate these images from the raw images in Fig. 3.

Download Full Size | PDF

 

Fig. 5. (a) and (b) are experimental B-scans formed using DSB; (c) is a maximum amplitude projection C-scan formed using DSB.

Download Full Size | PDF

We now calculate the theoretical view-limited lateral resolution of $1.4f_{\#}\lambda$ for comparison with experimental results. For $f_{\#}=1$, the expected theoretical lateral resolution is $1.4f_{\#}\lambda \approx 117\mathrm{\mu m}$. For this calculation, we used the mean received frequency (averaged with respect to power spectral density) of 17.7 MHz and a speed of sound of $c=1480\mathrm{m}/\mathrm{s}$.

It is a challenge to rigorously characterize the resolution of the SBR method. Using the full width at half-maximum (FWHM) is not appropriate as the uncertainty of the localization is not reflected in the FWHM of points in the SBR image. In addition, the FWHM varies significantly with the discretization of the dictionary and the value of $\tau$.

To address this challenge, we used photographic knowledge of straight wire trajectories in an attempt to provide validation of our SBR resolution estimates, and we used the smallest observable separations as our resolution criterion. We created a linear wire-trajectory model for each of the two wires past the intersection (for transverse locations larger than $\approx 2\mathrm{mm}$), creating the linear fit by using cross sections in the delay and sum beamformed image that clearly resolved the two points. Extrapolating these linear fits beyond the cross sections with points resolved by DSB provided a check on the SBR results. The comparison is shown in Fig. 6. Using these linear models, we estimated the resolving power of SBR as the separation between the linear fit estimated point locations on the last cross section where SBR provides two high-intensity peaks with a lower intensity middle region. This resulted in an estimate of 75 μm for SBR’s resolution.

 

Fig. 6. Comparison between the smoothed SBR image and the point locations expected (vertical lines) using a linear fit of wire trajectories generated from the DSB image. The linear fit vertical lines are at separations of 289, 209, 144, and 75 μm.

Download Full Size | PDF

Looking at cross-sectional slices of DSB B-scans, we found that DSB distinguished points to closer separations than the separations to which it was able to localize points well (see Fig. 7). If we required a localization error of less than 50 μm relative to the linear models, we found that DSB was only able to acceptably resolve the two points at separations of $\approx 240\mathrm{\mu m}$ or more. If we required DSB to only distinguish the two points (produce a cross-sectional line plot with two distinct peaks), then it was able to do this for separations of $\approx 200\mathrm{\mu m}$ or more.

 

Fig. 7. Comparison between the DSB image to the point locations expected (vertical lines) using a linear fit of wire trajectories generated from the DSB image. The linear fit vertical lines are at separations of 305, 241, and 201 μm.

Download Full Size | PDF

The DSB resolution (200–240 μm) observed is significantly worse than the theoretical limited-view value ($\approx 117\mathrm{\mu m}$), possibly due to the presence of lower frequencies. However, on the upper side of the intersection of the DSB C-scan (with transverse location less than 2 mm), we found that DSB resolved points to roughly 120 μm, close to the theoretical limited-view value.

Using SBR increased the localization precision relative to using DSB. The FWHM for the smoothed SBR image was about 25 and 10 μm in the lateral and axial directions, respectively. Note that these values are strongly dependent on the size of the Gaussian kernel used for blurring. DSB had much larger FWHM maximum values of 180 and 120 μm in the lateral and axial directions, respectively. These approximate FWHM values were obtained from B-scans where the two points were clearly distinguished, and the point-spread functions were regular in shape and small, corresponding to approximate best-case FWHM performance. The reduction in FWHM achieved by SBR should be taken in the context of Fig. 6, which suggests that the localization provided by the points in the smoothed SBR image is accurate to within 25 μm.

Using SBR reduced the intensity of artifacts outside high-intensity regions relative to using DSB. In the B-scans shown in Figs. 4(a) and 4(b), SBR has a mean signal level in regions below the half-maximum level of $-50$ and $-51\mathrm{dB}$, respectively. By contrast, in a subview of the same size as that used in Figs. 4(a) and 4(b), in the B-scans shown in Figs. 5(a) and 5(b), DSB has a mean signal level in regions below the half-maximum level of $-18$ and $-17\mathrm{dB}$, respectively.

We presented a photoacoustic imaging reconstruction strategy which incorporates spatial sparsity information during the beamforming process. SBR improved on the lateral resolution of DSB, with its 75 μm experimental lateral resolving power beating the 200–240 μm of DSB by more than a factor of 2.5. SBR’s resolution also beat the standard resolution limit for a view-limited transducer, which was 117 μm for our system. Additionally, SBR enjoyed greater precision in its localization and created images with fewer artifacts outside high-intensity regions. It should be noted that the approach presented should be scalable with transducer frequency and should be applicable to other transducer geometries.

Experimental results were obtained while imaging wires through a scattering medium, hinting at the possibility of applying SBR to extended objects such as vasculature in vivo. Future work should investigate the applicability of the proposed approach to the imaging of more complex objects with lower sparsity, and investigate the robustness of the proposed techniques to tissue heterogeneities and aberration. Future work could also involve adaptive or multi-scale grids, dictionary optimization for reduced coherence, and automation of the process for finding the optimal value of the weighting parameter $\tau$ (following the work in Refs. [15,16]).

This Letter represents, to the best of our knowledge, the first demonstration of SBR to achieve super-resolved PACT images deep in a scattering medium. This is a step towards future application of SBR to areas including immune cell and stem cell tracking, high-resolution imaging of microvasculature, and visualization of sparse neuron firing events.