Single laser-shot super-resolution photoacoustic tomography with fast sparsity-based reconstruction.

Recently, ℓ 1 -norm based reconstruction approaches have been used with linear array systems to improve photoacoustic resolution and demonstrate undersampled imaging when there is sufficient sparsity in some domain. However, such approaches have yet to beat the half-wavelength resolution limit. In this paper, the ability to beat the half-wavelength diffraction limit is demonstrated using a 5 MHz ring array photoacoustic tomography system and ℓ 1 -norm based reconstruction approaches. We used the array system to image wire targets at ≈ 2 - 3 cm depth in both intralipid scattering solution and water. The minimum observable separation was estimated as 70 ± 10 μ m , improving on the half-wavelength resolution limit of 145 μ m . This improvement was demonstrated even when using a random projection transform to reduce data by 99 % , enabling substantially faster reconstruction times. This is the first photoacoustic tomography approach capable of beating the half-wavelength resolution limit with a single laser shot.

Introduction

Photoacoustic tomography is a relatively recent imaging modality that provides optical absorption contrast with acoustic resolution. While optical resolution photoacoustic imaging provides micron-scale resolution at superficial depths, ultrasound diffraction limits spatial resolution at depths beyond an optical transport mean free path. Abbe's diffraction limit of , where is the received sensing wavelength and is the acoustic numerical aperture, has long stood as a lower bound on resolution of wave-based imaging systems. In the case of half-view (or more) detector geometries this limit becomes , where is the center wavelength.

Significant widespread attention has been given to approaches capable of beating this diffraction limit. One approach to super-resolution is to modify target interrogation to sharpen the point spread function (PSF). Examples of this approach in a fluorescence microscopy setting include stimulated emission depletion (STED) microscopy [1] and structured illumination microscopy [2]. In an optical resolution photoacoustic microscopy setting, this approach generally requires additional laser excitation, and harnesses nonlinear effects such as optical absorption saturation [3], Grüneisen relaxation [4], or reversible photoswitching [5]. However, these approaches are only appropriate for superficial optical imaging.

Super-resolution can also be obtained by super-localizing point sources, without necessarily seeking to sharpen the PSF of the interrogation system. This is done by estimating the centroids of sufficiently separated signal sources. Examples of this approach in a fluorescence microscopy setting include stochastic optical reconstruction microscopy (STORM) [6] and photoactivated localization microscopy (PALM) [7]. Localization-based methods have also been explored in the context of ultrasound imaging [8], [9].

In the acoustic resolution photoacoustic imaging context, localization-based methods have achieved super-resolution by localization of time-varying point sources [10], [11], and by localization of flowing absorbing particles [12], [13]. These localization-based methods require that signal sources be sufficiently well separated. Consequently, these approaches generally require multiple acquisition events, with sufficient sparsity in each frame. Such approaches require long acquisition times and may be challenged by tissue motion.

Another approach being investigated to improve spatial resolution compared to diffraction-limited imaging is sparsity-based reconstruction (SBR). SBR poses reconstruction as an optimization problem and incorporates prior information about the phantom in a regularization term that promotes sparsity in a specified domain. The inclusion of prior information then makes it possible to surpass traditional resolution limits. In contrast to the super-resolution methods reviewed above, SBR does not necessarily require external agents, multiple acquisitions, or other target interrogation modifications. In addition, instead of requiring signal sources to be well-separated, as required by super-localization methods, SBR requires the target to be sufficiently sparse in a known domain. In an ultrasound setting, SBR has enabled improved resolution of point targets [14], although to our knowledge it has not yet been demonstrated to enable resolution below the limit. In a photoacoustic setting, application of SBR has primarily made use of sparsity different from point-target sparsity (e.g. gradient sparsity) and often seeks to improve reconstruction while undersampling [15], [16], [17], [18], [19], [20], [21], [22]. In addition, some work applying SBR to improve resolution begins by applying a backprojection or delay and sum algorithm [17], [23], a step which may lose information. For these reasons, much of the prior work on photoacoustic SBR tells us little about the limits of the technique to resolve point targets.

Some recent work has used a SBR approach optimized for resolving point targets but did not surpass the resolution limit. This recent work [24], [25], [26], [27], [28] makes use of linear arrays and a dictionary designed to sparsify point targets and performs reconstruction directly in the channel-data domain. The approach in [24] was able to reduce artifacts and obtain higher signal-to-noise ratio. Prior work from our group [25] used a high-frequency linear array transducer and a SBR approach to beat the aperture-limited diffraction limit. The work in [26] used a similar approach to surpass the aperture-limited diffraction limit with a focus on using a sparse linear array. The approach in [27] also utilized a similar reconstruction approach, with a focus on reducing the data used for reconstruction. The work in [28] incorporated attenuation compensation into the SBR algorithm. However, none of these approaches demonstrated resolution below the limit.

While SBR has been shown to improve spatial resolution, its high computational burden (and hence slow reconstruction speed) poses a barrier to practical application. The computational problem arises from the fact that, as a model-based reconstruction approach, SBR involves the manipulation of very large matrices. Various approaches have been explored for reducing the size of the data to be manipulated in model-based reconstruction, such as exploiting symmetries to reduce model size [29], splitting the reconstruction problem into several smaller problems [30], [31], using a GPU to perform matrix calculations in parallel and on-the-fly [32], [33], reducing sampling rate by discarding negative frequencies [34], and using hardware or software to obtain measurements corresponding to projections of scrambled versions of the original data to be measured [35], [27], [36], [37]. However, to our knowledge, the only work exploring data-reduction while obtaining super-resolved photoacoustic images is [26], which used as little as of the available data by using a subset of the transducer elements available. To date, the robustness of SBR to retain super-resolution capability even under the application of data reduction schemes remains largely uncharacterized.

In addition to computational burden, the potential presence of non-sparse signal sources is an additional challenge for the practical implementation of SBR. However, to our knowledge, none of the recent work applying SBR optimized for point-target reconstruction has sought to characterize reconstruction performance in the presence of less-sparse signal sources [24], [26], [27], [25]. Instead, testing of this SBR approach so far has been restricted to the imaging of point targets.

In this work we sought to demonstrate resolution surpassing the limit with a single laser shot, while taking initial steps towards characterizing the robustness of SBR to data reduction (for reconstruction acceleration) and to the presence of non-sparse signals. Given that a full-view tomography system can achieve a resolution close to the limit to begin with, we hypothesized that SBR methods could surpass the half-wavelength limit when applied to full-view tomography data. To explore the robustness of SBR to data reduction, we performed SBR using a randomly projected version of the channel data, and also performed SBR in several smaller quadrants partitioning the field of view. Finally, to explore robustness to the presence of less-sparse signals, we imaged a phantom containing a less-sparse target, as well as a sparse target. We found that SBR could surpass the half-wavelength diffraction limit by roughly a factor of two, even when reducing the size of the data by using a random projection approach. This enabled rapid reconstruction (, corresponding to -fold acceleration) of super-resolved images. Our results demonstrate resolution beyond the limit with a single laser shot and indicate at least some robustness of the method to the presence of less-sparse signal sources and to the application of data reduction techniques.

Theory

SBR reconstructs an image by optimizing an objective function that incorporates a priori sparsity information, making use of ideas from compressive sensing (see [38] for an overview). To define the optimization problem, we begin by modeling the imaging system as , where is a vector of optical absorption coefficients for a collection of points in space, is a system matrix (“dictionary”) with columns corresponding to photoacoustic responses from these same points in space, is a noise vector, and is the observed photoacoustic channel data (see Fig. 1). We include a spatial sparsity promoting -norm regularization term, and an -norm term which ensures an approximate match between the observed channel data and the expected channel data under the estimated absorber intensity vector . Finally, we introduce a parameter which controls the extent to which sparsity is encouraged in the reconstructed image. The SBR image is then obtained by solving the following optimization problem:

Fig. 1

In a noiseless setting, the received pressure-over-time channel data is the superposition of responses from individual absorbers. Generating an image corresponds to forming an estimate of the location and strength of these absorbers and can be done using methods including sparsity-based reconstruction and back-projection.

One drawback of the SBR implementation described above is that it requires working with a very large dictionary matrix , which increases rapidly in size as the reconstruction grid is refined or as the reconstruction area is increased. We hypothesized that we could reduce computational burden by modifying the system to a system , where is some random (fixed) projection matrix to a lower dimension. The optimization problem then becomes:

Calculating can be done rapidly, while on the other hand calculating is intensive, but only must be done once. This modification effectively changes the sensing matrix to , which is much smaller and consequently requires much less memory and computation to manipulate. We can justify this approach by noting that random projection matrices have excellent theoretical properties that allow for recovery of sparse signals after they have been applied (see e.g. [38]). We hypothesized that this projection approach should enable high-speed, low-memory, acceptable-quality SBR imaging.

Materials and methods

To estimate the minimum separation at which SBR can resolve point sources, we imaged a crossed-wire target at successive cross-sections (see Fig. 2), using a ring array to detect photoacoustic signals generated by pulsed illumination. Experimental channel data was collected using a pulsed nanosecond Nd:YAG laser (Surelite OPO Plus, Continuum), a Imasonic ring array (with 256 elements spread over ), and a programmable ultrasound system (Vantage 256, Verasonics, US). The crossed-wire target was constructed from aluminum wires in diameter (ALW-29S, Heraeus) angled at relative to each other. We first imaged the target through - water, and then repeated the experiment at a similar depth in intralipid tissue mimicking solution. Cross-sectional images of the target were taken at steps. To improve the signal-to-noise ratio, we collected data repeatedly at each cross-section ( per cross-section in water and per cross-section in intralipid) for selectable levels of signal averaging, including for reconstructions with no averaging.

Fig. 2

To characterize resolution, we imaged two converging wires over a range of cross-sections. This experiment was performed both in water and in intralipid.

As an initial step towards characterizing the applicability of SBR for obtaining super-resolution in less-sparse contexts, we also imaged a less spatially-sparse phantom, containing both in-plane (as opposed to through-plane) crossed wires and a single through-plane wire. The same laser, ring array, ultrasound system, and aluminum wire were used as in the resolution sub-experiment described above. We imaged the target through of water with times averaging.

To generate the dictionary matrix , we made use of the dictionary refinement procedure described in [25]. This procedure forms a dictionary by starting with an initial approximate dictionary (estimated using simulation), and then modifying this dictionary to sparsify estimated absorber locations (for an experimental reference image of a point in water) at a low value. The responses from other points are then estimated by delaying and scaling this calibrated point response. For numerical stability, we normalized to have the maximum element normalized to unity magnitude, so that after normalization . The location of the points corresponding to the dictionary entries were selected by taking points generated by a uniform random distribution, and then moving these selected points apart to avoid clustering. This was done by applying displacements proportional to the reciprocal of the squared distance between points until clustering was visually minimized. Future work could select these points using a regular grid, or utilize an automated approach such as perturbation of a regular grid to reduce clustering, such as in [39].

We solved the SBR optimization problem using the L1-homotopy package [40], as we found it to provide fast and high-quality reconstruction performance relative to other SBR packages available. To accelerate the speed of reconstruction and reduce memory usage, we broke this reconstruction into two steps. We began by performing an initial reconstruction using a coarse dictionary. Working from the resulting image, we then formed a second finer dictionary, omitting spatial locations more than some small radius from locations with estimated nonzero signal. The radius used was 50 in the experiment where we imaged two wires in cross-section in water, and the radius used was 25 in the other experiments. We refer to the discarded set of points far from the initial phantom reconstruction and having zero amplitude as the non-relevant point set. For reconstruction, the value of on each reconstruction step was selected manually (to maximize qualitative image quality) and held constant across all cross-sections. A final B-scan image was then formed by using a kernel smoothing operator (with respect to a fine regular grid) which calculated image intensity at a point as a distance-weighted average of nearby intensities within a specified radius of 25 . All reconstructions were performed on a system with a Ryzen 3 1300X Quad-Core 3.5 GHz Processor and 16 GB of RAM.

We now detail the specific parameter values used for reconstruction. In the first sub-experiment (with through-plane wires) when imaging in intralipid, we set and used a dictionary with 2500 points for initial reconstruction, and then set and used a dictionary with 22500 points (prior to the removal of the non-relevant point set) for finer reconstruction. When imaging in water we set and used a dictionary with 2500 points for initial reconstruction, and then set and used a dictionary with 22500 points (prior to the removal of the non-relevant point set) for finer reconstruction. In the second sub-experiment (with in-plane wires) when reconstructing the entire field of view at once, we set and used a dictionary with 2500 points for initial reconstruction, and then set and used a dictionary with 22500 points (prior to the removal of the non-relevant point set) for finer reconstruction. We also performed reconstruction in the second sub-experiment using one quarter of the field of view at a time, with the aim of further testing robustness to less-sparse background signals. For each reconstruction quadrant we used and a dictionary with 2500 points.

In addition to reconstructing the data using SBR, we also used a simple form of back-projection (BP) to provide a basis for comparison. The BP reconstruction algorithm used is the universal back-projection algorithm [41] with the derivative term omitted. In addition, to reduce the level of generated artifacts, we truncated negative absorbance values. It should be noted that further reduction of oscillatory BP artifacts can be obtained by applying deconvolution, although we did not do this in this paper.

To test the ability of the projection approach to provide acceptable-quality images at accelerated reconstruction rates, we applied the projection approach to a cross-section close to the SBR resolution limit in the intralipid sub-experiment. The random projection matrix was generated by drawing each entry from a standard normal distribution. We used to project from a space with 256256 dimensions to a space with 2078 dimensions. To reconstruct quickly in a phantom independent way, we used only the first step in the reconstruction process described above. That is, we omitted the second step that uses a dictionary determined using an initial rough reconstruction.

Results

The reconstructed experimental images of the through-plane crossed-wire phantom are shown in Fig. 3. The C-scan images were generated by taking maximum amplitude projections of a collection of B-scans, and then upsampling by a factor of 10. We observe that SBR resolved the two wires to a smaller separation than BP, as well as providing reduced background signal.

Fig. 3

In both intralipid (a–f) and water (g–l), sparsity-based reconstruction (SBR) (d–f, j–l) was able to resolve two wires to a closer separation than back-projection (BP) (a–c, g–i). All reconstructions in this figure were generated using the full data observed.

To quantify the performance of SBR, we compared its estimated wire separation to that estimated by BP, as shown in Fig. 4. The reported separation on a given frame was calculated as the distance between the two highest intensity points in the reconstructed image, while requiring that the second highest intensity point be more than or from the first point in the case of SBR and BP, respectively. The failure point of a method was assessed qualitatively by observing when the image generated no longer contained two distinct peaks. In the case of BP, we also report the “half-maximum resolution”. We define this as the BP peak-separation just before the two BP intensity peaks are no longer separated by a dip in intensity (in a maximum amplitude projection) to half the value of the average of the two peaks. The half-maximum resolution of SBR is not plotted, as it is roughly equal to the separation between the points at SBR failure.

Fig. 4

The wire separations reported by back-projection (BP) and sparsity-based reconstruction (SBR) in both intralipid (a) and water (b) correlate well, and follow the expected linear trend for the crossed-wire phantom. In each case SBR was able to resolve the two wires at separations below the half-wavelength limit.

In both the water and intralipid case, we observe in Fig. 4 a strong correlation between the separations reported by the two reconstruction approaches up until BP fails to resolve the two wires. Beyond the separations at which BP fails, SBR provides separation estimates that are consistent with a roughly linear rate of reduction in separation, which is reasonable for the crossed-wire target. The source of the jump in separations seen in Fig. 4 (a) is not currently well understood, but could simply be an experimental artifact due to vibrations causing small movements in the wire during the imaging procedure.

We estimated the final separation of the points prior to SBR failure by using a linear fit on BP-estimated separations up to BP failure. This approach produced point separation estimates of and prior to SBR failure in the intralipid and water experiments, respectively. In each case, these estimates for wire separation prior to SBR failure are substantially below the half-wavelength resolution limit corresponding to the center frequency assuming speed of sound in water.

We also wished to determine whether random projection could be used to accelerate reconstruction while preserving super-resolution. As illustrated in Fig. 5, we were able to use random projection to accelerate the reconstruction process by a factor of , while retaining the ability to separate targets closer than .

Fig. 5

We found that using a random projection matrix accelerated sparsity-based reconstruction while preserving super-resolution in the intralipid sub-experiment. The reconstruction times listed are those required to perform optimization after forming the dictionary matrix and calculating any projections.

To explore SBR resolution with variable levels of noise, we additionally used only , and averaging to reconstruct images of Fig. 3. These results are shown in Fig. 6.

Fig. 6

We explored the robustness of sparsity-based reconstruction resolution to variable levels of noise, achieved with (a), (b), and (c) averaging. In intralipid, super-resolution was obtained with a single laser shot.

SBR was able to successfully reconstruct the less-sparse phantom containing in-plane wires, as show in Fig. 7 (b). We note that SBR largely recovers the in-plane wires, albeit with some gaps, and that it also localizes the through-plane wire. This was achieved even when reconstructing one quarter of the field of view at a time (Fig. 7 (c)).

Fig. 7

Sparsity-based reconstruction (SBR) successfully imaged a less-sparse phantom consisting of in-plane crossed wires and a through-plane wire. SBR succeeded when used to reconstruct the entire field of view at once (b), and when used to reconstruct one quarter of the field of view at a time (c). SBR reduced ringing compared to back-projection (BP) (a) but failed to reconstruct portions of the in-plane wires.

Discussion

We found that SBR was able to resolve points separated by a distance roughly half of the resolution limit. To our knowledge, this is the first time that this limit has been surpassed in photoacoustic tomography without making use of super-localization of moving absorbers. This could be important for imaging structures where there is no motion (e.g. micro-metastases). The resolution improvement ratio, defined as the ratio of the obtained super-resolution to the diffraction resolution limit, was in water and in intralipid. These results using a ring array system offer a larger resolution improvement than that reported in prior work using linear arrays ( in [26], and in [25]).

The super-resolution results obtained here were achieved using sparse targets in an idealized imaging environment. Compressive sensing performance in general depends on the properties of the sensing matrix , on the sparsity level, and on the noise level of the system. So, the best resolution obtainable will be context specific. That is, we can expect it to vary with the geometry and impulse response of the imaging system used, the noise level, the level of sparsity actually present in the imaged target, the effective sparsity obtained by using approximate point-response estimates, the size of the projection matrix used for data reduction, and the locations in space selected for reconstruction. Future work could seek to precisely characterize the impact these parameters have on the resulting SBR images. Importantly for application to in vivo imaging, we can expect performance to degrade as the imaging target becomes less sparse or as our ability to characterize the (potentially spatially varying) impulse response of the imaging system decreases. For applications involving sources of signal that are not spatially sparse, it may be appropriate to reconstruct with respect to a different sparsifying prior. In this case, performance may be improved by using a total variation prior or a linear combination of priors that includes spatial sparsity. Future work could explore the relative performance of SBR and traditional reconstruction approaches in these more challenging contexts.

Besides noise, an additional source of uncertainty in our super-resolution measurements was the selected pseudorandom reconstruction locations generated by a user-supervised algorithm. We explored this briefly for the transverse slice shown in Fig. 3 (e). Using five different realizations of the reconstruction locations, we found that the estimated distance between the two wires in this transverse slice varied by 7.4 . This along with the uncertainty in localization over multiple noise realizations is contained within the uncertainty we report in the abstract.

The high degree of computational burden associated with SBR poses another challenge to its practical implementation. The computational burden increases with the fineness of the dictionary used and the size of the field of view. By using a random projection operator to reduce the size of the dictionary matrix by , we reduced reconstruction time to per frame, while preserving super-resolution capability. This corresponded to a speed up by a factor of compared to when random projection was not used. While results did depend on the projection used, even an aggressive random projection allowed us to create super-resolved images, as illustrated in Fig. 5. We also found that we could reconstruct a larger field of view by performing several reconstructions independently over a collection of smaller areas. In the particular case shown in Fig. 7, each quadrant took to reconstruct with a dictionary containing 2500 point responses, but a dictionary (roughly 20 GB) would not fit in the RAM of the computer used for reconstruction. Both random projection and piecewise reconstruction may help enable faster SBR imaging, or help enable SBR imaging with larger fields of view or in three dimensions.

In regards to Fig. 6, note that the results in intralipid indicate the ability to achieve super-resolved images with a single laser shot. This is in contrast to super-localization approaches which typically require thousands of laser shots. It was seen that different realizations of the random projection matrix as well as different experimental noise realizations could impact image reconstruction and lead to failure of super-resolution and should be investigated in future work.

Conclusion

We implemented a photoacoustic tomography ring-array system with sparsity-based reconstruction to demonstrate super-resolution imaging with a single laser shot. We found we were able to experimentally resolve points with separation of roughly half the half-wavelength resolution limit ( vs. ). By making use of a random projection matrix, we were able to accelerate reconstruction to per frame while preserving super-resolution. In addition, we found that our SBR implementation optimized for point targets was able to generate a reasonable image in the presence of a less-sparse target, and even when reconstructing only a quarter of the field of view at a time. This suggests SBR has some robustness to non-sparse background signals. Both data reduction approaches explored (sub-region piecewise reconstruction and random projection) may help enable SBR imaging in contexts with larger field of view or in three dimensions. Future work may explore whether SBR can achieve super-resolution in less ideal contexts, where it is harder to form a high-quality dictionary of point responses and there is substantial non-sparse background signal.

Declaration of competing interest

Roger Zemp is co-founder and shareholder of illumiSonics Inc. and CliniSonix Inc., which, however, did not support this work.

Declaration of Competing Interest

The authors report no declarations of interest.