Super-resolution photoacoustic imaging via flow induced absorption fluctuations

TL;DR

This paper introduces a novel super-resolution photoacoustic imaging method that leverages flow-induced absorption fluctuations, surpassing traditional acoustic resolution limits and enabling detailed imaging of micro-vessels without contrast agents.

Contribution

The authors adapt super-resolution optical fluctuation imaging (SOFI) techniques to photoacoustic imaging, generalizing the mathematical framework for complex signals to achieve super-resolution.

Findings

Achieved super-resolved images free from oscillations caused by band-limited detection.

Demonstrated resolution enhancement using naturally fluctuating red blood cells.

Applicable to conventional photoacoustic systems without additional hardware.

Abstract

In deep tissue photoacoustic imaging the spatial resolution is inherently limited by the acoustic wavelength. We present an approach for surpassing the acoustic diffraction limit by exploiting temporal fluctuations in the sample absorption distribution, such as those induced by flowing particles. In addition to enhanced resolution, our approach inherently provides background reduction, and can be implemented with any conventional photoacoustic imaging system. The considerable resolution increase is made possible by adapting notions from super-resolution optical fluctuations imaging (SOFI) developed for blinking fluorescent molecules, to flowing acoustic emitters. By generalizing SOFI mathematical analysis to complex valued signals, we demonstrate super-resolved photoacoustic images that are free from oscillations caused by band-limited detection. The presented technique holds potential for contrast-agent free micro-vessels imaging, as red blood cells provide a strong endogenous source of naturally fluctuating absorption.

Figures

Equations

$${\alpha }_k(\mathbf{r})=f(\mathbf{r})\times g_k(\mathbf{r})$$

$$P{A}_k(\mathbf{r})=\left[I_0\times f(\mathbf{r})\times g_k(\mathbf{r})\right]\ast h(\mathbf{r})$$

$$C_n(\mathbf{r})={\mu }_n(\mathbf{r})-\sum _{k=1}^{n-1}\left(\begin{array}{c}n-1\\ k-1\end{array}\right)C_k(\mathbf{r}){\mu }_{n-k}(\mathbf{r})$$

$$K_{p,q}(\mathbf{r})=m_{p,q}(\mathbf{r})-\sum _{u=1}^p\sum _{v=1}^{q-1}\left(\begin{array}{c}p\\ u\end{array}\right)\left(\begin{array}{c}q-1\\ v\end{array}\right)K_{p-u,q-v}(\mathbf{r})m_{u,v}(\mathbf{r})$$

$$-\sum _{v=1}^{q-1}\left(\begin{array}{c}q-1\\ v\end{array}\right)K_{p,q-v}(\mathbf{r})m_{0,v}(\mathbf{r})-\sum _{u=1}^p\left(\begin{array}{c}p\\ u\end{array}\right)K_{p-u,q}(\mathbf{r})m_{u,0}(\mathbf{r})$$

$$K_{p,0}(\mathbf{r})=m_{p,0}(\mathbf{r})-\sum _{u=1}^{p-1}\left(\begin{array}{c}p-1\\ u-1\end{array}\right)K_{u,0}(\mathbf{r})m_{p-u,0}(\mathbf{r})$$

Full text

1 Introduction

Optical microscopy is an invaluable tool in biomedical investigations and clinical diagnostics. However, its use is limited to depths of not more than a few hundreds micrometers inside tissue due to light scattering. This fundamental limitation originates from the irremediable optical wavefront distortions that are induced by light scattering, preventing diffraction-limited optical focusing at larger depths. While deeper investigations are possible using either tomographic techniques based on optical diffusion or by combining non-optical modalities such as ultrasound, these deeper-penetrating techniques do not allow optical diffraction-limited resolution, and generally do not permit microscopic studies of e.g. cellular structures at depths of more than a millimeter [1]. One of the leading approaches for deep tissue ultrasound-aided optical imaging is photoacoustic (PA) imaging [2, 3]. Photoacoustic imaging relies on the generation of ultrasonic waves by light absorption in a target structure under pulsed optical illumination. The ultrasonic waves that are produced in the absorbing parts of the target via thermo-elastic stress generation, propagate in soft tissue almost without scattering, and are detected by an external ultrasonic detector. PA is thus a noninvasive imaging technique which provides images of optical absorption at large depth with a spatial resolution limited by acoustic diffraction. Ultimately, the ultrasound resolution in soft tissue is limited by the acoustic wavelength used, which is limited by the attenuation of high frequency ultrasonic waves. As a result, the depth-to-resolution ratio of deep-tissue PA imaging is approximately 200 in practice [2]. For example, at a depth of 5 mm a resolution of around 25 $\mu$m at best may be expected, two orders of magnitude above the optical diffraction limit [4]. Such a resolution is insufficient for the important tasks of deep-tissue cellular imaging and for the identification of angiogenic microvessels in the microvasculature [5].

Recently, it was demonstrated that the conventional acoustic diffraction-limit can be surpassed by exploiting temporal fluctuations in PA signals that originates from dynamic optical speckle illumination [6]. Inspired by super-resolution optical fluctuation imaging (SOFI) [7], where blinking fluorescent molecules are used to induce temporal fluctuations of fluorescence, multiple random speckle illuminations were exploited as a source of PA fluctuations [6]. The acoustic diffraction limit was surpassed by statistical analysis of a set of PA images obtained under such set of random, unknown optical speckle patterns. While impressive results based on this approach have been reported in several recent works using tissue phantoms [8, 9], a major practical hurdle prevents the use of optical speckle-based techniques for deep tissue PA imaging. This limit is set by the low contrast of the PA fluctuations amplitude, resulting from the large difference between the small optical speckle grain size to the acoustic wavelength used at depths larger than a millimeter. At such depths the optical speckle grain dimensions are of the optical wavelength scale (of the order of one micrometer), orders of magnitude smaller than the ultrasound wavelength (typically 100 µm for a typical frequency of about 15 MHz). This mismatch in dimensions leads to a small amplitude of fluctuations with respect to the mean value in each image pixel, since the PA signal at each image pixel is the sum of a large number of uncorrelated fluctuating speckle grains enclosed within the large acoustic point spread function (PSF) [10, 11]. Resolving the small fluctuations over the large mean signal is challenging under common signal to noise (SNR) conditions [11].

In this work, we present a novel PA technique that provides a resolution beyond the acoustic diffraction limit without relying on any structured or coherent illumination. By exploiting fluctuations in the sample itself rather than external illumination-induced fluctuations, it allows to overcome many of the difficulties and limitations of the state of the art approaches. Specifically, we demonstrate that fluctuations in the PA signals that are induced by the motion of absorbers leads to super-resolved PA images. In biological tissue, naturally moving endogenous absorbers may consist of red blood cells flowing through the microvasculature, an important imaging target in e.g. the evaluation of drugs that affect angiogenesis [5]. In this context of flowing absorbers, flowing exogenous contrast agents have been very recently exploited for super-resolved imaging in medical ultrasound, either by localization of sparse distribution of absorbers [12], or by adapting a SOFI-based statistical analysis [13]. In photoacoustics, sparsely distributed flowing particles have recently been proposed [14] as a replacement for speckle illumination [11] for improving visibility and contrast in limited-view PA tomography. Here, by employing a SOFI-based statistical analysis framework, we demonstrate that non-sparse flowing absorbers can be exploited for super-resolved background-reduced PA imaging. Furthermore, we extend the SOFI mathematical framework to allow processing of complex-valued acoustic signals, which is an important crucial step in order to discard the PA imaging artifacts that are caused by the typical axial oscillations of the PA band-limited PSF. By exploiting optical absorbers, our simple super-resolution approach possesses the advantage of optical and potentially endogenous contrast, such as the optical absorption of red blood cells, which is not available in ultrasound only imaging. In this context it is important to note that the SOFI-based processing enables super-resolved imaging even with non-sparse distributions of sources [7]. Here we demonstrate these advantages using conventional PA images of flowing absorber, which can be acquired with commercially available systems.

2 Principle and numerical simulations

To explain the source of increased resolution in our approach, we consider the following model of PA imaging: in its simplest form, PA imaging with short pulse illumination may be described as a convolution of the spatial distribution of absorbed optical energy $E_{\mathrm{abs}}(\mathbf{r})=I(\mathbf{r})\times\alpha(\mathbf{r})$ (with $I(\mathbf{r})$ the illumination distribution and $\alpha(\mathbf{r}$) the distribution of optical absorption) with the PSF $h(\mathbf{r})$ of the imaging system. For simplicity, we consider a constant and uniform illumination $I(\mathbf{r})=I_{0}$ for each laser shot, and a temporally varying random distribution of absorption, given for each of the $k=1..N$ laser shots by

[TABLE]

In the above model, $g_{k}(\mathbf{r})$ describes the random distribution of moving optical absorbers at the time of the $k^{th}$ laser shot and $f(\mathbf{r})$ the geometry of the structure through which optical absorbers are moving, i.e. the target structure for imaging. For PA imaging of blood vessels, $f(\mathbf{r})$ would represent the vascular structure (Fig. 1(a), white) and $g_{k}(\mathbf{r})$ the random distribution of red blood cells in the vascular structure at laser shot $k$ (Fig. 1(a), orange). Thus, for each laser shot $k$ the PA image $PA_{k}(\mathbf{r})$ may be described by the following equation:

[TABLE]

Our approach consists of estimating the target structure $f(\mathbf{r})$ beyond the resolution limit of the acoustic PSF $h(\mathbf{r})$, by exploiting the temporal fluctuations in the unknown distribution of absorbers $g_{k}(\mathbf{r})$ flowing through the structure at the $k^{th}$ laser shot. Following the original SOFI approach [7], an $n^{th}$-order statistical fluctuation analysis is performed on the acquired image set to obtain a $\sqrt{n}$ resolution increase without deconvolution. This is accomplished by calculating the $n^{th}$-order auto-cumulant $C_{n}(\mathbf{r})$ for each pixel in the PA image set. $C_{n}(\mathbf{r})$ is calculated recursively using the previous order cumulants and the $k=1..n^{th}$-order moments $\mu_{k}(\mathbf{r})=\dfrac{1}{N}\sum_{k=1}^{N}\left[PA_{k}(\mathbf{r})\right]^{m}$ [7]:

[TABLE]

where $C_{1}(\mathbf{r})=\mu_{1}(\mathbf{r})$ is the mean image, and $C_{2}(\mathbf{r})$ provides the variance image. Following Dertinger et al. [7], the $n^{th}$-order cumulant provides images that involve a convolution only with the $n^{th}$ power of the PSF rather than the PSF itself, and thus provide a $\sqrt{n}$ resolution increase for a Gaussian PSF, without deconvolution. Deconvolution could further increase the resolution [7, 6].

As a first study of the expected resolution increase, we performed numerical simulations of a dynamic random absorption distribution in a complex microvasculature structure. Fig. 1(a) shows a single distribution of the absorbers (in orange) ($g_{k}(\mathbf{r})$) within the vessel-like target (in white) ($f(\mathbf{r})$). For each random absorbers-distribution, a PA image $PA_{k}(\mathbf{r})$ was simulated by convolving $f(\mathbf{r})\times g_{k}(\mathbf{r})$ with a typical, experimentally measured PA imaging PSF of a linear ultrasound array (Fig. 1(c), inset), following Eq.2. Fig. 1(c) shows the mean image, which is an estimate of $h(\mathbf{r})\ast f(\mathbf{r})$, and provides the conventional PA image. As can be appreciated the corresponding resolution is too low to resolve the vessel-like structure, as is the situation in many practical microvessel imaging tasks [5]. Fig. 1(d), (e) and (f) respectively show the $2^{nd}$-, $3^{rd}$- and $6^{th}$-order cumulant images calculated from a total of $N=90000$ simulated PA images of optical absorbers randomly distributed inside the vessel structure. The inset in each image illustrates the corresponding $n^{th}$-power of the PSF $h(\mathbf{r})$. These results clearly illustrate how high-order PA fluctuation analysis can provide a significant resolution increase, with the resolution increasing with the analyzed order $n$, as is illustrated in the one-dimensional cross-section through the images, shown in Fig. 1(b) and (g). While it is impossible to resolve the individual vessels even in the variance image in Fig. 1(g), the $6^{th}$-order cumulant allows resolving the individual vessels. As a quantitative measure, we have computed the full width at half maximum (FWHM) of the images of an isolated thin vessel (Fig. 1(b)) for the different analyzed cumulant orders, and compared it to the transverse FWHM of the PSF (inset in Fig. 1(c)). The resulting FWHM is 9.8 pixels for the PSF, 7.3 pixels for the variance (blue trace), and 4.2 pixels for the $6^{th}$-order cumulant (orange trace), in line with the theoretically expected $\sqrt[]{n}$ resolution increase [7].

3 Experimental results

To experimentally validate our approach we have performed a proof-of-principle experiment whose results are presented in Fig. 3. A schematic of the experimental setup is shown in Fig. 2. A Capacitive Micromachined Ultrasonic Transducer (CMUT) array (L22-8v, Verasonics, Kirkland, Washington, USA) connected to a multi-channel acquisition electronics (High Frequency Vantage 256, Verasonics, Kirkland, Washington, USA) was used to produce two-dimensional PA images of channels crossing the imaging plane. PA signals with center frequency around 15 MHz were acquired simultaneously on 128 elements of the array (pitch 100 $\mu$m, elevation focus 15 mm) for each laser shot. The sample was illuminated with 5 ns laser pulses at 20 Hz repetition rate ($\lambda$ = 532 nm, fluence = 3.0 mJ/cm2) from a frequency-doubled Nd:YAG laser (Spitlight DPSS 250, Innolas Laser GmbH, Krailling, Germany). Channels of a microfluidic circuit filled with a flowing suspension of absorbing particles were used as a model of red blood cells flowing inside vessels. The circuit was built in PDMS (polydimethylsiloxane) with standard soft-lithography manufacturing technology  [15]. The circuit plane was located at the (X-Y)-plane, 18 mm below the ultrasound transducer. The imaged part of the circuit sample consisted of 5 identical parallel channels, passing perpendicularly through the imaging plane (only 4 are represented on Fig. 2). The 5 parallel channels were identical circuit sections arranged in series along the main circuit. To insure temporal decorrelation of the random stream of particle between the 5 sections, as required by the SOFI analysis, intermediate circuit elements with branching were used to shuffle the flow of absorbing particles. The center-to-center distance between neighbouring channels was 180 $\mu$m, and each channel was 40-$\mu$m wide (along the X direction) and 50-$\mu$m high (along the Z direction), (Fig. 3(b)). The circuit was filled with pure water loaded with a suspension of monodisperse spherical 10 $\mu$m diameter absorbing particles (microparticles GmbH, Berlin, Germany). The particle concentration was approximately 9100 particles per mm3, corresponding to about 13 particles per channel on average in each PA image, given by the 720 $\mu$m thick imaging plane defined by the full elevational width at half maximum of the PA PSF.

For each nanosecond laser pulse, corresponding to one "frozen" random distribution of the particles flowing through the channels, the PA radio-frequency (RF) signals were acquired. The laser jitter and pulse energy variations were compensated using a fixed reference absorber located outside the region of interest. Images were reconstructed from the raw recorded acoustic signals, using a conventional delay-and-sum algorithm, assuming a homogeneous sound speed of $1500m/s$. Electromagnetic noise from the laser results in parasitic signals appearing mostly synchronously on all the elements. To enhance the signal-to-noise ratio, the mean RF value across the elements was subtracted on the RF data, for each time and frame before beamforming. The thin (200 $\mu$m-thick) PDMS layer overlaying the channels was neglected in the reconstruction algorithm. To enhance SNR and reduce PA background features, clutter filtering based on spatio-temporal singular value decomposition [16] was performed on multiple series of 500 images ensembles keeping singular values from 2 to 8. Statistical analysis was performed on 4000 beamformed images, and its results are presented in Fig. 3(c-g,i-l). At the probe center frequency of 15 MHz, the corresponding lateral resolution (along the X direction) was about 140 $\mu$m. As result of the thin PDMS layer ultrasound absorption, the lateral resolution (along x direction) was degraded to approximately 200 µm, too coarse to resolve the individual channels (Fig. 3(c),(h)).

The inability to resolve the target structures with conventional ultrasound imaging is illustrated in Fig. 3(h), where a conventional pulse-echo image of the 5 channels filled with water obtained with a single plane-wave emission is shown. The unresolved channels appear as a continuous line at z = 17.96 mm (the saturated line at z = 17.8 mm corresponds to the strong reflection caused by the PDMS/water interface). This pulse-echo image corresponds to the optimal situation in term of SNR and has the exact same resolution as each of the conventional PA images. The mean clutter filtered PA image Fig. 3(c) shows features with strong PA signals, which do not correspond to channels: visual inspection of the circuit after the experiment suggested that various optical absorbers (possibly dust or clusters of absorbing beads) were fixed to the circuit and produced the strong PA background features that were reduced by clutter filtering but not totally suppressed. The $n^{th}$-order cumulant images for n=2, 3 and 6 are shown in Fig. 3(d), (e) and (f), with the corresponding one-dimensional profiles across the line of channels shown in Fig. 3(g), and validate experimentally the expected resolution increase. In addition to the resolution increase, these results also demonstrate that high-order fluctuation PA imaging is a powerful approach in suppressing background artifacts from stationary distributions of absorbers: the strong background features observed on the mean image (Fig. 3(c)), which preclude the observation of the channels, are suppressed in the cumulants images. In the context of biological tissue imaging, this may provide an effective approach to selectively image blood flow with super-resolution. Since the $n^{th}$-order cumulant are non-linearly related to the absorption profile [7], Fig. 3(i-k) presents the $n^{th}$-root of the $n^{th}$-order cumulant image. The $n^{th}$-root compensate for the nonlinear over-weighting of the high order cumulants for the objects with the strongest signals, and a quantity linearly related to the absorption of the object is recovered. It is important to note that even when the $n^{th}$-root is presented, the channels are increasingly better-resolved at higher orders (Fig. 3(l)). This enhanced resolution with increasing order depends on the $n^{th}$-order statistics of the flowing particles and can be more or less effective in practice depending on the actual statistics of the flowing particles.

Inspecting the PA images presented in Fig. 3 (i-k), one can clearly notice the substantial artefacts in the form of anisotropic axial oscillations in the PA signal. These oscillations results from the bandwidth limited bipolar impulse response of the ultrasound signals, and thus the acoustic PSF. In conventional PA and ultrasonic imaging this typical oscillating artefacts are avoided by displaying the envelope of the complex analytical signal or magnitude of complex-valued images. Such complex-value processing is a crucial step in avoiding any misinterpretation of PSF oscillations as real structures. To retrieve a similar envelope- or magnitude-related information in the presented SOFI-based fluctuation analysis, one can not simply perform a Hilbert transform on the cumulant images, which are nonlinearly related to the raw signals. For example cumulants of even orders are by definition non-negative. Thus we develop below an analytical signal processing that is based on a generalization of the real value moments and cumulants definitions for complex random variables [17]. Complex-valued PA images were obtained by beamforming the complex analytical temporal signals (obtained with a Hilbert transform of the raw real-valued RF signals). We defined $n^{th}$ order complex moments as: $m_{p,q}(\mathbf{r})=\dfrac{1}{N}\sum_{k=1}^{N}\left[PA_{k}(\mathbf{r})]^{p}.[PA_{k}(\mathbf{r})^{*}\right]^{q}$ where $q=n-p$. The complex cumulants of order $n=p+q$, $K_{p,q}$, are linked to the complex moments, $m_{p,q}$ by the following formula derived from [17]. For $q>0$:

[TABLE]

If $q=0$, the following formula applies:

[TABLE]

where $p+q$ is the order of the cumulant. There are $n+1$ cumulants for a given order and it can be noted that $K_{p,q}=K_{q,p}*$.

For both simulation (Fig. 4(a,c)) and experimental results (Fig. 4(b,d)), we compare the $3^{rd}$-order real cumulants Fig. 4(a-b) to the complex central ones Fig. 4(c-d), meaning $p=q$ if $n$ is even and $p=q\pm 1$ if $n$ is odd. By taking the magnitude of the complex cumulants, axial oscillations are discarded, as for conventional ultrasound imaging. Indeed, $K_{p,q}(\mathbf{r})$ is related to $h(\mathbf{r})^{p}.h(\mathbf{r})^{*q}$ thanks to the additivity and homogeneity properties [17]. Then $|K_{p,q}(\mathbf{r})|$ is related to $|h(\mathbf{r})^{p}.h(\mathbf{r})^{*q}|=|h(\mathbf{r})|^{p+q}=|h(\mathbf{r})|^{n}$. Cumulants moduli with different $(p,q)$ values for a given $n^{th}$-order give different statistical estimates of the object to the $n^{th}$-power convoluted with the magnitude of the PSF to the $n^{th}$-power. This analysis therefore provides artefact-free sub-acoustic resolution PA images.

4 Conclusion

We have demonstrated numerically and experimentally that flow-induced fluctuations in PA signals can be harnessed to provide resolution beyond the acoustic diffraction-limit. Our technique can be directly implemented on any conventional PA imaging system and with any illumination. The technique may be used to image blood micro-vessels that are otherwise blurred, by exploiting the natural flow of endogenous absorbers such as red blood cells. Unlike previous approaches [6, 11], the technique does not require coherent illumination, and is considerably less demanding in SNR requirement for fluctuations analysis, in cases where the absorbers dimensions are larger than the speckle grain size (i.e. the optical wavelength). In addition to the improved resolution, the technique increases the visibility of invisible structures in limited-view PA tomography [11, 14]. Unlike localization based super-resolution techniques [12], in SOFI, the signal emitters distribution may be allowed to be dense, enabling a smaller required number of images and thus a reduction in acquisition time, which has to be tuned to reach a fair trade-off between spatial resolution, signal-to-noise ratio and time resolution [7].

The basic statistical autocumulants analysis performed here can be potentially improved by calculating cross-cumulants and fluctuations at non-zero time-lags, which were shown to eliminate noise-originated correlations [7, 13]. Advanced image reconstruction algorithms [8, 9] can be employed to significantly improve image quality and spatio-temporal resolution, by taking into account knowledge of the imaging system response and any inherent structural sparsity of the target or correlation in the particle flow. In this respect, it is interesting to note that independent fluctuations occur between different blood vessels, but the fluctuations are spatio-temporally correlated between pixels of the same vessel, which can be used to smooth-out features along the vessel and to algorithmically improve reconstruction fidelity [8]. We have performed the SOFI statistical fluctuations analysis on both real- and complex-valued beamformed images, and shown that the latter is an effective approach to discard the typical axial oscillations of the acoustic PSF. It may be possible to combine several cumulants of the same order for reducing even further imaging and statistical artefacts. Finally, we note that the presented technique is based on statistical analysis of flowing absorbers, as such it may be combined with a Doppler based analysis [18, 19] and related statistical analysis [20] to provide simultaneously information on flow velocity and particle concentration, in addition to the significantly improved resolution, from the same image dataset.

Funding. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grants no. 677909, 681514), and Human Frontiers Science Program. O.K. was supported by an Azrieli Faculty Fellowship.

Acknowledgment. We thank Prof. Yonina Eldar and Oren Solomon for helpful discussions. We thank Mehdi Inglebert, Gwennou Coupier, Danielle Centani, Philippe Marmottant and especially Sylvie Costrel for their help in the design and fabrication of the microfluidic samples.