1. INTRODUCTION

Photoacoustic (PA) imaging is a widely spread biomedical imaging modality that makes use of ultrasound waves resulting from light absorption [1]. In acoustic-resolution PA imaging, arrays of ultrasound detectors are usually used to record PA signals generated by absorbing structures. PA reconstruction algorithms (such as delay-and-sum or backpropagation algorithms) are then used to provide images from a set of PA signals. PA images may, however, present artefacts that arise from coherence of PA waves and characteristics of the detection system such as geometry and frequency bandwidth. Limited-view artefacts occur when strongly anisotropic absorbing structures (such as blood vessels) coherently emit PA waves in some preferential directions, and these waves cannot be detected due to the finite aperture and/or directivity of the ultrasound detectors. In such a case, some parts of the absorbing structure are not visible on the reconstructed PA image. Visibility artefacts also occur with resonant ultrasound detectors, which filter out low-frequency components of PA waves emitted by large absorbers (large as compared to the detection wavelength range). For instance, for large blood vessels, only the vessel boundaries may be visible in PA images.

Various solutions have been proposed to overcome the limited-view problem with finite-aperture detectors. One approach consists in enhancing the detection aperture by performing multiple PA acquisitions at different angles. This can be done by rotating the ultrasound probe around the object [2] or by spinning the object itself [3]. The detection aperture can also be augmented by placing acoustic reflectors at edges of the imaging zone [4,5]. Alternatively, the imaged object and the ultrasound probe can be placed inside a reverberating cavity [6–8]. This allows detecting PA waves emitted in all possible directions, even when only a single-element probe is used [8].

All the mentioned techniques require that the whole range of imaging angles is accessible for detection. However, this is not necessarily feasible in a real clinical environment. In addition, mechanical scanning of the detectors or sample increases the acquisition time and, therefore, degrades the temporal resolution. Methods that do not require detecting waves emitted in all directions have also been proposed. Subwavelength sparsely distributed absorbing particles act individually as isotropic PA sources. This can been exploited to remove limited-view artefacts: Dean-Ben et al. demonstrated that visibility in the limited-view scenario can be improved by means of the localization approach [9] or by using a nonlinear combination of tomographic reconstructions representing sparsely distributed moving particles [10]. Unless sparsely distributed contrast agents are used [11], these techniques remain limited in practice, as blood consists of a concentrated suspension of red blood cells. One more approach consists in heating tissue locally with a focused ultrasound beam and thus generating artificial PA sources via the temperature dependence of the Gruneisen coefficient [12]. By scanning the focused ultrasound beam across the sample and accumulating the resulting PA images, the whole object is reconstructed. However, this approach is limited in practice by the safety thresholds of focused ultrasound. Similar to all other aforementioned approaches, this method requires a very long acquisition time leading to a low temporal resolution. Recently, deep-learning-based approaches have also been proposed as a way to palliate visibility artefacts from limited-view detection in photoacoustic imaging [13–15]. While powerful, deep-learning-based imaging methods also have limitations. These include the availability of a training set and the performances inherently limited to unseen samples close enough to those of the training set, which is a severe practical limitation.

An approach based on multiple-speckle illumination was proposed by Gateau et al. [16]. In this experimental work, a random intensity distribution of speckle patterns, which changed from pulse to pulse, induced fluctuations in each pixel of the corresponding series of PA images. It was demonstrated experimentally with a free-space-generated series of speckle patterns that a second-order fluctuation image provided a faithful representation of the absorbing distribution, free of limited-view and limited-bandwidth artefacts. However, exploiting optical speckle illumination for imaging tissue at depth turns out to be challenging since the fluctuation signal is expected to be very small when the speckle grain size is of the order of the optical wavelength. As a consequence, this technique has so far never been demonstrated in a close-to-clinical environment. Moreover, the work by Gateau et al. [16] provided no clear theoretical explanation for the visibility enhancement in fluctuation imaging with multiple-speckle illumination. Since that work, there have been other studies based on PA signal fluctuations. In particular, such fluctuations were deployed to obtain super resolution in PA imaging [17–20]. Most generally, PA fluctuation imaging utilizes some randomness in PA generation to provide enhanced images as compared to conventional PA imaging. While it has been successfully demonstrated experimentally for super-resolution imaging (induced either from multiple-speckle illumination [17–19] or from random distributions of absorbing particles [20]) and for visibility enhancement (with multiple-speckle illumination [16]), there is, however, no comprehensive theoretical description of fluctuation PA imaging to date.

In the first part of this work, we propose a unified theoretical framework relevant to both fluctuations induced by random illumination patterns and fluctuations induced by random distributions of absorbing particles. The theoretical considerations are validated in Monte Carlo finite-difference time-domain (FDTD) simulations of PA generation in a random medium. In particular, we explain why visibility artefacts are absent in second-order fluctuation images. This property turns out to be completely independent of the origin of fluctuations. In addition, our theoretical results provide a quantitative comparison between fluctuations from multiple-speckle illumination and fluctuations from random distributions of absorbing particles. In the second part, we demonstrate experimentally that fluctuations induced by blood flow can be exploited under physiological conditions to produce PA fluctuation images free of visibility artefacts, as opposed to conventional PA imaging. Two-dimensional images of vessel-mimicking phantoms flowing with blood at physiological concentrations are obtained as a first proof-of-concept demonstration. Finally, the approach is validated in vivo by performing 3D imaging of a chicken embryo.

2. THEORY

In this section, we propose a unified theoretical framework for PA fluctuation imaging applicable to both fluctuations induced by random illumination patterns and fluctuations induced by random distributions of absorbing particles. We first provide theoretical expressions for the mean PA image and the fluctuation PA image. These expressions are then used to explain why the fluctuation images do not possess the visibility artefacts observed on the mean images (equivalent to conventional PA images). Finally, we discuss the dependence of the amplitude of the fluctuation image on characteristic statistical properties (mean, variance, characteristic size) of the random process inducing fluctuations (illumination patterns or distribution of absorbing particles). In particular, we compare the case of randomly generated multiple-speckle illumination to that of randomly distributed red blood cells, and we discuss quantitatively why flowing red blood cells lead to PA fluctuations much larger than those produced by multiple-speckle illumination deep in tissue.

A. Analytical Predictions

1. Theoretical Framework

PA imaging provides images of the absorbed energy density $\alpha$ that we describe in this work as

From the most general perspective, ${g_k}({\bf r}) \times f({\bf r})$ thus describes the spatial distribution of absorbed energy. We assume that ${g_k}({\bf r})$ is a spatially stationary and isotropic random process. Its relevant statistical properties in the context of our work are the mean $\eta = \langle {g_k}({\bf r}{) \rangle _k}$ and the variance $\sigma _g^2 = \langle ({g_k}({\bf r}) - \eta {)^2}{ \rangle _k}$ as first-order statistical properties and its autocovariance

For multiple-speckle illumination, ${V_g}$ (resp. ${D_g}$) is of the order of the characteristic volume (resp. size) of a speckle grain. For random distributions of monodispersed absorbing particles, ${V_g}$ (resp. ${D_g}$) is of the order of the particle characteristic volume (resp. size).

We assume that each PA reconstruction ${A_k}({\bf r})$ corresponding to laser shot $k$ can be expressed as the convolution between ${\alpha _k}({\bf r})$ and the point spread function (PSF) $h({\bf r})$ of the imaging system:

Complex-valued reconstruction provides a straightforward way to remove radio-frequency oscillations by using the modulus of the complex-valued reconstruction as the final image, and it is a widely used approach in medical ultrasound imaging. In the context of photoacoustic imaging, an illustration of the differences between real-valued-only images and complex-modulus images may be found in a previous publication from our group [20].

2. Mean Photoacoustic Image

We consider the mean PA image defined by $E[A]({\bf r}) = \langle {A_k}({\bf r}{) \rangle _k}$, where ${ \langle . \rangle _k}$ denotes an ensemble average (average over laser shots in practice). As $\eta = \langle {g_k}({\bf r}{) \rangle _k}$, $E[A]({\bf r})$ can be straightforwardly expressed as

3. Photoacoustic Fluctuation Image

In the context of this work, we define the PA fluctuation image as a second-order fluctuation image, defined as the square root of the variance image ${\sigma ^2}[A]({\bf r}) = \langle [{A_k}({\bf r}) - E[A]({\bf r}{)]^2}{ \rangle _k} = \langle [{A_k}({\bf r}) - E[A]({\bf r})] \times {[{A_k}({\bf r}) - E[A]({\bf r})]^*}{ \rangle _k}$, for both real-valued and complex-valued reconstructions. Higher-order PA fluctuation images may also be defined [20], but they are beyond the scope of this work. The variance image ${\sigma ^2}[A]({\bf r})$ can be expressed by use of the autocovariance function $C({{\bf r}_1},{{\bf r}_2})$:

B. Consequences on Visibility Artefacts

Equation (5) indicates that the fluctuation image involves a convolution with the square of the PSF. This is a well-known result from the Stochastic Optical Fluctuation Imaging (SOFI) approach initially proposed for super-resolution imaging relying on blinking fluorophores. The SOFI approach permits obtaining super-resolved images based on the fact that the PSF to the second power is sharper than the PSF itself [21]. Our group adapted the SOFI approach to super-resolution PA imaging, with fluctuations induced from either multiple-speckle illumination [17] or moving absorbers [20]. Although this convolution with the squared PSF resulting in the second-order fluctuation image has been known both for optical and PA imaging, its effect was analyzed only from the resolution perspective [17,20]. However, as we will now discuss, it also has a direct impact in the context of both the limited-view and resonant-transducer problems.

1. Illustration with Simulation Results

To illustrate our discussion, we first use two-dimensional (2D) numerical simulations showing how fluctuation imaging removes visibility artefacts. The numerical simulations were carried out with the 2D version of a FDTD freely available code developed by our group [22]. The simulations’ output consisted of a set of RF signals detected by a linear transducer array similar to that used in the experiments (central frequency 15 MHz, bandwidth 80%, $\text{N} = {128}$ elements, pitch 100 µm). The input of each simulation consisted of a source image defined on a Cartesian grid with a spatial step of 2.5 µm.

Three different structures $f({\bf r})$ were used to illustrate how PA fluctuation imaging removes artefacts of conventional PA imaging as illustrated in Fig. 1 (top row). Similar structures used in the phantom experiments are represented in Fig. 4. The simulation input for each realization of the random process consisted of a random source image ${g_k}({\bf r}) \times f({\bf r})$. The distribution ${g_k}({\bf r})$ was represented by an image of super-pixels (group of $q \times q$ pixels on the simulation grid) that was randomly assigned binary values 0 or 1 following a Bernouilli distribution with mean value $\eta$. The simulations permitted production of the full range of $\eta = \langle {{g_k}} \rangle = 0\% -100\%$. In this model, $\eta$ represents the surface fraction of the absorbed light and can be varied from 0% to 100%.

 

Fig. 1. Two-dimensional simulation results from complex-valued images. (${a_1} - {c_1}$) Typical realization for each of the three different structures $f({\bf r})$ with $\eta = 50\%$. For all three objects, the probe consists of a linear array ($\text{N} = {128}$ elements, 15 MHz central frequency, bandwidth 80%, pitch 100 µm) located 15 mm away above the objects. The corresponding experimental configurations are shown in Fig. 4. (${a_2} - {c_2}$) Modulus of the mean photoacoustic reconstruction obtained for each type of structure. (${a_3} - {c_3}$) Corresponding fluctuation images, free of visibility artefacts. Mean and fluctuation images are estimated from $\text{N} = {100}$ realizations.

Download Full Size | PDF

Figure 1 shows some example results obtained in our Monte Carlo simulations. To obtain these results, the surface fraction of absorbed energy was set to $\eta = 50\;\%$ and the size of each individual light absorbing patch was ${20}\;{\unicode{x00B5}\text{m}} \times {20}\;{\unicode{x00B5}\text{m}}$ ($16 \times 16$ super-pixel). Figure 1 shows the modulus images obtained from complex-valued signals. Figure 1 (column a) illustrates the limited-view problem: conventional PA imaging (middle row) only shows features parallel to the probe, while PA fluctuation imaging (bottom row) renders the full structure correctly. Figure 1 (column b) illustrates the limited-bandwidth problem: conventional PA imaging (middle row) only shows the upper and lower boundaries (high spatial frequencies) of the imaged object and does not reveal its inner content (low spatial frequencies), while PA fluctuation imaging (bottom row) correctly renders the full structure. Figure 1 (column c) illustrates both problems and confirms the outperformance of fluctuation imaging (bottom row). Such outperformance have first been demonstrated experimentally in our group several years ago in the specific context of multiple-speckle illumination [16], but without any clear theoretical explanation.

To illustrate the principle of the proposed approached, the simulations results presented in Fig. 1 were obtained without any additional sources of fluctuations (such as detection noise or laser pulse energy fluctuations). Additional simulation results presented in Supplement 1 (Section 4) confirm that our approach also works very efficiently under a noisy environment: pulse laser energy fluctuations may be compensated via singular value decomposition (SVD) filtering; the fluctuation of interest may be detected on top of a constant noise background provided that a sufficiently large number N of measurements is used to estimate the fluctuation image. Our simulation results with noise show that N may range in practice from a few units to a few hundreds (see Fig. S3 of Supplement 1), depending on the ratio between the fluctuation of interest and the noise level (FNR, see Section 3 of Supplement 1.)

2. Mechanism of Visibility Enhancement in Fluctuation Imaging

Based on the expression of the fluctuation image involving a convolution with the squared PSF, we now explain why PA fluctuation imaging removes visibility artefacts. The PSF corresponding to the images shown in Fig. 1 was computed in the same simulation setting by placing a point source in the center of the field of view (FOV) at the distance of $z = 15\;\text{mm} $ from the probe. We note ${h_r}({\bf r})$ the real-valued PSF and ${h_c}({\bf r})$ its complex-valued counterpart for the considered 1D linear transducer array (aligned along the horizontal direction as in Fig. 4). Figure 2 shows the PSF $h({\bf r})$ and its square modulus $|h{|^2}({\bf r})$ along with the corresponding spatial Fourier transforms.

 

Fig. 2. Two-dimensional simulation results: typical PSF $h({\bf r})$ and ${h^2}({\bf r})$ for limited-view photoacoustic imaging with resonant transducer, for both the real-valued PSF (${h_r}({\bf r})$) and the complex-valued PSF (${h_c}({\bf r})$). Top row: spatial domain. Bottom row: Fourier domain. All functions are normalized to 1. The PSF was derived from FDTD-computed signals, with the same 1D linear array as that used for Fig. 1, with a point source 15 mm away from the transducer.

Download Full Size | PDF

The bipolar nature of ${h_r}({\bf r})$ and ${h_c}({\bf r})$ for resonant transducers leads to low-spatial-frequency components of the objects being filtered out, which can be explained with Figs. 2(${a_2}$) and 2(${c_2}$) showing the PSF in the Fourier space. Moreover, the anisotropy of ${h_r}({\bf r})$ and ${h_c}({\bf r})$ leads to a selective detection of only spatial frequencies within the numerical aperture of the transducer. Previous analyses of $|h{|^2}({\bf r})$ in the context of super-resolution PA imaging used its extended support in the Fourier space as compared to that of $h({\bf r})$ to explain the obtained super resolution (see the Supplement 1 in Refs. [17,20]). However, it also follows from the support of $h_r^2({\bf r})$ and $|{h_c}{|^2}({\bf r})$ in the Fourier space [see Figs. 2(${c_2}$) and 2(${d_2}$)] that imaging with the squared modulus of the PSF preserves low spatial frequencies in all directions as opposed to $h({\bf r})$. As a consequence, the fluctuation image does not suffer from the limited-view and resonant-bandwidth artefacts observed in conventional imaging.

Figures 2(${c_1}$) and 2(${d_1}$) also illustrate how using the modulus of the complex-valued PSF eliminates the oscillatory behavior of ${h_r}({\bf r})$ [Figs. 2(${a_1}$) and 2(${b_1}$)]. Interestingly, the conclusions above drawn from the analysis in the Fourier space may also be obtained in the physical space by considering the difference between coherent and incoherent summation. Since $h({\bf r})$ is a bipolar function, its convolution with absorbing objects that are always positive may lead to destructive interference, which in practice occurs when dealing with large objects or those that emit PA waves away from the detectors. As opposed to $h({\bf r})$, ${h^2}({\bf r})$ is positive only and can therefore not lead to destructive interference. In other words, conventional imaging involves coherent summation of bipolar PA waves from positive sources, which may lead to destructive interference, whereas fluctuation imaging involves summation of positive fluctuations that can never interfere destructively. One can draw a parallel with optics where amplitudes sum up coherently while intensities sum up incoherently.

C. Amplitude of the Fluctuation Image

The mechanism discussed and illustrated above with Monte Carlo simulations and a simple random model is independent of the nature of fluctuations. These fluctuations may arise either from random heterogeneous illumination patterns (such as speckle patterns) or from random distributions of optical absorbers (such as red blood cells). However, the amplitude of the fluctuation image fundamentally depends on the statistical properties of the relevant random process (mean $\eta$, standard deviation $\sigma _g^2$, typical length scale ${D_g}$) as indicated by Eq. (5). To illustrate our discussion, we consider the fluctuation amplitude expected from structures much larger than the PSF. In this case, $\sqrt {{f^2}({\bf r}) * |h{|^2}({\bf r})} \simeq \sqrt {\int_{{\bf r}^\prime} |h{|^2}({\bf r}^\prime)\text{d}{\bf r}^\prime} = \Vert h{\Vert_2}$, and the fluctuation amplitude is given by

We first demonstrate by use of 2D numerical simulations that the fluctuation amplitude depends on $\eta$ and ${V_g} = D_g^2$ as predicted by our theory. To this end, we shall use the simple random model already introduced in Section 2.B.1. Then we provide more specific expressions for the fluctuation amplitude in the case of multiple-speckle illumination and in the case of random distributions of absorbing particles. We finally discuss why fluctuations from absorbing particles such as red blood cells are expected to be much larger than those from multiple-speckle illumination deep in tissue. We recall that Eqs. (5) and (6) are valid only for characteristic lengths ${D_g}$ much smaller than the shortest dimension of the PSF (about 50 µm as half a wavelength in the context of our work; see Fig. 2).

1. Validation with Simulations

To validate our theory with simulations, we consider a simple model of a random distribution of absorbed energy for which ${V_g}$ (and thus ${D_g}$) is independent of $\eta$. To do so, we use the model already introduced in Section 2.B.1, for which ${g_k}$ is defined on a Cartesian grid with pixel/voxel size $D$. Each pixel of the grid may acquire a random binary value 0 or 1 from a Bernoulli distribution with mean value of $\eta$. We recall that $\eta$ defined as such corresponds to the average surface/volume fraction $\eta = \langle {g_k}({\bf r}{) \rangle _k}$ of the corresponding random medium. The variance for the Bernoulli distribution is given by $\sigma _g^2 = \eta \times (1 - \eta)$. Moreover, it can be shown straightforwardly for this random medium that $C({{\bf r}^\prime} ,{{\bf r}^{\prime \prime}}) = \eta \times (1 - \eta)\prod\nolimits_{i = 1}^n {\Lambda _D}({x_{{i}}}^\prime - {x_{i}}^{\prime \prime})$, with $n $ being the space dimension ($n = {2}$ or $n = {3}$) and ${\Lambda _D}$ being the 1D unit triangular function with support $2 \times D$. From Eq. (2), one obtains that ${V_g} = {D^n}$ (i.e., ${D_g} = D$), independently of $\eta$. The expression for the fluctuation amplitude for this random distribution of absorbed energy becomes

To validate the predictions of Eq. (7), we ran 2D simulations ($n = 2$) for the disk structure (much larger than the PSF) with different values of $\eta$ and $D$. The fluctuation amplitude was estimated by spatially averaging the fluctuation image over a circular region of interest (see ROI in Fig. 3). As illustrated in Fig. 3, the dependence of the fluctuation amplitude on $\eta$ (for fixed $D = 10\;\unicode{x00B5}\text{m}$) and on the pixel size $D$ (for fixed $\eta = 50\;\%$) measured from the FDTD simulations ($N = 100$ random realizations for each set of parameters) exactly matches the one following from Eq. (7). In particular, the dependence on $\eta$ follows $\sqrt {\eta \times (1 - \eta)}$, and $\eta = 50\%$ is the value that maximizes the fluctuation amplitude. As expected, the size dependence diverges from the theory when $D$ approaches half a wavelength (50 µm).

 

Fig. 3. Two-dimensional FDTD simulation results. The fluctuation amplitude $\sigma [A]$ was estimated over a circular region of interest (ROI) for the disk configuration. The absorbed energy is distributed over square pixels of size $D$, with probability $\eta$. (a) Dependence of $\sigma [A]$ as a function of $D$ (for a fixed value $\eta = 50\%$). (b) Dependence of $\sigma [A]$ as a function of $\eta$ (fixed size $D = 20\;\unicode{x00B5}\text{m}$). The theoretical dependences are those given by Eq. (7) with $n = 2$.

Download Full Size | PDF

2. Multiple-Speckle Illumination

We consider here a homogeneously absorbing structure $f({\bf r})$ (constant absorption coefficient ${\mu _0}$), illuminated by random multiple-speckle illumination described by the fluence ${F_k}({\bf r}) = {F_0} \times {g_k}({\bf r})$. We further assume that the speckle is fully developed, so the statistical properties of the fluence distribution are well known: the intensity/fluence follows an exponential distribution, with ${\sigma _F} = \langle F \rangle $, where $\langle F \rangle = \eta {F_0}$ is the average intensity/fluence. ${D_g} = {D_s}$ is a measure of the speckle grain size, and it is independent of the average fluence $\eta {F_0}$. ${V_s} = D_s^n$ is the speckle grain area (i.e., $n = 2$) or volume (i.e., $n = 3$), depending on the relevant imaging configuration. The fluctuation amplitude in this case is thus given by

While it has long been acknowledged that detecting fluctuations from speckle illumination required large enough speckle grains, Eq. (8) provides for the first time a quantitative theoretical prediction of the expected fluctuation amplitude as a function of all relevant parameters, including the geometry of the imaged object and the PSF of the imaging system. For a fixed average fluence, Eq. (8) shows that the fluctuation amplitude is proportional to the square root of the area/volume of the speckle grain (depending on the relevant dimensionality).

3. Random Distribution of Absorbing Particles

We now consider the randomness that arises from the sample rather than from the illumination. In this case, the fluence is constant (${F_0}$), and the absorption distribution is described by ${\mu _k}({\bf r}) = {\mu _0} \times {g_k}({\bf r})$. More specifically, we consider randomness induced by random distributions of absorbing particles of a given characteristic volume ${V_p}$ (characteristic size ${D_p} = V_p^{1/n}$) inside the structure to image $f({\bf r})$, such as red blood cells flowing in blood vessels. In this case, ${g_k}({\bf r})$ is a binary function [${g_k}({\bf r}) = 1$ inside absorbing particles, ${g_k}({\bf r}) = 0$ outside], and $\eta = \langle {g_k}({\bf r}{) \rangle _k}$ simply represents the average volume fraction of absorbers. For blood, ${\mu _0}$ is the absorption coefficient of pure hemoglobin (i.e., of a single red blood cell), and $\eta {\mu _0}$ is the average absorption coefficient of the whole blood as a suspension of red blood cells. For ${g_k}({\bf r})$, which is a binary random variable with $\eta = \langle {g_k}({\bf r}{) \rangle _k}$ (Bernoulli distribution), the variance is $\sigma _g^2 = \eta \times (1 - \eta)$, regardless of the spatial structure of the random medium. In this case, the expression for the fluctuation image becomes

4. Multiple-Speckle Illumination versus Particle Fluctuations

As pointed out earlier, it is the ratio between the amplitude of the fluctuations of interest and the amplitude of other fluctuations (laser pulse energy fluctuations, thermal noise) that determines the ability to detect the relevant fluctuations experimentally. The theoretical expressions of the fluctuation images obtained for the case of multiple-speckle illumination and random distributions of absorbers provide a quantitative comparison between the two situations. From Eqs. (8) and (9), one gets

D. Conclusions from the Theory

In this part, we provided a unified theoretical framework to predict the amplitude of PA fluctuation images, whether fluctuations arise from random illuminations or from random distributions of absorbers. As a first result of our theoretical investigation, we explained why fluctuation PA imaging palliates the visibility artefacts emerging in conventional limited-view and resonant-bandwidth PA imaging. As a second and major result, we provided with Eq. (5) a general expression of the fluctuation amplitude as a function of statistical properties of the random process inducing fluctuations. We consequently demonstrated that fluctuations from red blood cells randomly distributed in blood vessels are expected to be 1 or 2 orders or magnitude larger than those expected from multiple-speckle illumination of blood vessels at depth in tissue. We recall that this conclusion holds because optical speckle grains at depth are much smaller (half the optical wavelength) than both the acoustic PSF [$h({\bf r})$] and the vessel $[f(\text{r})]$ as the main assumption of our theory. These theoretical results quantitatively explain why it has recently been possible to experimentally detect fluctuations from flowing red blood cells (in the context of super-resolution imaging [20]), whereas the use of speckle-induced fluctuations imaging remain limited to proof-of-concept experiments where speckle grains could be made large enough via free-space propagation. As a general conclusion, our theoretical results suggest that fluctuations from flowing red blood cells at physiological concentrations (around 50%) may be exploited to remove the visibility artefacts of conventional PA imaging. This has never been demonstrated experimentally so far, and is the objective of the second part of our work.

3. EXPERIMENTAL RESULTS

A. Phantom Experiments

As a first experimental demonstration of the proposed approach, we obtained cross-sectional images of the three different types of structures introduced in the simulation part. In each experiment, PA signal fluctuations were produced by a controlled flow of blood at physiological concentration.

1. Experimental Setup

The three types of experimental configurations are shown in Fig. 4. The first experiment consisted in imaging a C-shaped structure formed by a polycarbonate (PC) capillary (inner diameter $D = 100\;\unicode{x00B5}\text{m}$, wall thickness $w = 20\;\unicode{x00B5}\text{m}$, Paradigm Optics Inc., Vancouver, USA), while the second and the third experiments consisted in imaging a cross section of a glass tube (inner diameter $D = 1\;\text{mm} $, wall thickness $w = 100\;\unicode{x00B5}\text{m}$, Capillary Tube Supplies Ltd, UK) aligned within the imaging plane (secondexperiment) and perpendicularly to the imaging plane (thirdexperiment). Both samples and the probe were held in a water tank. For each experiment, blood flow was induced through the tubes and capillary via a syringe pump (KDS Legato 100, KD Scientific, Holliston, MA, USA). For all experiments, the blood velocity was approximately 1 cm/s.

 

Fig. 4. Imaging configurations used for the phantom experiments: (a) in-plane C-shaped capillary, (b) glass tube aligned in-plane, and (c) glass tube perpendicular to the imaging plane.

Download Full Size | PDF

The samples were illuminated in an oblique manner from the top of the probe, with $\tau = 5$ ns laser pulses [$\lambda = 680\;\text{nm} $, fluence $\approx 3\,\text{mJ}/{\text{cm}^2}$, pulse repetition frequency $(\text{PRF}) = {100}\;\text{Hz}$] by a frequency-doubled Nd:YAG laser (Spitlight DPSS 250, Innolas Laser GmbH, Krailling, Germany). At each laser shot, PA signals were recorded with a capacitive micromachined ultrasonic (CMUT) array (L22-8v, Verasonics, USA: 128 elements, pitch ${\sim}100\; \unicode{x00B5}\text{m}$, center frequency ${f_c} \sim 15\;\text{MHz} $, bandwidth 80%) connected to multichannel acquisition electronics (High Frequency Vantage 256, Verasonics, USA).

For each PA acquisition, the raw data was available as a RF frame containing the signals received by each transducer element of the CMUT array measured with a sampling frequency ${f_s} = 62.5\;\text{MHz}$. The total number of acquired RF frames was $M{1 = 1000}$ (C-shaped capillary), $M{2 = 1000}$ (glass tube, parallel orientation), and $M3 = 10,\!000$ (glass tube, perpendicular orientation). Given the 100 Hz PRF, the total acquisition times for each experiment were ${T_1} = {T_2} = 10\;\text{s}$ and ${T_3} = 100\;\text{s}$. A longer acquisition time was chosen for the third experiment to ensure that the number of independent red blood cell (RBC) configurations induced by the blood flow was similar in all three experiments. In contrast to the first and the second experiments, the blood flow in the third experiment was oriented perpendicular to the imaging plane. However, the characteristic PSF size in the direction perpendicular to the imaging ($x - z$) plane was about 10 times larger (of the order of 1 mm) than the dimension of the PSF within the imaging plane (of the order of 100 µm). With the chosen blood velocity of 1 cm/s, the RBC distribution over 100 µm (respectively 1 mm) distance is typically renewed every 10 ms (respectively 100 ms). Consequently, ${T_1}$, ${T_2}$, and ${T_3}$ each correspond to 1000 independent configurations during the measurement time.

2. Data Processing

For each of the three configurations, $N = 1000$ independent frames over the total acquisition time were used for data processing. For each experiment, $N = 1000$ complex-valued RF frames were computed from the $N = 1000$ real-valued frames via a Hilbert transform along the time axis. Complex-valued PA reconstructions ${A_k}(x,z)$ (with $\text{k} = {1}\cdots\text{N}$) were computed by applying standard delay-and-sum beamforming reconstruction. The mean reconstruction was computed as $E[A](x,z) = \frac{1}{N}\sum\nolimits_{k = 1}^M {A_k}(x,z)$, and the mean image was defined as $|E[A](x,z)|$. While the fluctuations of interest were those induced by the blood flow, other sources of fluctuations were present in the experiments, including laser pulse energy fluctuations and detection noise. Because laser pulse energy fluctuations identically affect all PA sources, they may be compensated by pulse-to-pulse energy monitoring. However, this requires additional hardware and signal processing.

Here we used spatiotemporal SVD as a key processing step to filter out laser pulse energy fluctuations (PEF). The SVD approach was first introduced in ultrasound imaging to discriminate tissue and blood motion [30] and was used by our group [20] later on to extract relevant fluctuations in fluctuation-based super-resolution PA imaging. In brief, SVD decomposes the initial data into a basis of spatiotemporal singular vectors. By choosing the singular vectors corresponding to relevant fluctuations, one can suppress signals with different spatiotemporal behaviour such as tissue motion, laser and electronics noise, etc. In our experiments, the first singular vectors (with the highest singular values) corresponded to laser PEF. In Supplement 1, we illustrate in both the simulation and experimental results how filtering out the first singular values indeed provides a very efficient way to remove the effect of parasitic laser PEF (see Supplement 1, Section 4.A, Fig. S2). We also illustrate the effect of SVD filtering as a function of the range of selected singular values (see Supplement 1, Section 4.C, Fig. S5). For all experimental images, the captions indicate the range of singular values that were selected (noted SVD cutoffs from now on).

For each set of $\text{N} = {1000}$ SVD-filtered PA images $A_k^{\text{SVD}}(x,z)$, the standard deviation image was computed through its variance as

3. Results

The experimental results are shown in Fig. 5. Row 1 shows the conventional photoacoustic image. Row 2 shows the fluctuation images that are obtained without the SVD step. Row 3 shows the fluctuation images obtained after SVD filtering. These results demonstrate that fluctuations induced by blood flow at physiological concentration can be exploited to suppress visibility artefacts of conventional PA imaging. Both limited-view and resonant-bandwidth artefacts are indeed well suppressed, in agreement with the simulation results shown in Fig. 1. The residual signal (clutter) observed at the bottom of images ${b_2}$ and ${c_2}$ is likely due to acoustic reverberation induced by the tube made of glass. This artefact is thus related to the sample rather than the proposed method. The slight inhomogeneity of the fluctuation amplitude within the C-shaped capillary (image ${a_2}$) may be caused by aberration induced by the capillary and/or an imperfect alignment of the structure in the imaging plane.

 

Fig. 5. Experimental results for the three configurations (a), (b), and (c) introduced in Figs. 1 and 4. Row 1: conventional photoacoustic images (mean). Row 2: fluctuation images (standard deviation) without SVD filtering. Row 3: fluctuation images (standard deviation) after SVD filtering. SVD cutoffs: a3 [31–61], b3 [11–381], and c3 [15–101]. Fluctuations were induced by blood at physiological concentrations flowing in the structures with a velocity around 1 cm/s. All images were obtained from processing $N = 1000$ complex-valued RF frames.

Download Full Size | PDF

As confirmed by simulations (see Fig. S2 of Supplement 1), the laser pulse energy fluctuations (PEF, measured with a powermeter around 3% relative rms value) completely overwhelms the fluctuation image if not compensated for. Our results illustrate that removing the first singular values efficiently filters out the effect of PEF and avoid demanding hardware-based monitoring and compensation of those fluctuations. Figure S5 of Supplement 1 further illustrates that once the low singular values are removed to filter out the PEF, the fluctuation images very marginally depend on the choice of SVD upper cutoffs.

For the results presented in Fig. 5, $\text{N} = {1000}$ frames were used to compute the fluctuation images. In Supplement 1, we present phantom experimental results obtained from different values of N (Fig. S3, bottom row), which show that depending on the desired image quality, N could be reduced from a few hundred to a few tens of images, increasing the temporal resolution by the same factor. In addition, the 1D profiles through the fluctuation image presented in Fig. S3 also clearly illustrate that the fluctuation of interest appears on top of a noise fluctuation background caused by detection noise. On the principle, provided that N is sufficiently large, any fluctuation of interest can be detected on top of the noise background, thanks to the additivity of the variance of independent noise sources. Importantly, we note that quantitative measurements (not considered here) require that the noise background can be estimated and substracted from the variance image. In all the results presented in this paper, the color map was set such that black corresponded to the background noise on the fluctuation image.

The dependence on N illustrates the trade-off between acquisition time/temporal resolution and quality of the image fluctuation, which in practice will depend on the level of fluctuation of interest as compared to noise. For the phantom experiments shown here, the peak-to-noise ratio in the signal domain was around 25–30, but the fluctuation of interest (from blood motion) were about half the noise rms value, which corresponded to a quite realistic situation. Thanks to both the beamforming step (which averages out RF noise by summation over the 128 channels) and the statistical estimation of the fluctuation along N realizations, the fluctuation of interest can in the end be detected in the image domain on top of the noise background, even if it is weaker.

B. In vivo Experiments

The proposed approach is further illustrated with a preliminary demonstration of 3D imaging of a live chicken embryo. The chorioallantoic membrane (CAM) of the chicken embryo is composed of small blood vessels whose role is to capture oxygen through the shell and to supply it to the developing embryo. Vessels are spread on a large surface and can have 3D orientation in the vicinity of the embryo, and they therefore provide a relevant model to assess the feasibility and performance of the proposed approach.

1. Experimental Setup and Methods

The experimental configuration is illustrated on Fig. 6. A custom 256-channel spherical transducer array (8 MHz, bandwidth 80%, $f = 35\;\text{mm}$, $f/D = 0.7$) was designed by our group and fabricated by Imasonics (Voray-sur-L’Oignon, France). This transducer was connected to the same Verasonics acquisition system used in our phantom experiments. An in-house coupling cone, filled with water and closed with a latex membrane, was used to couple the transducer to the sample. Fertilized eggs were obtained from a local farm and placed in an incubator for 6 days. After making a small hole through the shell using a scalpel, 1.5 mL of egg white was removed with a syringe, and the top part of the shell was cut with scissors. Warm phosphate-buffered saline (PBS) solution was poured into the shell to ensure acoustic contact with the latex membrane. The sample was held by a hammock made of food wrap. The hammock was attached to a bowl containing water maintained at 36°C using a hot plate. The sample was illuminated with $\lambda = 730\;\text{nm} $ light via a custom fiber bundle (Ceramoptec, Germany) going through the transducer, resulting in a fluence at the top surface of the sample of approximately $3\;{\text{mJ/cm}^{2}}$. $M = 768$ frames were acquired at 100 Hz and saved to the computer for offline data processing. The associated average power density, $300\;{\text{mW/cm}^{2}}$, was of the order of the ANSI limit of $230\;{\text{mW/cm}^{2}}$. Offline processing was performed exactly as in the 2D phantom experiments, with delay-and-sum beamforming, SVD filtering, and statistical computations.

 

Fig. 6. Schematic of the experimental setup for 3D in vivo imaging of a chicken embryo.

Download Full Size | PDF

2. Results

The top row of Fig. 7 shows the maximum intensity projections (MIPs) of the mean image in the three principal directions. These images present strong deterministic background features usually called clutter. In our experiment, this clutter is likely due to the limited number of transducers ($N = 256$) sparsely distributed over the ultrasound probe. Despite the relatively large detection aperture ($f/D = 0.7$), several structures oriented predominantly along the $z$ direction (upward and downward vessels) are missing or very weak in this conventional reconstruction because of the limited-view problem. In contrast, the MIPs of the 3D fluctuation images (bottom row of Fig. 7) reveal vessels with all possible orientations. This preliminary experiment also suggests that, in addition to removing the visibility artefacts, extracting the flow-induced fluctuations also allows removing the clutter that appears presumably due to the sparse nature of our probe (inducing a 3D PSF with grating lobes). For the results presented in Fig. 7, all the $\text{N} = {768}$ frames (measurement time 7.68 seconds) were used to compute the fluctuation images, with SVD $\text{cutoffs} = [{25 - 220}]$. In Supplement 1, we present results obtained from different values of N (Fig. S3) and also illustrate the effect of SVD filtering as a function of singular values cutoffs (Fig. 4(C)). As observed for the phantom experiments, filtering out the first singular values is essential to efficiently isolate the fluctuations from blood from those from the laser PEF. Figure S3 shows that $\text{N}=128$ (measurement time 1.28 s) still provides a relatively good image quality and illustrates the trade-off between acquisition time/temporal resolution and quality of the fluctuation image.

 

Fig. 7. In vivo images obtained experimentally with the chicken embryo CAM model. Mean intensity projections (MIP) of the 3D images reconstructions along the $x$, $y$, and $z$ directions. Top row: conventional PA imaging. Bottom row: PA fluctuation imaging. The mean and fluctuation images were obtained from $\text{N} = {768}$ frames, corresponding to a measurement time of 7.68s. SVD cutoffs: [25–220].

Download Full Size | PDF

The purpose of this preliminary experiment was to illustrate the performance of the method in a biologically relevant situation with realistic blood flow. A detailed description and investigation of the performance of our 3D imaging setup are being carried out and were out of the scope of the work presented here.

4. DISCUSSION AND CONCLUSION

In conclusion, we demonstrated both theoretically and experimentally that fluctuations induced by blood flow can be exploited to palliate visibility artefacts of conventional limited-view and resonant-bandwidth PA imaging. We first propose a theoretical framework that provides a quantitative expression of the fluctuation image as a function of the system PSF and statistical properties of the random medium at stake. Our theoretical framework applies to fluctuations induced either by fluctuating illumination patterns or by random distributions of absorbing particles (including flowing red blood cells). Our theoretical predictions suggested in particular that visibility problems could be overcome through exploiting the flow of blood at physiological concentration (overcoming visibility problems had been demonstrated previously only experimentally, and those experiments relied on multiple-speckle illuminations under unrealistic conditions).

We first provided two-dimensional images obtained with vessel-mimicking structures flowing with blood at physiological concentration. Finally, we presented a preliminary experiment demonstrating 3D fluctuation imaging in vivo of the vasculature of a chicken embryo CAM model. We emphasize that the proposed approach does not require additional hardware and can be implemented on conventional PA equipment since it relies on simple statistical processing.

The proposed technique has two main potential limitations inherent to its very principle, which are fundamentally related, as one may be overcome at the expense of the other: temporal resolution and sensitivity. As discussed in Section 3.A.3 and illustrated in Fig. S3 of Supplement 1, any fluctuation of interest can in principle be detected on top of any noise background (sensitivity), provided that the number N of acquisitions is made large enough, at the cost of temporal resolution.

What we illustrate in Supplement 1 (see Sections 3 and 4) is that the relevant number to quantify the level of fluctuation is the fluctuation-to-noise ratio (FNR), rather than the conventional signal-to-noise ratio on the PA signals. As a rule of thumb, as illustrated by our phantom experiments and additional simulation results provided in Section 4 of Supplement 1 (Fig. S2), for 128 channels, a FNR on the order of 1 requires $N \sim$ a few hundreds, while a FNR ${\sim}10$ requires $N \sim$ a few tens. Importantly, this estimate is valid for N uncorrelated frames, which was the case for our simulations and phantom experiments. For high enough laser pulse repetition rate, or equivalently for slow enough blood flow, there might be correlations between RF frames, in which case the robustness in the statistical estimation is affected. In particular, it should be kept in mind that samples with different dynamics (such as vascular networks with different vessel sizes and different blood velocities) may be reconstructed with a variable robustness. For instance, a robust reconstruction of blood microcapillaries may require an acquisition time longer than for large vessels, if the laser repetition period is short enough to resolve the dynamics in the microcapillaries.

In the case of our proof-of-concept in vivo experiment, the chicken embryo was mostly transparent, which is a relatively favorable situation as compared to imaging deep into tissue. The SNR ($\text{SNR} = \text{peak} - \text{to} - \text{noise}$ ratio in the signals domain) in this experiment was on the order of 50, about twice as high as that for the phantom experiments. However, the fluence value was approximately $3\;{\text{mJ/cm}^{2}}$, whereas this could be increased by a factor of 10 for deep-tissue experiments. In Section 6 of Supplement 1, we predict that the amplitude of the photoacoustic fluctuation image is typically 40 times (for our 15 MHz linear array) to 100 times (for our 8 MHz spherical 2D array) lower than the amplitude of the conventional photoacoustic image, when imaging 25 µm (linear array) and 50 µm blood vessels (spherical 2D array). Converted into the signal domain, this predicts that a FNR about 1 requires a SNR of about 20 for our linear array and a SNR of about 50 for our 2D spherical array, which corresponds to the cases of both our phantom experiments and in vivo experiments.

Although drawing firm conclusions relative to imaging deep in tissues must involve further in vivo investigations in small animals or humans, our results on the chicken embryo suggest that the technique has the potential to be applied for deep-tissue imaging, provided that a sufficient FNR can be reached, i.e., that the number N of acquisitions is sufficiently large. Our estimates above indicate that a SNR around a few tens should provide blood flow with a FNR about 1, which should provide very good images for N of about a few hundreds.

When compared to deep-learning-based approaches to solve the visibility problem, which can be implemented in a single-shot mode, our approach has a lower temporal resolution. However, it has the merit of extreme simplicity, and it can be implemented straightforwardly as opposed to deep-learning approaches, which require the availability of a training set and some quite specific know-how.

The proposed theoretical framework could also be potentially extended to investigate spatiotemporal cross-correlations of fluctuation images, which carry quantitative information on the flow of absorbing particles and thus may contribute to the field of PA Doppler measurements. Finally, although we focused here on the application to the visibility problem, our theoretical prediction of the second-order fluctuation image will also be useful to quantitatively exploit fluctuation-based super-resolution imaging.