# Photoacoustic fluctuation imaging: theory and application to blood flow
imaging.

Sergey Vilov  Univ. Grenoble Alpes, CNRS, LIPhy, 38000 Grenoble, France
Guillaume Godefroy  Univ. Grenoble Alpes, CNRS, LIPhy, 38000 Grenoble, France
Bastien Arnal  Univ. Grenoble Alpes, CNRS, LIPhy, 38000 Grenoble, France
Emmanuel Bossy  Univ. Grenoble Alpes, CNRS, LIPhy, 38000 Grenoble, France
Corresponding author<EMAIL_ADDRESS>
###### Abstract

Photoacoustic fluctuation imaging, which exploits randomness in photoacoustic
generation, provides enhanced images in terms of resolution and visibility, as
compared to conventional photoacoustic images. While a few experimental
demonstrations of photoacoustic fluctuation imaging have been reported, it has
to date not been described theoretically. In the first part of this work, we
propose a theory relevant to fluctuations induced either by random
illumination patterns or by random distributions of absorbing particles. The
theoretical predictions are validated by Monte Carlo finite-difference time-
domain simulations of photoacoustic generation in random particle media. We
provide a physical insight into why visibility artefacts are absent from
second-order fluctuation images. In the second part, we demonstrate
experimentally that harnessing randomness induced by the flow of red blood
cells produce photoacoustic fluctuation images free of visibility artefacts.
As a first proof of concept, we obtain two-dimensional images of blood vessel
phantoms. Photoacoustic fluctuation imaging is finally applied in vivo to
obtain 3D images of the vascularization in a chicken embryo.

##  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 can not 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, 7, 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 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. Sub-wavelength 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, 14, 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 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 signals fluctuations. In
particular, such fluctuations were deployed to obtain super-resolution in PA
imaging [17, 18, 19, 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, 18, 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 $2^{nd}$-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 a 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
flown 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 discuss quantitatively why flowing red blood cells lead to PA fluctuations
much larger than those produced by multiple speckle illumination deep in
tissue.

###  2.1 Analytical predictions

####  2.1.1 Theoretical framework

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

$\alpha_{k}(\mathbf{r})=\mu_{0}F_{0}[g_{k}(\mathbf{r})\times f(\mathbf{r})]$
(1)

In the expression above, $F_{0}$ has the dimension of a light fluence,
$\mu_{0}$ has the dimension of an absorption coefficient; $f(\mathbf{r})$ is a
dimensionless binary function that describes the structure containing optical
absorbers (for instance blood vasculature with red blood cells) :
$f(\mathbf{r})=1$ inside the imaged structure, $f(\mathbf{r})=0$) outside;
$g_{k}(\mathbf{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 [16], $F_{0}\times g_{k}(\mathbf{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 $\mu_{0}$
is considered homogenous in the absorbing structure $f(\mathbf{r})$. For
random distributions of absorbing particles, such as red blood cells,
$\mu_{0}\times g_{k}(\mathbf{r})$ is the associated heterogeneous absorption
coefficient for each laser shot $k$ while $g_{k}(\mathbf{r})$ is a binary
function equal to 1 inside the particles, and 0 elsewhere.

From the most general perspective, $g_{k}(\mathbf{r})\times f(\mathbf{r})$
thus describes the spatial distribution of absorbed energy. We assume that
$g_{k}(\mathbf{r})$ is a spatially stationary and isotropic random process.
Its relevant statistical properties in the context of our work are the mean
$\eta=<g_{k}(\mathbf{r})>_{k}$ and the variance
$\sigma_{g}^{2}=<(g_{k}(\mathbf{r})-\eta)^{2}>_{k}$ as first-order statistical
properties and its autocovariance

$C(\mathbf{r}_{1},\mathbf{r}_{2})=C(\left\|\mathbf{r}_{1}-\mathbf{r}_{2}\right\|)=<(g_{k}(\mathbf{r}_{1})-\eta)(g_{k}(\mathbf{r}_{2})-\eta)>_{k}$

as a second-order statistical property. We note that
$C(\boldsymbol{0})=\sigma_{g}^{2}$. We further assume that the normalized
autocovariance
$C(\left\|\mathbf{r}_{1}-\mathbf{r}_{2}\right\|)/C(\boldsymbol{0})$ is a
function with a finite volume $V_{g}$ and characteristic width $D_{g}$ (with
$D_{g}^{n}=V_{g}$, n=2 or n=3 depending on the relevant dimensionality):

$\int_{\mathbb{R}^{n}}C(\left\|\mathbf{r}\right\|)d\mathbf{r}=\sigma_{g}^{2}V_{g}=\sigma_{g}^{2}D_{g}^{n}$
(2)

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}(\mathbf{r})$ corresponding to
laser shot $k$ can be expressed as the convolution between
$\alpha_{k}(\mathbf{r})$ and the point spread function (PSF) $h(\mathbf{r})$
of the imaging system:

$A_{k}(\mathbf{r})=\Gamma\cdot\alpha_{k}(\mathbf{r})\ast h(\mathbf{r}),$ (3)

where $\Gamma$ 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(\mathbf{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. 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.1.2 Mean photoacoustic image

We consider the mean PA image defined by
$E[A](\mathbf{r})=<A_{k}(\mathbf{r})>_{k}$, where $<.>_{k}$ denotes an
ensemble average (average over laser shots in practice). As
$\eta=<g_{k}(\mathbf{r})>_{k}$, $E[A](\mathbf{r})$ can be straightforwardly
expressed as:

$E[A](\mathbf{r})=\Gamma F_{0}\mu_{0}\eta[f(\mathbf{r})\ast h(\mathbf{r})]$
(4)

In the case of random distributions of absorbing particles,
$\eta=<g_{k}(\mathbf{r})>_{k}$ corresponds to the volume fraction of
absorbers, and $\eta\mu_{0}$ is the average absorption coefficient of the
imaged structure. In the case of random illumination patterns and a
homogeneously absorbing structure, $<g_{k}(\mathbf{r})F_{0}>_{k}=\eta F_{0}$
corresponds to the average light fluence inside the imaged structure, and
$\mu_{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(\mathbf{r})$. Eq. (4) is independent of the real-valued or complex-valued
nature of $h(\mathbf{r})$. For real-valued reconstruction, $E[A](\mathbf{r})$
directly represents the mean PA image, whereas for complex-valued
reconstruction, the mean PA image is rather represented by the modulus
$\left|E[A](\mathbf{r})\right|$.

####  2.1.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](\mathbf{r})=<[A_{k}(\mathbf{r})-E[A](\mathbf{r})]^{2}>_{k}=<[A_{k}(\mathbf{r})-E[A](\mathbf{r})]\times[A_{k}(\mathbf{r})-E[A](\mathbf{r})]^{*}>_{k}$,
for both real-valued and complex-valued reconstructions. Higher-order PA
fluctuation images may also be defined [20], but are beyond the scope of this
work. The variance image $\sigma^{2}[A](\mathbf{r})$ can be expressed by use
of the autocovariance function $C(\mathbf{r}_{1},\mathbf{r}_{2})$:

$\begin{split}&\sigma^{2}[A](\mathbf{r})=\Gamma^{2}\mu_{0}^{2}F_{0}^{2}<\int_{\mathbf{r}_{1}}(g_{k}(\mathbf{r}_{1})-\eta)f(\mathbf{r}_{1})h(\mathbf{r}-\mathbf{r}_{1})\mathrm{d}\mathbf{r}_{1}\times\\\
&\int_{\mathbf{r}_{2}}(g_{k}(\mathbf{r}_{2})-\eta)f(\mathbf{r}_{2})h^{*}(\mathbf{r}-\mathbf{r}_{2})\mathrm{d}\mathbf{r}_{2}>_{k}=\\\
&=\Gamma^{2}\mu_{0}^{2}F_{0}^{2}\iint
C(\mathbf{r}_{1},\mathbf{r}_{2})f(\mathbf{r}_{1})f(\mathbf{r}_{2})h(\mathbf{r}-\mathbf{r}_{1})h^{*}(\mathbf{r}-\mathbf{r}_{2})d\mathbf{r}_{1}d\mathbf{r}_{2}\end{split}$

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(\mathbf{r})$ and $f(\mathbf{r})$.
In the context of random distributions of red blood cells, $D_{g}$ is of the
order of 10 $\mu$m. In the context of speckle illumination, $D_{g}$ is the
typical size of a speckle grain and is therefore not larger than 1 $\mu$m for
visible or near-infrared light. Moreover, the shortest characteristic length
scale of $h(\mathbf{r})$ for resonant transducers is the wavelength $\lambda$
at the central frequency. Our assumption is therefore fully verified for
vessels at least several tens of $\mu$m in diameter and central frequencies up
to several tens of MHz (acoustic wavelengths being not smaller than several
tens of $\mu$m), which corresponds to the most of practical situations of
interest. Under this assumption, $C(\mathbf{r}_{1},\mathbf{r}_{2})$ may be
replaced in the integral by
$\left[\int_{\mathbb{R}^{p}}C(\left\|\mathbf{r}^{\prime}\right\|)d\mathbf{r}^{\prime}\right]\times\delta(\mathbf{r}_{1}-\mathbf{r}_{2})=\sigma^{2}_{g}V_{g}\delta(\mathbf{r}_{1}-\mathbf{r}_{2})$:

$\begin{split}\sigma^{2}[A](\mathbf{r})&=\Gamma^{2}\mu_{0}^{2}F_{0}^{2}\times\\\
&\iint\sigma^{2}_{g}V_{g}\delta(\mathbf{r}_{1}-\mathbf{r}_{2})f(\mathbf{r}^{\prime})f(\mathbf{r}_{2})h(\mathbf{r}-\mathbf{r}_{1})h^{*}(\mathbf{r}-\mathbf{r}_{2})\mathrm{d}\mathbf{r}_{1}\mathrm{d}\mathbf{r}_{2}\end{split}$

This leads to the following final expression for the PA fluctuation image,
valid for both real-valued and complex-valued PSF:

$\sigma[A](\mathbf{r})=\Gamma\mu_{0}F_{0}\sigma_{g}\sqrt{V_{g}}\sqrt{f^{2}(\mathbf{r})\ast|h|^{2}(\mathbf{r})}$
(5)

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 [17, 20],
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 $\eta$, its standard deviation $\sigma^{2}_{g}$ and the volume of its
normalized autocovariance $V_{g}$). In particular, Eq. (5) remains valid when
the characteristic size of $f(\mathbf{r})$ is smaller than that of
$h(\mathbf{r})$, which is of interest for super-resolution imaging: as opposed
to our previous works on fluctuation-based super-resolution imaging [17, 20],
we here provide a quantitative prediction of the second-order super-resolution
photoacoustic fluctuation image, which is relevant for quantitative analysis
of such images. In this work, we however focus on the consequences of Eq. (5)
on the visibility problem.

###  2.2 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 Super-resolution
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 $2^{nd}$
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 $2^{nd}-$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.

####  2.2.1 Illustration with simulation results

Figure 1: Two-dimensional simulation results from complex valued images.
$a_{1}-c_{1}$: typical realization for each of the three different structures
$f(\mathbf{r})$ with $\eta=50\ \%$. For all three objects, the probe consists
of a linear array (N = 128 elements, 15 MHz central frequency, bandwidth
$80\%$, pitch 110 $\mu$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 N = 100
realizations.

.

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
[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\%$, N = 128 elements, pitch 110 $\mu$m). The
input of each simulation consisted of a source image defined on a cartesian
grid with a spatial step of 2.5 $\mu$m.

Three different structures $f(\mathbf{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 phantoms
experiments are represented in Fig.4. The simulation input for each
realization of the random process consisted of a random source image
$g_{k}(\mathbf{r})\times f(\mathbf{r})$. The distribution $g_{k}(\mathbf{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 producing the full range of $\eta=<g_{k}>$ = 0 - 100%. In this
model, $\eta$ 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 to
$\eta=50\ \%$ and the size of each individual light absorbing patch was $20\
\mu m\times 20\ \mu m$ ($16\times 16$ 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 [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 the Supplementary Information
(section 4) confirms that our approach also works very efficiently under noisy
environment: pulse laser energy fluctuations may be compensated for via 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
are used to estimate the fluctuations image. Our simulation results with noise
show that N may range in practice from a few units to a few hundred (see Fig.
S3 of the Supplementary Information), depending on the ratio between the
fluctuation of interest and the noise level (FNR, see section 3 of the
Supplementary Information.

####  2.2.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$ mm from the probe. We note
$h_{r}(\mathbf{r})$ the real-valued PSF and $h_{c}(\mathbf{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(\mathbf{r})$ and its square modulus $|h|^{2}(\mathbf{r})$ along with the
corresponding spatial Fourier transforms.

The bipolar nature of $h_{r}(\mathbf{r})$ and $h_{c}(\mathbf{r})$ for resonant
transducers leads to low spatial frequency components of the objects being
filtered out, which can be explained with Fig.2.($a_{2}$) and Fig.2.($c_{2}$)
showing the PSF in the Fourier space. Moreover, the anisotropy of
$h_{r}(\mathbf{r})$ and $h_{c}(\mathbf{r})$ leads to a selective detection of
only spatial frequencies within the numerical aperture of the transducer.
Previous analyses of $|h|^{2}(\mathbf{r})$ in the context of super-resolution
PA imaging used its extended support in the Fourier space as compared to that
of $h(\mathbf{r})$, to explain the obtained super-resolution (see
Supplementary Information in [17, 20]). However, it also follows from the
support of $h_{r}^{2}(\mathbf{r})$ and $|h_{c}|^{2}(\mathbf{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(\mathbf{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}(\mathbf{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(\mathbf{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(\mathbf{r})$, $h^{2}(\mathbf{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.

Figure 2: Two-dimensional simulation results: typical PSF $h(\mathbf{r})$ and
$h^{2}(\mathbf{r})$ for limited-view photoacoustic imaging with resonant
transducer, for both the real-valued PSF ($h_{r}(\mathbf{r})$) and the
complex-valued PSF ($h_{c}(\mathbf{r})$). Top row: spatial domain. Bottom row:
Fourier domain. All functions are normalized to one. 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.

###  2.3 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^{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,
$\sqrt{f^{2}(\mathbf{r})\ast|h|^{2}(\mathbf{r})}\simeq\sqrt{\int_{\mathbf{r}^{\prime}}|h|^{2}(\mathbf{r}^{\prime})\mathrm{d}\mathbf{r}^{\prime}}=||h||_{2}$,
and the fluctuation amplitude is given by:

$\sigma[A]=\Gamma\mu_{0}F_{0}\sigma_{g}\sqrt{V_{g}}||h||_{2}$ (6)

while the average PA value is expected to be zero ($E[A]\propto
f(\mathbf{r})\ast
h(\mathbf{r})\propto\sqrt{\int_{\mathbf{r}^{\prime}}h(\mathbf{r}^{\prime})\mathrm{d}\mathbf{r}^{\prime}}=0$),
as a consequence of the resonant bandwidth artefact (bipolar nature of
$h(\mathbf{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 ($\eta$, $\sigma_{g}$ and $V_{g}$).

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.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\
\mu\mathrm{m}$ as half a wavelength in the context of our work, see Fig. 2).

####  2.3.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.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 $\eta$. We recall that
$\eta$ defined as such corresponds to the average surface/volume fraction
$\eta=<g_{k}(\mathbf{r})>_{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(\mathbf{r}^{\prime},\mathbf{r}^{\prime\prime})=\eta\times(1-\eta)\prod_{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 :

$\sigma[A](\mathbf{r})=\Gamma\mu_{0}F_{0}\sqrt{\eta(1-\eta)D^{n}}||h||_{2}$
(7)

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\ \mu 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\ \mu\mathrm{m}$).

Figure 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\ \mu m$). The theoretical dependence are
those given by Eq.7 with $n=2$.

####  2.3.2 Multiple speckle illumination

We consider here a homogeneously absorbing structure $f(\mathbf{r})$ (constant
absorption coefficient $\mu_{0}$), illuminated by random multiple speckle
illumination described by the fluence $F_{k}(\mathbf{r})=F_{0}\times
g_{k}(\mathbf{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}=<F>$,
where $<F>=\eta 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
$\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:

$\sigma[A]_{\mathrm{speckle}}=\Gamma\mu_{0}<F>\sqrt{V_{s}}||h||_{2}$ (8)

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).

####  2.3.3 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 $\mu_{k}(\mathbf{r})=\mu_{0}\times
g_{k}(\mathbf{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(\mathbf{r})$, such as red blood cells flowing in blood vessels. In
this case, $g_{k}(\mathbf{r})$ is a binary function ($g_{k}(\mathbf{r})=1$
inside absorbing particles, $g_{k}(\mathbf{r})=0$ outside) and
$\eta=<g_{k}(\mathbf{r})>_{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}(\mathbf{r})$ which is a binary random variable with
$\eta=<g_{k}(\mathbf{r})>_{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:

$\sigma[A]_{\mathrm{particles}}=\Gamma\mu_{0}F_{0}\sqrt{\eta(1-\eta)V_{g}(\eta)}||h||_{2}$
(9)

The characteristic volume $V_{g}(\eta)$, defined by Eq.2, is of the order of
$V_{p}$, but its exact value generally depends on the volume fraction $\eta$,
especially at high volume fraction where the positions of neighbors may become
correlated. The dependence of $\sigma[A](\mathbf{r})$ to the volume of the
absorber (for a fixed $\eta$) 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
$\eta$ that maximizes the fluctuation amplitude is given by the maximum of
$\eta(1-\eta)V_{g}(\eta)$, whose exact form depends on $V_{g}(\eta)$. In the
simplest case where $V_{g}$ is independent of $\eta$, our theory predicts that
the fluctuation amplitude is proportional to $\sqrt{\eta\times(1-\eta)}$,
which is maximized for $\eta=50\ \%$. It is out of the scope of this work to
investigate the dependance of $V_{g}(\eta)$ for various types of media.
Nonetheless, we demonstrate in the Supplementary information that
$V_{g}(\eta)$ may be expressed as a function of $V_{p}$ through the so-called
packing factor:

$V_{g}(\eta)=\frac{W(\eta)}{1-\eta}V_{p}$ (10)

The packing factor $W(\eta)$ has been used in the field of ultrasound imaging
to study the non linear dependence of echogeneicity as a function of $\eta$ in
particles media [23, 24, 25, 26], including blood [24, 27, 26, 28]. By use of
the expression above, the PA fluctuation may thus also be expressed as

$\sigma[A]_{\mathrm{particles}}=\Gamma\mu_{0}F_{0}\sqrt{\eta
W(\eta)V_{p}}||h||_{2}$ (11)

which interestingly has the same dependence on $\eta\times W(\eta)$ as that
reported for the backscattering coefficient [26], and agrees with recent
experimental results [29] (see Supplementary Inf.)

####  2.3.4 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 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:

$\frac{\sigma[A]_{\mathrm{particles}}(\mathbf{r})}{\sigma[A]_{\mathrm{speckle}}(\mathbf{r})}=\sqrt{\frac{\eta
W(\eta)V_{p}}{V_{s}}}=\sqrt{\eta
W(\eta)\left(\frac{D_{p}}{D_{s}}\right)^{n}},$ (12)

where it is assumed that the average fluence and the absorption coefficient
are identical in both cases. Applied to red blood cells ($V_{p}\sim 100\ \mu
m^{3}$, $\eta\sim 50\ \%$, $W(50\%)\sim 16\ \%$ [26] and Supp. Inf) and and
speckle grains in the multiple scattering regime deep in tissue
($D_{s}\sim\frac{\lambda}{2}\sim 0.25\ \mu 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.

###  2.4 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 one or two 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(\mathbf{r})$) and the vessel size ($f(\mathbf{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

###  3.1 Phantoms 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.

####  3.1.1 Experimental setup

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

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 $\mu m$, wall thickness $w$ = 20 $\mu
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 mm, wall thickness $w$ = 100 $\mu m$, Capillary Tube Supplies
Ltd, UK) aligned within the imaging plane ($2^{nd}$ experiment) and
perpendicularly to the imaging plane ($3^{rd}$ experiment). 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.

The samples were illuminated in an oblique manner from the top of the probe,
with $\tau$ = 5 ns laser pulses ($\lambda$ = 680 nm, fluence $\approx$ 3
mJ/cm2, pulse repetition frequency (PRF) = 100 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\mu m$, center frequency $f_{c}\sim 15$ 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 radio-frequency (RF)
frame containing the signals received by each transducer element of the CMUT
array measured with a sampling frequency $f_{s}=62.5\ \mathrm{MHz}$. 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 $3rd$ experiment to ensure that the number
of independent red blood cells (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 $3rd$ 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\ \mu m$). With the chosen blood velocity of 1
cm/s, the RBC distribution over $100\ \mu 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.

####  3.1.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 k=1…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_{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 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) 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 the Supplementary Information, we illustrate on both 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
Supplementary information, section 4.A, Fig. S2). We also illustrate the
effect of SVD filtering as a function of the range of selected singular values
(see Supplementary information, section 4C, 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 N = 1000 SVD-filtered PA images, $A^{\mathrm{SVD}}_{k}(x,z)$,
the standard deviation image was computed through its variance as :

$\sigma[A](x,z)=\sqrt{\frac{1}{N-1}\sum_{k=1}^{N}\Bigl{\lvert}A^{\mathrm{SVD}}_{k}(x,z)-E[A^{\mathrm{SVD}}_{k}](x,z)\Bigr{\rvert}^{2}}$
(13)

We note that the expression above is similar to the one used in ultrasound
power Doppler (in which particular case the mean value of the SVD filtered
images is zero) [30].

Figure 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.

####  3.1.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 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.

As confirmed by simulations (see Fig. S2 of the supplementary material), 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. Fig. S5 of
the supplementary material further illustrates that once the low singular
values are removed to filter out the PEF, the fluctuation images very
marginally depends on the choice of SVD upper cutoffs

For the results presented on Fig.5, N = 1000 frames were used to compute the
fluctuation images. In the supplementary information, we present phantoms
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 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) requires that the noise background can be estimated and substracted from
the variance image. In all the results presented in this paper, the colormap
was set such that black corresponded to the background noise on the
fluctuation image.

The dependence on N illustrate the trade-off between acquisition time/temporal
resolution and quality of the fluctuation image, 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 realisations, the fluctuation of interest can in the end be detected in the
image domain on top of the noise background even if weaker.

###  3.2 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.

####  3.2.1 Experimental setup and methods

Figure 6: Schematic of the experimental setup for 3D in vivo imaging of
chicken embryo.

The experimental configuration is illustrated on Fig.6. A custom 256-channel
spherical transducer array ($8$ MHz, bandwidth $80\%$, $f=35\ \mathrm{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^{\circ}C$ using a hot
plate. The sample was illuminated with $\lambda$ = $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\
\mathrm{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\ \mathrm{mW/cm^{2}}$, was of the order of the ANSI limit of $230\
\mathrm{mW/cm^{2}}$. Offline processing was performed exactly as in the 2D
phantom experiments, with delay-and-sum beamforming, SVD filtering and
statistical computations.

####  3.2.2 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 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, 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 on Fig.7, all the N = 768 frames
(measurement time 7.68 seconds) were used to compute the fluctuation images,
with SVD cutoffs = [25-220]. In the supplementary information, 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. 4C). As
observed for the phantoms experiments, filtering out the first singular values
is essential to isolate efficiently the fluctuations from blood from those
from laser PEF. Fig. S3 shows that N = 128 (measurement time 1.28 seconds)
still provides a relatively good image quality, and illustrate the tradeoff
between acquisition time/temporal resolution and quality of the fluctuation
image.

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 was out of the scope of the work presented
here.

Figure 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 N =
768 frames, corresponding to a measurement time of 7.68s. SVD cutoffs:
[25-220].

##  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 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.

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 on Fig. S3 of the Supplementary Information,
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 the temporal resolution.

What we illustrate in the supplementary material (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 the
Supplementary Information (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 phantoms experiments. For laser pulse
repetition rate high enough, or equivalently for blood flow slow enough, 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, an robust reconstruction of blood
microcapillaries may require a longer acquisition time that for larger
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 (SNR = peak to noise ratio in the signals
domain) in this experiment was on the order of 50, about twice higher than for
the phantom experiments. However, the fluence value was approximately $3\
\mathrm{mJ/cm^{2}}$, whereas this could be increased by a factor ten for deep
tissue experiments. In section 6 of the Supplementary material, 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 blood vessels $25\ \mu$m (linear array) and $50\ \mu$m (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 tissue must
involve further in vivo investigations in small animal 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 fluctuation to
noise ratio (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 a FNR about 1, which should provide very good images
for N 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. It however has the merit of extreme simplicity, and
can be implemented straightforwardly, as opposed to deep-learning approaches
which requires the availability of a training set, and some quite specific
know-how.

The proposed theoretical framework could also be potentially extended to
investigate spatio-temporal 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 be also be useful to
quantitatively exploit fluctuation-based super-resolution imaging.

##  5 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).

##  6 Disclosures

The authors declare no conflicts of interest.

## References

  * [1] P. Beard, “Biomedical photoacoustic imaging,” Interface focus 1, 602–631 (2011).
  * [2] D. Yang, D. Xing, S. Yang, and L. Xiang, “Fast full-view photoacoustic imaging by combined scanning with a linear transducer array,” Optics express 15, 15566–15575 (2007).
  * [3] R. A. Kruger, W. L. Kiser Jr, D. R. Reinecke, and G. A. Kruger, “Thermoacoustic computed tomography using a conventional linear transducer array,” Medical physics 30, 856–860 (2003).
  * [4] B. Huang, J. Xia, K. I. Maslov, and L. V. Wang, “Improving limited-view photoacoustic tomography with an acoustic reflector,” Journal of biomedical optics 18, 110505 (2013).
  * [5] G. Li, J. Xia, K. Wang, K. Maslov, M. A. Anastasio, and L. V. Wang, “Tripling the detection view of high-frequency linear-array-based photoacoustic computed tomography by using two planar acoustic reflectors,” Quantitative Imaging in Medicine and Surgery 5 (2014).
  * [6] B. T. Cox, S. R. Arridge, and P. C. Beard, “Photoacoustic tomography with a limited-aperture planar sensor and a reverberant cavity,” Inverse Problems 23, S95–S112 (2007).
  * [7] R. Ellwood, E. Zhang, P. Beard, and B. Cox, “Photoacoustic imaging using acoustic reflectors to enhance planar arrays,” Journal of Biomedical optics 19, 126012 (2014).
  * [8] B. Cox and P. Beard, “Photoacoustic tomography with a single detector in a reverberant cavity,” The Journal of the Acoustical Society of America 125, 1426–1436 (2009).
  * [9] X. L. Dean-Ben and D. Razansky, “Localization optoacoustic tomography,” Light: Science & Applications 7, 18004 (2018).
  * [10] X. L. Deán-Ben, L. Ding, and D. Razansky, “Dynamic particle enhancement in limited-view optoacoustic tomography,” Opt. Lett. 42, 827–830 (2017).
  * [11] P. Zhang, L. Li, L. Lin, J. Shi, and L. V. Wang, “In vivo superresolution photoacoustic computed tomography by localization of single dyed droplets,” Light: Science & Applications 8, 1–9 (2019).
  * [12] L. Wang, G. Li, J. Xia, and L. V. Wang, “Ultrasonic-heating-encoded photoacoustic tomography with virtually augmented detection view,” Optica 2, 307–312 (2015).
  * [13] M. W. Kim, G.-S. Jeng, I. Pelivanov, and M. O’Donnell, “Deep-learning image reconstruction for real-time photoacoustic system,” IEEE Transactions on Medical Imaging (2020).
  * [14] N. Davoudi, X. L. Deán-Ben, and D. Razansky, “Deep learning optoacoustic tomography with sparse data,” Nature Machine Intelligence 1, 453–460 (2019).
  * [15] G. Godefroy, B. Arnal, and E. Bossy, “Compensating for visibility artefacts in photoacoustic imaging with a deep learning approach providing prediction uncertainties,” Photoacoustics p. 100218 (2020).
  * [16] J. Gateau, T. Chaigne, O. Katz, S. Gigan, and E. Bossy, “Improving visibility in photoacoustic imaging using dynamic speckle illumination,” Optics letters 38, 5188–5191 (2013).
  * [17] 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).
  * [18] 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,” Optics express 25, 4875–4886 (2017).
  * [19] 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).
  * [20] T. Chaigne, B. Arnal, S. Vilov, E. Bossy, and O. Katz, “Super-resolution photoacoustic imaging via flow-induced absorption fluctuations,” Optica 4, 1397–1404 (2017).
  * [21] T. Dertinger, R. Colyer, G. Iyer, S. Weiss, and J. Enderlein, “Fast, background-free, 3d super-resolution optical fluctuation imaging (sofi),” Proceedings of the National Academy of Sciences 106, 22287–22292 (2009).
  * [22] E. Bossy, “Simsonic, a fdtd ultrasound simulation freeware,” .
  * [23] V. Twersky, “Acoustic bulk parameters in distributions of pair-correlated scatterers,” The Journal of the Acoustical Society of America 64, 1710–1719 (1978).
  * [24] K. K. Shung, “On the ultrasound scattering from blood as a function of hematocrit,” IEEE Transactions on Sonics and Ultrasonics 29, 327–330 (1982).
  * [25] V. Twersky, “Low-frequency scattering by correlated distributions of randomly oriented particles,” The Journal of the Acoustical Society of America 81, 1609–1618 (1987).
  * [26] P. A. Bascom and R. S. Cobbold, “On a fractal packing approach for understanding ultrasonic backscattering from blood,” The Journal of the Acoustical Society of America 98, 3040–3049 (1995).
  * [27] N. Berger, R. Lucas, and V. Twersky, “Polydisperse scattering theory and comparisons with data for red blood cells,” The Journal of the Acoustical Society of America 89, 1394–1401 (1991).
  * [28] D. Savery and G. Cloutier, “A point process approach to assess the frequency dependence of ultrasound backscattering by aggregating red blood cells,” The Journal of the Acoustical Society of America 110, 3252–3262 (2001).
  * [29] T. M. Bücking, “Acoustic resolution photoacoustic flowmetry,” Ph.D. thesis, UCL (University College London) (2019).
  * [30] C. Demené, T. Deffieux, M. Pernot, B.-F. Osmanski, V. Biran, J.-L. Gennisson, L.-A. Sieu, A. Bergel, S. Franqui, and J.-M. Correas, “Spatiotemporal clutter filtering of ultrafast ultrasound data highly increases doppler and fultrasound sensitivity,” IEEE transactions on medical imaging 34, 2271–2285 (2015).

sample