A unified framework for photoacoustic fluctuation imaging. Application to visibility enhancement with fluctuations induced by blood flow

Photoacoustic fluctuation imaging relies on exploiting some randomness in photoacoustic generation to provide fluctuation-based enhanced images, as compared to conventional photoacoustic images. While this approach has been successfully demonstrated experimentally for super-resolution imaging (induced either from multiple speckle illumination or from the flow of random distributions of absorbing particles) and for visibility enhancement (with multiple speckle illumination only), to date there is no comprehensive theoretical description of fluctuation photoacoustic imaging. 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 predictions are then validated by Monte-Carlo finite-difference time-domain (FDTD) simulations of photoacoustic generation in a random medium. In particular, we explain why visibility artefacts are absent on 2nd-order fluctuation images. In the second part, we demonstrate experimentally that the flow of blood under physiological conditions can be exploited to produce photoacoustic fluctuation images free of the visibility artefacts of conventional photoacoustic imaging. Two-dimensional images of blood vessel phantoms flown with blood at physiological concentrations are obtained in vitro as a first proof-of-concept demonstration. Finally, the approach is applied in vivo to obtain 3D images of a chicken embryo.


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 can not be not 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 embryo.

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 discuss quantitatively why flowing red blood cells lead to PA fluctuations much larger than those produced by multiple speckle illumination deep in tissue.

Theoretical framework
PA imaging provides images of the absorbed energy density α that we describe in this work as: In the expression above, F 0 has the dimension of a light fluence, µ 0 has the dimension of an absorption coefficient; f (r) is a dimensionless binary function that describes the structure containing optical absorbers (for instance blood vasculature with red blood cells) : f (r) = 1 inside the imaged structure, f (r) = 0) outside; g k (r) is the k th realization of the random process inducing fluctuations: it may either describe a spatial distribution of light intensity or a spatial distribution of absorbers. For random illumination patterns, such as speckle patterns [13], F 0 × g k (r) (with g k being a continuous real-valued function) is the associated random fluence corresponding to each laser shot k, and the absorption coefficient µ 0 is considered homogenous in the absorbing structure f (r). For random distributions of absorbing particles, such as red blood cells, µ 0 × g k (r) is the associated heterogeneous absorption coefficient for each laser shot k while g k (r) is a binary function equal to 1 inside the particles, and 0 elsewhere. From the most general perspective, g k (r) × f (r) thus describes the spatial distribution of absorbed energy. We assume that g k (r) is a spatially stationary and isotropic random process. Its relevant statistical properties in the context of our work are the mean η =< g k (r) > k and the variance σ 2 g =< (g k (r) − η) 2 > k as first-order statistical properties and its autocovariance as a second-order statistical property. We note that C(0) = σ 2 g . We further assume that the normalized autocovariance C( r 1 − r 2 )/C(0) is a function with a finite volume V g and characteristic width D g (with D n g = V g , n=2 or n=3 depending on the relevant dimensionality): 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 mono-dispersed 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 (r) corresponding to laser shot k can be expressed as the convolution between α k (r) and the point spread function (PSF) h(r) of the imaging system: where Γ is the GrÃĳneisen parameter. The convolution in Eq. 3 implies a space-invariant PSF, which is a valid assumption within the small fields of view that we deal with in this work. We consider the most general framework where the PSF h(r) is either a real-valued function or a complex-valued function. A real-valued PSF corresponds to reconstruction based on real-valued RF signals, whereas a complex-valued PSF corresponds to reconstruction based on complex-valued signals. 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 is a widely used approach in medical ultrasound imaging.

Mean photoacoustic image
We consider the mean PA image defined by E[A](r) =< A k (r) > k , where < . > k denotes an ensemble average (average over laser shots in practice). As η =< g k (r) > k , E[A](r) can be straightforwardly expressed as: In the case of random distributions of absorbing particles, η =< g k (r) > k corresponds to the volume fraction of absorbers, and ηµ 0 is the average absorption coefficient of the imaged structure. In the case of random illumination patterns and a homogeneously absorbing structure, < g k (r)F 0 > k = ηF 0 corresponds to the average light fluence inside the imaged structure, and µ 0 is the absorption coefficient of the homogeneously absorbing imaged structure. Most generally, Eq. (4) corresponds to conventional PA reconstruction of a homogeneously absorbing structure described by f (r). Eq. (4) is independent of the real-valued or complex-valued nature of h(r). For real-valued reconstruction, E[A](r) directly represents the mean PA image, whereas for complex-valued reconstruction, the mean PA image is rather represented by the modulus |E[A](r)|.

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 (r)] * > k , for both real-valued and complex-valued reconstructions. Higher-order PA fluctuation images may also be defined [17], but are beyond the scope of this work. The variance image σ 2 [A](r) can be expressed by use of the autocovariance function C(r 1 , r 2 ): As the fundamental and major assumption of our work, we now assume that the characteristic size D g of the random distribution is much smaller than the shortest characteristic sizes of both h(r) and f (r). In the context of random distributions of red blood cells, D g is of the order of 10 µm.
In the context of speckle illumination, D g is the typical size of a speckle grain and is therefore not larger than 1 µm for visible or near-infrared light. Moreover, the shortest characteristic length scale of h(r) for resonant transducers is the wavelength λ at the central frequency. Our assumption is therefore fully verified for vessels at least several tens of µm in diameter and central frequencies up to several tens of MHz (acoustic wavelengths being not smaller than several tens of µm), which corresponds to the most of practical situations of interest. Under this assumption, C(r 1 , r 2 ) may be replaced in the integral by This leads to the following final expression for the PA fluctuation image, valid for both real-valued and complex-valued PSF: While it has already been shown by our group that the PA fluctuation image could be written as a convolution involving the square of the PSF [14,17], Eq. (5) for the first time explicitly links the amplitude of the fluctuation image to statistical properties of the random process at stake (through its mean η, its standard deviation σ 2 g and the volume of its normalized autocovariance V g ).

Consequences on visibility artefacts
Eq. (5 ) indicates that the fluctuation image involves a convolution with the square of the PSF. This is a well-known result from the 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 2 nd power is sharper than the PSF itself [18]. Our group adapted the SOFI approach to super-resolution PA imaging, with fluctuations induced from either multiple-speckle illumination [14] or moving absorbers [17]. Although this convolution with the squared PSF resulting in the 2 nd −order fluctuation image has been known both for optical and PA imaging, its effect was analyzed only from the resolution perspective [14,17]. However, as we will now discuss, it also has a direct impact in the context of both the limited-view and resonant-transducer problems.

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 finite-difference time-domain (FDTD) freely-available code developed by our group [19]. 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, N = 128 elements, pitch 110 µ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 (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 in vitro experiments are represented in Fig.4. The simulation input for each realization of the random process consisted of a random source image g k (r) × f (r). The distribution g k (r) was represented by an image of super-pixels (group of q × q pixels on the simulation grid), that was randomly assigned binary values 0 or 1 following a Bernouilli distribution with mean value η. The simulations permitted producing the full range of η =< g k > = 0 -100%. In this model, η represents the surface fraction of the absorbed light and can be varied from 0 to 100 %. Fig.1 shows some example results obtained in our Monte-Carlo simulations. To obtain these results, the surface fraction of absorbed energy was set super-pixel). Fig.1 shows the modulus images obtained from complex-valued signals. Fig.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. Fig.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. Fig.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 [13], but without any clear theoretical explanation.

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 mm from the probe. We note h r (r) the real-valued PSF and h c (r) its complex-valued counterpart for the considered 1D linear transducer array (aligned along the horizontal direction as in Fig.4). Fig.2 shows the PSF h(r) and its square modulus |h| 2 (r) along with the corresponding spatial Fourier transforms. The bipolar nature of h r (r) and h c (r) for resonant transducers leads to low spatial frequency components of the objects being filtered out, which can be explained with showing the PSF in the Fourier space. Moreover, the anisotropy of h r (r) and h c (r) leads to a selective detection of only spatial frequencies within the numerical aperture of the transducer. Previous analyses of |h| 2 (r) in the context of super-resolution PA imaging used its extended support in the Fourier space as compared to that of h(r), to explain the obtained super-resolution (see Supplementary Information in [14,17]). However, it also follows from the support of h 2 r (r) and |h c | 2 (r) in the Fourier space (see Figs.2.(b 2 ) and Fig.2.(d 2 )) that imaging with the squared modulus of the PSF preserves low spatial frequencies in all directions, as opposed to h(r). As a consequence, the fluctuation image does not suffer from the limited-view and resonant-bandwidth artefacts observed in conventional imaging. Fig.2.(c 1 ) and (d 1 ) also illustrate how using the modulus of the complex valued PSF eliminates the oscillatory behaviour of h r (r) (Fig.2.(a 1 ) and (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(r) is a bipolar function, its convolution with absorbing objects which 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(r), h 2 (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, that 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.

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 η, standard deviation σ 2 g , 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, f 2 (r) * |h| 2 (r) ∫ r |h| 2 (r )dr = ||h|| 2 , and the fluctuation amplitude is given by: while the average PA value is expected to be zero (E[A] ∝ f (r) * h(r) ∝ ∫ r h(r )dr = 0), as a consequence of the resonant bandwidth artefact (bipolar nature of h(r)). In practice, PA fluctuation amplitude must be larger than other fluctuations such as the the detection noise to be measurable. The fluctuation amplitude is in part determined by the sensitivity of the system (quantified through ||h|| 2 ), and by the statistical properties of the random distribution of absorbed energy (η, σ g and V g ). We first demonstrate by use of 2D numerical simulations that the fluctuation amplitude depends on η and V g = D 2 g as predicted by our theory. To this end, we shall use the simple random model already introduced in section 2.2.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 equations 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).

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 η. To do so, we use the model already introduced in section 2.2.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 η. We recall that η defined as such corresponds to the average surface/volume fraction η =< g k (r) > k of the corresponding random medium. The variance for the Bernoulli distribution is given by σ 2 g = η × (1 − η). Moreover, it can be shown straightforwardly for this random medium that C(r , r ) = η × (1 − η) n i=1 Λ D (x i − x i ), with n being the space dimension (n=2 or n=3) and Λ D being the 1D unit triangular function with support 2 × D. From Eq.2, one obtains that V g = D n (i.e. D g = D), independently of η. 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 η 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 η (for fixed D = 10 µm) and on the pixel size D (for fixed η = 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 η follows η × (1 − η) and η = 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).

Multiple speckle illumination
We consider here a homogeneously absorbing structure f (r) (constant absorption coefficient µ 0 ), illuminated by random multiple speckle illumination described by the fluence F k (r) = F 0 × g k (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 σ F =< F >, where < F >= ηF 0 is the average intensity/fluence. D g = D s is a measure of the speckle grain size, and is independent of the average fluence ηF 0 . V s = D n s is the speckle grain area (if n = 2) or volume if (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).

Random distribution of absorbing particles
We now consider 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 µ k (r) = µ 0 ×g k (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 1/n p ) inside the structure to image f (r), such as red blood cells flowing in blood vessels. In this case, g k (r) is a binary function (g k (r) = 1 inside absorbing particles, g k (r) = 0 outside) and η =< g k (r) > k simply represents the average volume fraction of absorbers. For blood, µ 0 is the absorption coefficient of pure hemoglobin (i.e. of a single red blood cell), and ηµ 0 is the average absorption coefficient of the whole blood as a suspension of red blood cells. For g k (r) which is a binary random variable with η =< g k (r) > k (Bernoulli distribution), the variance is σ 2 g = η × (1 − η), regardless of the spatial structure of the random medium. In this case, the expression for the fluctuation image becomes: The characteristic volume V g (η), defined by Eq.2, is of the order of V p , but its exact value generally depends on the volume fraction η, especially at high volume fraction where the positions of neighbors may become correlated. The dependence of σ[A](r) to the volume of the absorber (for a fixed η) is exactly the same as the dependence to the size of the speckle grain: the larger the particles, the larger the amplitude of the fluctuation image, within the limit that the particle size should be small compared to the PSF. For a given particle size and shape, the value of η that maximizes the fluctuation amplitude is given by the maximum of η(1 − η)V g (η), whose exact form depends on V g (η). In the simplest case where V g is independent of η, our theory predicts that the fluctuation amplitude is proportional to η × (1 − η), which is maximized for η = 50 %. It is out of the scope of this work to investigate the dependance of V g (η) for various types of media. Nonetheless, we demonstrate in the Supplementary information that V g (η) may be expressed as a function of V p through the so-called packing factor: The packing factor W(η) has been used in the field of ultrasound imaging to study the non linear dependence of echogeneicity as a function of η in particles media [2,3,6,22], including blood [3,[5][6][7]. By use of the expression above, the PA fluctuation may thus also be expressed as which interestingly has the same dependence on η × W(η) as that reported for the backscattering coefficient [6], and agrees with recent experimental results [8] (see Supplementary Inf.)

Multiple speckle illumination versus particles fluctuations
As pointed out earlier, it is the ratio between the amplitude of the fluctuations of interest and the amplitude of other fluctuations (laser shot-to-shot intensity 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: where it is assumed that the average fluence and the absorption coefficient are identical in both cases. Applied to red blood cells (V p ∼ 100 µm 3 , η ∼ 50 %, W(50%) ∼ 16 % [6] and Supp. Inf) and and speckle grains in the multiple scattering regime deep in tissue (D s ∼ λ 2 ∼ 0.25 µm in the visible range), the relation above (with n = 3) indicates that fluctuations arising from red blood cells as absorbing particles are at least one order of magnitude larger than fluctuations expected from multiple speckle illumination.

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 resonantbandwidth 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 one or two orders or magnitude larger than those expected from multiple speckle illumination of blood vessels at depth in tissue. 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 [17], 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.

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

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 µm, wall thickness w = 20 µ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 The total number of acquired RF frames was M 1 = 1,000 (C-shaped capillary), M2 = 1,000 (glass tube, parallel orientation) and M3 = 10,000 (glass tube, perpendicular orientation). Given the 100 Hz PRF, the total acquisition times for each experiment was T 1 = T 2 = 10 s and T 3 = 100 s. A longer acquisition time was chosen for the 3r d experiment to ensure that the number of independent RBC configurations induced by the blood flow was similar in all the three experiments. In contrast to the 1st and the 2nd experiments, the blood flow in the 3r d experiment was oriented perpendicular to the imaging plane. However, the characteristic PSF size in the direction perpendicular to the imaging (xz) plane was about ten 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, each of T 1 , T 2 and T 3 corresponds to 1000 independent configurations during the measurement time. Because shot-to-shot energy fluctuations affect identically all PA sources, they may be compensated by pulse-to-pulse energy monitoring. However, this requires additional hardware and signal processing. Here, we used spatio-temporal singular value decomposition (SVD) to filter out both shot-to-shot laser fluctuations and the detection noise. The SVD approach was first introduced in ultrasound imaging to discriminate tissue and blood motion [27] and was used by our group [17] 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 electronic noise, etc. In our experiments, the first singular vectors (with the highest singular values) corresponded to laser fluctuations, while noise was predominant in the singular vectors with the lowest singular values. The singular vectors corresponding to relevant fluctuations were determined empirically, and corresponded to singular values from 22 to 70, out of the total M=1000 values. For each set of M=1000 SVD-filtered PA images, A SVD k (x, z), the standard deviation image was computed as as the mean value is removed by SVD filtering along with the laser fluctuations. The standard deviation image is therefore given by the following expression: We note that the expression above is similar to what is used in ultrasound power doppler [27].

Results
The experimental results are shown in Fig. 5. These results demonstrate that fluctuations induced by a 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 the 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.

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, therefore, provide a relevant model to assess the feasibility and performance of the proposed approach.

Experimental setup and methods
The experimental configuration is illustrated on Fig.6. A custom 256-channel spherical transducer array (8M Hz, f = 35 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 in vitro 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 λ = 730nm 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 mJ/cm 2 . M = 768 frames were acquired at 100 Hz and saved to the computer for off-line data processing. The associated average power density, 300 mW/cm 2 , was of the order of the ANSI limit of 230 mW/cm 2 . Offline processing was performed exactly as in the 2D in vitro experiments, with standard delay-and-sum beamforming, SVD filtering and statistical computations. For the example shown here, SVD values between 11 and 60 out of the 768 values were conserved.

Results
The top row of Fig.7 shows the maximum intensity projections (MIP) of the mean image in the three principal directions. These images present strong deterministic background features usually 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. The purpose of this preliminary experiment was limited to illustrating the performance of the method in a biologically-relevant situation. A detailed description and investigation of the performance of our 3D imaging setup are being carried out and was out of the scope of the work presented here.

Conclusion
In conclusion, we demonstrated both theoretically and experimentally that fluctuations induced by a blood flow can be exploited to palliate visibility artefacts of conventional limited-view and resonant-bandwidth PA imaging. We first proposed a unified theoretical framework that provided 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 RBCs). 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 in vitro for vesselmimicking structures flown 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. It should, however, be kept in mind that the proposed approach requires many PA acquisitions in order to compute the standard deviation image, which is done at the cost of temporal resolution. Finally, the proposed theoretical framework could potentially be extended to investigate spatio-temporal cross-correlation of fluctuation images, which carry quantitative information on the flow of absorbing particles, and thus may contribute to the field of PA Doppler measurements.

Funding Information
This project has received funding from the European Research Council (ERC) under the European UnionâĂŹs Horizon 2020 research and innovation program (grant agreement No 681514-COHERENCE).

Disclosures
The authors declare no conflicts of interest.