In deep-tissue photoacoustic imaging, optical-contrast images of deep-lying structures are formed by recording acoustic waves that are generated by optical absorption. Although photoacoustics is perhaps the leading technique for high-resolution deep-tissue optical imaging, its spatial resolution is fundamentally limited by the acoustic wavelength, which is orders of magnitude longer than the optical diffraction limit. Here, we present an approach for surpassing the acoustic diffraction limit in photoacoustics by exploiting inherent temporal fluctuations in the photoacoustic signals due to sample dynamics, such as those induced by the flow of absorbing red blood cells. This was achieved using a conventional photoacoustic imaging system by adapting concepts from super-resolution fluorescence fluctuation microscopy to the statistical analysis of acoustic signals from flowing acoustic emitters. Specifically, we experimentally demonstrate that flow of absorbing particles and whole human blood yields super-resolved photoacoustic images, and provides static background reduction. By generalizing the statistical analysis to complex-valued signals, we demonstrate super-resolved photoacoustic images that are free from common photoacoustic imaging artifacts caused by band-limited acoustic detection. The presented technique holds potential for contrast-agent-free microvessel imaging, as red blood cells provide a strong endogenous source of naturally fluctuating absorption.
© 2017 Optical Society of America
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 hundred micrometers inside tissue due to light scattering. This fundamental limitation originates from the optical wavefront distortions that are induced by scattering, preventing diffraction-limited optical focusing at larger depths. While deeper investigations are possible using either optical-diffusion-based tomography or by combining nonoptical modalities such as ultrasound, the resolution of these deeper-penetrating techniques does not allow microscopic studies at depths larger than a millimeter . 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 thermoelastic stress generation propagate in soft tissue without appreciable scattering and can be detected by an external ultrasonic detector. PA imaging is thus a noninvasive imaging technique that provides images of optical absorption at large depths with a spatial resolution limited by acoustic diffraction. Ultimately, the ultrasound resolution in soft tissues is determined 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 100–200 in practice . For example, at a depth of 5 mm, a resolution of around 25 μm at best may be expected . This resolution is two orders of magnitude inferior compared to the optical diffraction limit and is insufficient for the important tasks of cellular imaging and identification of angiogenic microvessels in the microvasculature .
Recently, it was demonstrated that the acoustic diffraction limit can be surpassed in PA imaging by exploiting temporal fluctuations in PA signals induced by dynamic optical speckle illumination . Inspired by super-resolution optical fluctuation imaging (SOFI) , where blinking fluorescent molecules are used to induce temporal fluctuations of fluorescence, multiple random speckle illuminations were exploited as a source of PA fluctuations . In this approach, the acoustic diffraction limit was surpassed by the statistical analysis of a set of PA images obtained using randomly varying optical speckle patterns.
While impressive proof-of-principle results based on this approach have been reported using tissue phantoms [8,9], it has not yet been applied in practical deep-tissue imaging. Two main reasons for the difficulty in applying speckle-based approaches for deep-tissue imaging are the requirement for a long-coherence laser source and the inherent low amplitude of the PA fluctuations compared to the PA mean signal at large depth, which results from the large difference between the small optical speckle grain size and the acoustic wavelength: at a depth larger than a millimeter, the optical speckle grain dimensions are of the optical wavelength scale (of the order of 1 μm), orders of magnitude smaller than the ultrasound wavelength (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 signal value since the PA signal at each image pixel is the sum of a very large number of uncorrelated fluctuating speckle grains enclosed within the acoustic point spread function (PSF) [10,11]. Resolving the small fluctuations over the large mean signal may be challenging under common signal-to-noise ratio (SNR) conditions .
Here, we present a novel PA imaging technique that provides a resolution beyond the acoustic diffraction limit using a conventional PA imaging system and without relying on any structured or coherent illumination. By exploiting fluctuations in the sample itself rather than external-illumination-induced fluctuations, our approach allows us to overcome many of the difficulties and limitations of the state-of-the-art speckle-based approaches. Specifically, we demonstrate that fluctuations in PA signals that are induced by the natural motion of absorbers lead to super-resolved PA images. In biological tissues, naturally moving endogenous absorbers may consist of red blood cells (RBCs) flowing through the microvasculature, which is an important imaging target in, e.g., the evaluation of drugs that affect angiogenesis . In this context of flowing absorbers, flowing exogenous contrast agents have been recently exploited for super-resolved imaging in medical ultrasounds, either by localization of the sparse distribution of absorbers  or by adapting a SOFI-based statistical analysis . In photoacoustics, sparsely distributed flowing particles have recently been proposed  as a replacement for speckle illumination  for improving visibility and contrast in limited-view PA tomography, and very recently, localization-based analysis has been employed on such a sparse distribution of absorbers for super-resolved PA imaging [15,16]. However, these approaches require a very sparse distribution of absorbers, which severely limits their application, in particular, in contrast-agent-free scenarios. Here, by adapting an advanced SOFI-based statistical analysis framework, we demonstrate that even a nonsparse natural distribution of flowing absorbers such as RBCs can be exploited for super-resolved PA imaging, as well as for simultaneous background reduction, using a conventional PA imaging system. Furthermore, we extend the SOFI mathematical framework to allow processing of complex-valued acoustic signals, which is a crucial step in eliminating PA imaging artifacts that are caused by band-limited acoustic detection and lead to the typical axial oscillations of PA images. By relying on optical absorption, our super-resolution approach possesses the advantage of optical and potentially endogenous contrast of RBCs, which are not available in ultrasound echography.
2. PRINCIPLE AND NUMERICAL ILLUSTRATION
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 [where is the illumination spatial profile and is the distribution of optical absorption] with the PSF of the imaging system. For simplicity, we consider a constant and uniform illumination for each laser shot, and a temporally varying random distribution of absorption, given for each of the laser shots by
In the above model, describes the random distribution of moving optical absorbers at the time of the th laser shot and the geometry of the structure through which optical absorbers are moving, i.e., the target structure for imaging. For PA imaging of blood vessels, would represent the vascular structure [Fig. 1(a), white] and the random distribution of RBCs in the vascular structure at the time of the th laser shot [Fig. 1(a), orange]. Thus, for each laser shot , the acquired PA image may be described by the following equation:
Our approach consists of estimating the target structure beyond the resolution limit of the acoustic PSF by exploiting the temporal fluctuations in the unknown distribution of absorbers flowing through the structure at the th laser shot, without employing deconvolution. Following the original SOFI approach , an th-order statistical temporal fluctuation analysis is performed on the acquired images dataset to obtain a resolution increase without deconvolution. This is accomplished by calculating the th-order temporal autocumulant for each pixel in the PA image set.
is calculated recursively using the previous-order cumulants and the th-order moments :7], the th-order cumulant provides images that involve a convolution only with the th power of the PSF rather than the PSF itself and thus provide a resolution increase for a Gaussian PSF, without deconvolution. Deconvolution could further increase the resolution [6,7] (see Section 5).
As a first illustration of the expected resolution increase of our approach under optimal conditions, we performed numerical simulations of a dynamic random absorbers distribution inside a complex microvasculature structure. Figure 1(a) shows a single distribution of the absorbers (in orange) within the vessel-like target (in white) . For each random absorbers distribution, a PA image was simulated by convolving with a typical, experimentally measured PA imaging PSF of a linear ultrasound array [Fig. 1(c), inset], following Eq. (2). The vessel diameters in the simulations are 15–40 μm, where each absorber has a diameter of 10 μm, matching a typical RBC, and the PSF transverse and axial dimensions are respectively 49 and 61 μm full width at half-maximum (FWHM), representing a transducer with a center frequency of 36 MHz, 70% bandwidth, and an -number of 1.2. The absorbers’ concentration is 50% (surface fraction). Figure 1(c) shows the mean image, which is an estimate of , 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 deep-tissue microvessel imaging tasks . Figures 1(d), 1(e), and 1(f) respectively show the 2nd-, 3rd-, and 5th-order cumulant images calculated from a total of simulated PA images of randomly distributed optical absorbers. The inset in each image displays the corresponding th-power of the PSF . These results illustrate how high-order PA fluctuation analysis can provide a significant resolution increase, with the resolution increasing with the analyzed cumulant order [Figs. 1(b) and 1(g)]. While it is impossible to resolve the individual vessels, even in the variance image in Fig. 1(g), the 5th-order cumulant allows for resolving the individual vessels. As a quantitative measure, we have computed the FWHM of a cross section through a thin isolated vessel [Fig. 1(b)] for the different analyzed cumulant orders and compared it to the resolution of the conventional PA imaging given by the transverse FWHM of the PSF [inset in Fig. 1(c)]. The resulting transverse FWHM is 49 μm for the PSF (black trace), 36.5 μm for the variance analysis (red trace), and 21 μm for the 5th-order cumulant (blue trace), in line with the theoretically expected resolution increase . A similar resolution increase is obtained along the axial dimension, as is quantitatively analyzed in Fig. 4 below.
3. EXPERIMENTAL RESULTS
To experimentally validate our approach, we have performed two proof-of-principle experiments using a microfluidic phantom schematically depicted in Fig. 2. In the first experiment (Figs. 3 and 4), the microfluidic channels were flown with a low concentration solution of absorbing microbeads. In a second experiment, the channels were flown with whole human blood at physiological concentration, as is detailed below (Fig. 5). A schematic of the experimental setup used for both experiments is shown in Fig. 2. A capacitive micromachined ultrasonic transducer array (L22-8v, Verasonics, USA) connected to a multichannel acquisition system (Vantage 256 High Frequency, Verasonics, USA) was used to acquire two-dimensional PA images of flowing absorbing microbeads in microfluidic channels crossing the imaging plane. The sample was illuminated with 5 ns laser pulses at 20 Hz repetition rate (, ) from a frequency-doubled Nd:YAG laser (Spitlight DPSS 250, Innolas Laser GmbH, Germany). For each laser shot, PA signals with a center frequency around 15 MHz were acquired simultaneously on 128 elements of the array (pitch 100 μm, elevation focus 15 mm).
As a first controlled experiment, a flowing suspension of absorbing particles in the microfluidic channels was used as a model of RBCs flowing inside vessels. The circuit was built in polydimethylsiloxane (PDMS) with standard soft-lithography manufacturing technology . The circuit plane was located at the () plane, 18 mm below the ultrasound transducer (Fig. 2). The imaged part of the circuit sample consisted of five identical parallel channels, passing perpendicularly through the imaging plane (only four are sketched in Fig. 2). The center-to-center distance between neighboring channels was 180 μm, and each channel was 40 μm wide (along the direction) and 50 μm high (along the direction); see [Fig. 3(b)]. The circuit was filled with pure water loaded with a suspension of monodisperse spherical 10 μm diameter absorbing particles (microparticles GmbH, Germany). The particle concentration was approximately 9100 particles per cubic millimeter (0.5% volume fraction), corresponding to about 13 particles per channel on average in each PA image, given by the 720 μm thick imaging plane defined by the full elevational width at half-maximum of the PA PSF. This relatively low particle concentration was chosen to obtain a controlled flow in the first proof-of-principle validation experiment, without clogging the microfluidic circuit, while providing a good-enough signal-to-noise ratio. The particle flow was monitored directly with a camera and can be viewed in the supplementary movie Visualization 1. The flow speed within the channels was approximately 25 mm/s. Following the first successful validation experiment with flowing microbeads (Figs. 3 and 4), we performed additional experiments with whole human blood, as is detailed below (Fig. 5).
Experimentally, 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 1500 m/s [neglecting the thin (200 μm thick) PDMS layer overlaying the channels]. Electromagnetic interference pick-up noise from the -switched laser circuitry results in parasitic signals appearing mostly synchronously on all the elements. To reduce such constant PA background features, the mean RF value across the elements was subtracted on the RF data for each time and frame before beamforming. Additional signal postprocessing on the stack of beamformed images was performed to separate the fluctuations of interests (originating from flowing particles) from other remaining sources of fluctuations (such as electronic noise on the RF signals) by use of singular-value decomposition-based filtering. The principle of this processing step, which has been used, for instance, in Ref. , is explained in detail in Section 5 of Supplement 1.
The results of employing our technique using 4000 beamformed images are presented in Figs. 3(c)–3(g) and 3(i)–3(l). Examples for individual raw acquired beamformed images are presented in Fig. S4 of Supplement 1. At the probe center frequency of 15 MHz, the theoretical lateral resolution (along the direction) is about 140 μm, given the theoretical aperture of the probe. Experimentally, this lateral resolution was approximately 200 μm—too coarse to resolve the individual channels [Figs. 3(c) and 3(h)]. The inability to resolve the target structures also with conventional ultrasound echographic imaging is illustrated in Fig. 3(h), where a conventional pulse-echo image of the five channels filled with air obtained with a single plane-wave emission is shown. The unresolved channels appear as a continuous line at (the saturated line at corresponds to the strong reflection caused by the PDMS/water interface). This pulse-echo image corresponds to the optimal situation in terms of SNR and has the exact same resolution as the conventional PA images [Fig. 3(c)].
In contrast to the conventional resolution images, the th-order cumulant analysis [Figs. 3(d)–3(f)] clearly resolves the individual channels, as is also evident in the corresponding one-dimensional profiles shown in Fig. 3(g). These results validate experimentally the resolution increase from flow-induced fluctuations. The lateral resolution was evaluated from the far-left channel to be 132, 108, 90, 84, and 72 μm, respectively, for orders from 2 to 6, closely following the theoretical resolution improvement . 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 do not originate from 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 th-order cumulant are nonlinearly related to the absorption profile, they overweight structures with stronger signals . To compensate for this nonlinear weighting, Figs. 3(i)–3(k) present the th root of the th-order cumulant image. The th root compensates for the nonlinear overweighting of the high-order cumulants and provides a quantity that is linearly related to the absorption profile. Interestingly, even when the 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 exact statistics of the flowing particles, which can be influenced by several parameters, such as the diameters of the blood vessels or the mechanical properties of the RBCs and vessel walls .
Inspecting the PA fluctuation images presented in Figs. 3(i)–3(k), one can clearly see the substantial artifacts in the form of anisotropic axial oscillations. These oscillations result from the bandwidth-limited bipolar impulse response of the ultrasound signals, and thus they are inherent to the acoustic PSF [Fig. 1(c), inset]. In conventional PA and ultrasonic imaging, these typical oscillating artifacts are avoided by displaying the envelope of the complex analytical signal, i.e., the magnitude of the complex-valued images, obtained with a Hilbert transform. Such complex-valued processing is a crucial step in avoiding any misinterpretation of PSF oscillations as real structural patterns. To retrieve similar envelope- or magnitude-related information in the SOFI-based fluctuation analysis, one can not simply perform a Hilbert transform on the cumulant images since they are nonlinearly related to the raw signals. For example, cumulants of even orders are, by definition, nonnegative. Thus, we develop below an analytical signal processing framework that is based on a generalization of the real-valued moments and cumulants to complex random variables . To perform this analysis, complex-valued PA images were obtained by beamforming the complex analytical temporal signals (obtained via a Hilbert transform of the raw real-valued RF signals). Following Eriksson et al. , we then define the th-order complex moment as , where . The complex cumulants of order , , are linked to the complex moments, , by the following formula derived from . For (see the derivation in Supplement 1),
If , the following formula applies:
To demonstrate the superiority of the complex analysis framework, we present a comparison of the complex analysis to the conventional SOFI analysis for both simulation [Figs. 4(a) and 4(c)] and experimental results [Figs. 4(b) and 4(d)]. 3rd-order real and complex cumulants are compared, where the central complex cumulants are taken (i.e., for the odd th order used). By taking the magnitude of the complex cumulants, axial oscillations are effectively discarded, as desired. As detailed in Section 3 of Supplement 1, thanks to the additivity and homogeneity properties of complex cumulants , the complex cumulants of the photoacoustic images involve a convolution of cumulants related to the object with , i.e., with what can be considered an th-order complex PSF. In particular, for even order and , involves a convolution with , which is a nonoscillating real th-order PSF, and the cumulant images obtained are the closest to the object. For odd order, central cumulants with give the best representation of the object. For both even and odd orders, the more differs from , the more grainy the result of the convolution with the complex PSF (see Fig. S1 of Supplement 1). As a quantitative measure for the axial resolution increases in both experiments and simulations, we have measured the FWHM of an isolated thin simulated [yellow line in Fig. 4(c)] and a single microfluidic channel [Fig. 4(d)] for the different analyzed cumulant orders. These widths are reported in Fig. 4(f) and Fig. 4(h), and both follow very closely the theoretically expected resolution increase, as validated with fitted curves. In summary, our SOFI analysis with complex cumulants provides oscillation-free super-resolved PA images.
Our experiment with a relatively diluted suspension of flowing absorbing microbeads (Figs. 3 and 4) was used as an experimental proof-of-principle of the proposed super-resolution approach. While the presented approach is general, an important potential application is imaging of microscopic blood vessels. In this discussion, we address the important parameters that dictate the applicability of the approach to in vivo imaging, explain the limitations of our proof-of-principle experiments, and present preliminary experimental results obtained with flowing whole human blood, proving super-resolution imaging using a microfluidic phantom. The key ingredient of the proposed technique is the estimation of high-order statistics via the computation of th-order cumulants from a stack of photoacoustic images. The principal requirements in order to implement this approach are therefore to be able to measure the fluctuations of interest with a sufficient signal-to-noise ratio and to acquire a large enough number of images to compute sufficiently accurate statistical estimates of the th-order cumulant. As is discussed in Section 2, the characteristic number of images, , required to compute a 3rd–6th-order cumulant is typically of the order of a few thousands to tens of thousands [7,21–23]. While this number can be reduced by using advanced reconstruction algorithms , as a general consequence of SOFI, this sets a practical limitation on the acquisition time and the temporal resolution of the technique, which are both dictated by . The time required to acquire uncorrelated images is dictated by the flow rate and/or the laser repetition rate. The minimum time between consecutive frames that would produce completely uncorrelated pixel values is given by , where is the dimension of the PSF along the flow direction and is the flow speed. In the specific case of our experiments with flowing beads, with the flow direction (with ) perpendicular to the imaging plane and an elevation resolution thickness of about 1 mm, this correlation time was approximately 1/20 s. Thus, the laser pulse repetition rate had to be lowered to 20 Hz to ensure decorrelation from one image to the next (although the maximum available PRF of our laser was 100 Hz). The acquisition time was therefore approximately 3 min to obtain the images shown in Fig. 3. Considering a best-case scenario with a high-frequency transducer with and a flow speed of , characteristic for venules of 50 μm diameter, the acquisition rate of uncorrelated frames would be limited to about 100 Hz. The typical acquisition time to obtain a SOFI image from a few thousand acquisitions is therefore of the order of at least a few tens of seconds. This lack of high temporal resolution is intrinsic and general to the SOFI technique but is, in principle, better than the limit of other ultrasound-based super-resolution techniques  since higher concentrations of emitters can be used. When capillaries of 5–10 μm diameter are considered, the flow speed can be an order of magnitude lower, further limiting the acquisition rate but allowing averaging several pulses for a higher SNR. As a general limitation, the presented technique is thus restricted to structural imaging rather than fast dynamics. Besides temporal consideration, the major practical limitation to implement the technique in vivo is the ability to measure fluctuations induced by flow and to separate those from other sources of fluctuations and fixed absorption. The ability to measure fluctuations that are induced by blood flow depends on the amplitude of these fluctuations, which depends both on the optical absorption coefficient and the concentration and distribution of absorbing RBC as well as on the measurement noise, the dynamic range of the acquisition system, and any static and dynamic background fluctuations. One key point for the applicability of our proposed approach to image blood vessels, which could not be validated using the experiments with absorbing beads, is the ability to exploit the fluctuations induced by real blood flow for super-resolution. It is important to note that the ability to measure fluctuations in a photoacoustic signal due to the flow of RBC in blood vessels is not questionable, as it forms the basis for photoacoustic flowmetry and velocimetry techniques. In particular, the specific question of the spatiotemporal fluctuations induced by flowing RBC has been theoretically and experimentally studied in several recent works [24–27]. Going beyond these recent works, and in order to prove that real blood flow can induce photoacoustic fluctuations that enable super-resolution imaging using a standard photoacoustic system, we have performed a preliminary experiment with the microfluidic circuit flown with whole human blood, having a physiological volume fraction of absorbers (hematocrit ) that is 100 times higher than the microparticles flown in the experiments of Figs. 3 and 4. Figure 5 shows the results of these preliminary experiments. Figure 5(b) shows an image taken with an optical microscope while blood was flown in the channels. Some granularity due to the random RBC distribution can be observed, as is more clearly visible in the higher magnification in the top-right inset. The lower inset shows static RBCs located at the input of the circuit and confirms the expected shape and concentration of RBC within the circuit. A movie of blood flowing through the circuit is available as supplementary information (see Fig. S3(b) of Supplement 1, Visualization 2, and Visualization 3). The blood flow was set to approximately 25 mm/s, which is a realistic value for blood flow in vivo. The laser repetition rate was turned to its maximal value of 100 Hz. As compared to the experiments with microbeads, which allowed well controlled and reproducible results, it was considerably more difficult to insure a stable and controllable behavior of the blood flowing through the channels. The quality of the reconstructed SOFI images were highly dependent on the specific dataset and signal processing—in particular, the choice of the relevant SVD filtering for static and other dynamic sources of fluctuations (see Section 6 of Supplement 1). However, even though the results of these preliminary experiments were sensitive to the experimental parameters, the results shown in Fig. 5 demonstrate that by analyzing a specific dataset containing more than 15,000 images (total acquisition time of 2.5 min), it was possible to visualize the five microchannels flown with whole blood well beyond the acoustic diffraction limit [Figs. 5(c)–5(f)]. These preliminary results thus provide a first proof that real blood flow can induce fluctuations that can be exploited for super-resolution imaging and suggest that it may be possible to translate the presented approach to in vivo imaging. The main challenge in applying the technique in vivo is expected to be the background motion-induced fluctuations of living samples and more complicated imaged structures. Here, we have considered complicated vascular structures in the numerical simulations of Fig. 1, and in-plane and out-of-plane static and dynamic fluctuations were naturally present in our experiments. Motion-induced fluctuations can be minimized mechanically and via careful SVD filtering. While our experiments did not involve optical scattering structures, these should not affect the performance of the technique, as long as they do not induce acoustic scattering or aberrations.
We have demonstrated that flow-induced fluctuations in PA signals can be harnessed to provide a 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 potentially used to image blood microvessels 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 less demanding in SNR requirements for fluctuation analysis, whenever the absorbers’ dimensions are larger than the optical speckle grain size. In addition to the improved resolution, the technique inherently increases the visibility of invisible structures in limited-view PA tomography [11,14] and provide background reduction. Unlike localization-based super-resolution techniques , in SOFI, the signal emitters’ distribution is allowed to be dense, enabling a smaller required number of images and thus a reduction in acquisition time.
The basic statistical autocumulants analysis performed here can be potentially improved by calculating cross cumulants and fluctuations at nonzero time lags, which were shown to eliminate noise-induced correlations [7,13]. We have demonstrated a resolution increase using th-order cumulant analysis. A greater resolution increase is attainable by employing deconvolution on the cumulant images, which could reach an -times resolution increase for an th-order cumulant . The resolution increase in SOFI is fundamentally limited by the highest-order cumulant that could be well estimated, given the number of acquired uncorrelated frames . The maximum available number of acquired frames, , is determined by the sample dynamics and total acquisition time (see Section 4). The number of frames required to achieve a sufficiently high imaging fidelity in SOFI is dictated by the convergence properties of the cumulant considered, which is dependent on many parameters, such as the cumulant order, the particle concentration, the sample complexity, the fluctuation statistics, and the system’s PSF. There is, to the best of our knowledge, no general and simple answer to the important question of the required number of frames to obtain a “satisfactory” th-order SOFI image, which is still an active field of research [21–23]. In the results presented here, the number of frames was determined heuristically to provide a high enough fidelity with the imaged structures. The dependence of the obtained SOFI images on the number of frames is presented in Figs. S4 and S5 of Supplement 1. A further increase in resolution and imaging fidelity and a considerable reduction in the number of required frames, beyond that allowed by SOFI, would be provided by employing advanced compressed-sensing-based reconstruction algorithms . This is made possible by taking into account any a priori knowledge of the imaging system response and any known structural sparsity of the target or statistical fluctuations or correlations.
In this respect, it is interesting to note that independent fluctuations occur between different blood vessels, but the fluctuations are spatiotemporally correlated between pixels along the same vessel, which may be used to smooth-out features along the vessel and to algorithmically improve reconstruction fidelity. 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 also be possible to combine several complex cumulants of the same order for reducing even further imaging and statistical artifacts. 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 [24,25] and related statistical analysis  to provide simultaneous information on flow velocity and particle concentration, in addition to significantly improved resolution, from the same image dataset. While an important potential impact of the technique is in deep-tissue blood-vessel imaging, data from any experiment where natural fluctuations of absorption are present, either due to flow or another source of fluctuations, could be analyzed under the SOFI framework for increased resolution. An example of such sources of fluctuations could be changes in the optical properties within the vessel (e.g., change in blood oxygenation) or in the surrounding tissue (e.g., global motion). In principle, the SVD filtering used for clutter filtering can be utilized to separate the multiple sources of fluctuations  and may lead to imaging specificity through the isolation of the fluctuations of interest.
H2020 European Research Council (ERC) (677909, 681514); Azrieli Foundation; Human Frontier Science Program (HFSP).
We thank Simon Christoph Stein, Prof. Joerg Enderlein, Prof. Shimon Weiss, Prof. Yonina Eldar, Oren Solomon, and Prof. Sylvain Gigan 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. We thank Sylvain Losserand for a fruitful discussion about microfluidic experiments with blood. O.K. was supported by an Azrieli Faculty Fellowship.
See Supplement 1 for supporting content.
1. V. Ntziachristos, “Going deeper than microscopy: the optical imaging frontier in biology,” Nat. Methods 7, 603–614 (2010). [CrossRef]
2. L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science 335, 1458–1462 (2012). [CrossRef]
3. P. Beard, “Biomedical photoacoustic imaging,” Interface Focus 1, 602–631 (2011). [CrossRef]
4. M. Omar, D. Soliman, J. Gateau, and V. Ntziachristos, “Ultrawideband reflection-mode optoacoustic mesoscopy,” Opt. Lett. 39, 3911–3914 (2014). [CrossRef]
5. D. M. McDonald and P. L. Choyke, “Imaging of angiogenesis: from microscope to clinic,” Nat. Med. 9, 713–725 (2003). [CrossRef]
6. T. Chaigne, J. Gateau, M. Allain, O. Katz, S. Gigan, A. Sentenac, and E. Bossy, “Super-resolution photoacoustic fluctuation imaging with multiple speckle illumination,” Optica 3, 54–57 (2016). [CrossRef]
7. T. Dertinger, R. Colyer, G. Iyer, S. Weiss, and J. Enderlein, “Fast, background-free, 3d super-resolution optical fluctuation imaging (SOFI),” Proc. Natl. Acad. Sci. USA 106, 22287–22292 (2009). [CrossRef]
8. E. Hojman, T. Chaigne, O. Solomon, S. Gigan, E. Bossy, Y. C. Eldar, and O. Katz, “Photoacoustic imaging beyond the acoustic diffraction-limit with dynamic speckle illumination and sparse joint support recovery,” Opt. Express 25, 4875–4886 (2017). [CrossRef]
9. T. W. Murray, M. Haltmeier, T. Berer, E. Leiss-Holzinger, and P. Burgholzer, “Super-resolution photoacoustic microscopy using blind structured illumination,” Optica 4, 17–22 (2017). [CrossRef]
10. J. W. Goodman, Speckle Phenomena in Optics: Theory and Applications (Roberts, 2007).
11. J. Gateau, T. Chaigne, O. Katz, S. Gigan, and E. Bossy, “Improving visibility in photoacoustic imaging using dynamic speckle illumination,” Opt. Lett. 38, 5188–5191 (2013). [CrossRef]
12. C. Errico, J. Pierre, S. Pezet, Y. Desailly, Z. Lenkei, O. Couture, and M. Tanter, “Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging,” Nature 527, 499–502 (2015). [CrossRef]
13. A. Bar-Zion, C. Tremblay-Darveau, O. Solomon, D. Adam, and Y. C. Eldar, “Fast vascular ultrasound imaging with enhanced spatial resolution and background rejection,” IEEE Trans. Med. Imaging 36, 169–180 (2017). [CrossRef]
14. X. L. Deán-Ben, L. Ding, and D. Razansky, “Dynamic particle enhancement in limited-view optoacoustic tomography,” Opt. Lett. 42, 827–830 (2017). [CrossRef]
15. S. Vilov, B. Arnal, and E. Bossy, “Overcoming the acoustic diffraction limit in photoacoustic imaging by localization of flowing absorbers,” Opt. Lett. 42, 4379–4382 (2017).
16. X. L. Dean-Ben and D. Razansky, “Localization optoacoustic tomography,” arXiv:1707.02145 (2017).
17. S. Tang and G. Whitesides, “Basic microfluidic and soft lithographic techniques,” in Optofluidics: Fundamentals, Devices and Applications (McGraw-Hill, 2010).
18. C. Demene, T. Deffieux, M. Pernot, B.-F. Osmanski, V. Biran, J.-L. Gennisson, L.-A. Sieu, A. Bergel, S. Franqui, J.-M. Correas, I. Cohen, O. Baud, and M. Tanter, “Spatiotemporal clutter filtering of ultrafast ultrasound data highly increases Doppler and ultrasound sensitivity,” IEEE Trans. Med. Imaging 34, 2271–2285 (2015). [CrossRef]
19. T. W. Secomb, “Blood flow in the microcirculation,” Annu. Rev. Fluid Mech. 49, 443–461 (2017). [CrossRef]
20. J. Eriksson, E. Ollila, and V. Koivunen, “Essential statistics and tools for complex random variables,” IEEE Trans. Signal Process. 58, 5400–5408 (2010). [CrossRef]
21. X. Wang, D. Chen, B. Yu, and H. Niu, “Statistical precision in super-resolution optical fluctuation imaging,” Appl. Opt. 55, 7911–7916 (2016). [CrossRef]
22. S. Geissbuehler, C. Dellagiacoma, and T. Lasser, “Comparison between sofi and storm,” Biomed. Opt. Express 2, 408–420 (2011). [CrossRef]
23. J. D. Müller, “Cumulant analysis in fluorescence fluctuation spectroscopy,” Biophys. J. 86, 3981–3992 (2004). [CrossRef]
24. H. Fang, K. Maslov, and L. V. Wang, “Photoacoustic Doppler effect from flowing small light-absorbing particles,” Phys. Rev. Lett. 99, 184501 (2007). [CrossRef]
25. J. Brunker and P. Beard, “Acoustic resolution photoacoustic Doppler velocimetry in blood-mimicking fluids,” Sci. Rep. 6, 20902 (2016). [CrossRef]
26. J. Brunker and P. Beard, “Velocity measurements in whole blood using acoustic resolution photoacoustic Doppler,” Biomed. Opt. Express 7, 2789–2806 (2016). [CrossRef]
27. T. Bücking, P. van den Berg, S. Balabani, W. Steenbergen, P. Beard, and J. Brunker, “Acoustic resolution photoacoustic Doppler flowmetry using a transducer array: optimising processing for velocity contrast,” Proc. SPIE 10064, 100642M (2017). [CrossRef]
28. Y. Zhou, J. Yao, K. I. Maslov, and L. V. Wang, “Calibration-free absolute quantification of particle concentration by statistical analyses of photoacoustic signals in vivo,” J. Biomed. Opt. 19, 037001 (2014). [CrossRef]
29. B. Arnal, J. Baranger, C. Demene, M. Tanter, and M. Pernot, “In vivo real-time cavitation imaging in moving organs,” Phys. Med. Biol. 62, 843–857 (2017). [CrossRef]