# Large field-of-view non-invasive imaging through scattering layers using
fluctuating random illumination

Lei Zhu  Laboratoire Kastler Brossel, ENS–Université PSL, CNRS, Sorbonne
Université, College de France, 24 Rue Lhomond, F-75005 Paris, France.  School
of Physics and Optoelectronic Engineering, Xidian University, Xi’an 710071,
China  Fernando Soldevila  Laboratoire Kastler Brossel, ENS–Université PSL,
CNRS, Sorbonne Université, College de France, 24 Rue Lhomond, F-75005 Paris,
France.  Claudio Moretti  Laboratoire Kastler Brossel, ENS–Université PSL,
CNRS, Sorbonne Université, College de France, 24 Rue Lhomond, F-75005 Paris,
France.  Alexandra d’Arco  Laboratoire Kastler Brossel, ENS–Université PSL,
CNRS, Sorbonne Université, College de France, 24 Rue Lhomond, F-75005 Paris,
France.  Antoine Boniface  Laboratoire Kastler Brossel, ENS–Université PSL,
CNRS, Sorbonne Université, College de France, 24 Rue Lhomond, F-75005 Paris,
France.  Xiaopeng Shao  School of Physics and Optoelectronic Engineering,
Xidian University, Xi’an 710071, China  Hilton B. de Aguiar  Laboratoire
Kastler Brossel, ENS–Université PSL, CNRS, Sorbonne Université, College de
France, 24 Rue Lhomond, F-75005 Paris, France.  Sylvain Gigan  Laboratoire
Kastler Brossel, ENS–Université PSL, CNRS, Sorbonne Université, College de
France, 24 Rue Lhomond, F-75005 Paris, France.  Corresponding author:
<EMAIL_ADDRESS>
###### Abstract

Non-invasive optical imaging techniques are essential diagnostic tools in many
fields. Although various recent methods have been proposed to utilize and
control light in multiple scattering media, non-invasive optical imaging
through and inside scattering layers across a large field of view remains
elusive due to the physical limits set by the optical memory effect,
especially without wavefront shaping techniques. Here, we demonstrate an
approach that enables non-invasive fluorescence imaging behind scattering
layers with field-of-views extending well beyond the optical memory effect.
The method consists in demixing the speckle patterns emitted by a fluorescent
object under variable unknown random illumination, using matrix factorization
and a novel fingerprint-based reconstruction. Experimental validation shows
the efficiency and robustness of the method with various fluorescent samples,
covering a field of view up to three times the optical memory effect range.
Our non-invasive imaging technique is simple, neither requires a spatial light
modulator nor a guide star, and can be generalized to a wide range of
incoherent contrast mechanisms and illumination schemes.

##  INTRODUCTION

Non-invasive optical imaging has important applications in various fields
ranging from biotechnology [1, 2] to optical detection [3]. However,
inhomogeneous samples, such as biological tissues, scatter light, which
results in a complex speckle pattern on the detector [4, 5]. With increasing
depth, separating the low amount of ballistic light from the scattered light
becomes a big challenge [6, 7]. Over the years, many approaches have been put
forward to overcome this problem by exploiting or suppressing the scattered
light. With the development of spatial light modulators (SLMs), multiple ways
to control and manipulate scattered light have been developed [8, 9]. Several
techniques have been proposed to focus light by making use of feedback signals
to optimize the incident wavefront to recreate a focus that is then used for
raster-scanning microscopy [10, 11]. These techniques require access to both
sides of the scattering layer to optimize the wavefront, which strongly limits
their application in real-case scenarios. To overcome this, other strategies
have been proposed based on wavefront shaping and various feedback signals
such as fluorescence or ultrasound signals [11, 12, 13, 14]. However, these
approaches either require long acquisition times or are limited to small
fields of view (FoV). On the other hand, several techniques exploiting the
angular speckle correlations, known as the optical memory effect (ME) [15, 16,
17], have also been proposed for imaging objects hidden behind scattering
media [18, 19]. While these approaches are fast, their FoV is still limited by
the ME range.

Linear fluorescence is widely used in biology and biomedical sciences [20, 21,
22]. It enables imaging of cellular, subcellular, or molecular components, and
has the advantages of high spatial resolution, contrast, and speed. Recent
advances have allowed both focusing and imaging through scattering media using
fluorescent light. Even so, these methods either rely on the use of guide
stars [11], are limited to the ME range [23], or need to characterize the
scattering medium [24].

Here, we present a robust approach that allows to non-invasively image through
scattering layers at depths far beyond the ME range, using fluctuating random
speckle illumination simply generated with a rotating diffuser. Once excited,
each fluorescent emitter generates a unique speckle pattern on the detector,
which constitutes its fingerprint. Each image captured by the camera is an
incoherent sum of the fingerprints from the emitters, with different relative
weights due to the variable random illumination patterns and sample structure.
To retrieve each individual fingerprint, we capture a set of images while
randomly changing the illumination, and use a Non-negative Matrix
Factorization (NMF) algorithm to demix the set of acquired frames. After that,
the fingerprints are used to reconstruct the final image by exploring the
correlations between them. To validate the technique, we experimentally
demonstrate our non-invasive approach both on fluorescent beads and on
continuous fluorescent objects.

##  RESULTS

Figure 1: Schematic of the experimental setup and reconstruction principle.
(a) Schematic view of experimental setup. A coherent light source illuminates
a rotating diffuser in order to excite the fluorescent object through a
scattering medium with a random modulated speckle pattern. Once excited, the
emitted signal from the fluorescent objects is recorded with a camera.
$I_{fluo}$ is a series of $t$ fluorescent speckles corresponding to different
random speckle illuminations. The fingerprints can be recovered from
$I_{fluo}$ by using NMF. (b-e) Fingerprint-based reconstruction. (b) Pairwise
deconvolution (labelled as $\circledast^{-1}$) between all the possible pairs
of emitter fingerprints is performed. (c) The result of each deconvolution
provides the relative position between one emitter and its neighbors. (d) By
adding the resulting images for each emitter, it is possible to recover a
partial image of the object centered at that emitter (see Eq.4). (e) All the
partial images can be merged into the final reconstruction according to the
relative position between neighbouring emitters. Dashed circle indicates the
optical memory range. Scale bar sizes are
$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$. RD: rotating diffuser, DM:
dichroic mirror, OB: objective, Scat.: scattering medium, Fluo. Obj.:
fluorescent object, SF: spectral filter, TL: tube lens.

The experimental setup is depicted in Fig.1a. A rotating holographic diffuser
modulates the incident light coming from a laser by adding a random phase when
the light propagates through it. Then, the modulated light travels through the
scattering medium and generates a random unknown speckle which illuminates the
object. The excitation speckle induces a fluorescent response from the object,
which propagates back through the scattering medium, generating an incoherent
sum of speckles that is measured by the camera. Although the captured images
are low-contrast, random, and seemingly information-less, they contain all the
fingerprints from the independent emitters of the object, but with time-
varying weights. Furthermore, independent emitters within the ME range will
produce correlated but shifted fingerprints on the camera [15], while emitters
outside the ME range will produce totally uncorrelated fingerprints. For a
given speckle illumination, the captured image, $I_{\textsl{fluo}}$, can be
expressed as a linear superposition of those fingerprints with different
weights. Thus, the camera image is given by:

$I_{fluo}(r,t)=\sum^{P}_{k=1}w_{k}(r)h_{k}(t),$ (1)

where $I_{fluo}(r,t)$ corresponds to a low contrast speckle for the $t^{th}$
illumination, $r$ is the spatial coordinate, $w_{k}(r)$ represents the
fingerprint of the $k^{th}$ independent emitter of the object, $h_{k}(t)$
stands for the amount of excitation light at the $k^{th}$ emitter during the
$t^{th}$ illumination, and $P$ is the number of independent emitters. Given
enough different random illuminations, a collection of frames can be used to
retrieve each individual fingerprint, $w_{k}(r)$, by using an NMF algorithm,
that we will now explain in detail.

### Fingerprint demixing procedure

After randomly exciting the object with a variety of $t$ random speckles, a
series of camera images, $I_{fluo}(r,t)$, are collected. It is possible to
retrieve the fingerprints, $w_{k}$, corresponding to each independent emitter,
from the measurements $I_{fluo}(r,t)$ by finding the solution to the
minimization problem [25]:

$\min_{W>0,H>0}\Arrowvert I-WH\Arrowvert^{2},$ (2)

This minimization problem can be formulated as a low rank factorization, where
the matrix $I\in\mathbb{R}_{+}^{r\times t}$ contains all the $I_{fluo}(r,t)$,
can be approximated with two real positive matrices
$W\in\mathbb{R}_{+}^{r\times\rho}$ (the fingerprints) and
$H\in\mathbb{R}_{+}^{\rho\times t}$ (the temporal evolutions), where $r$ are
the pixels, $\rho$ is the estimated rank of $I$ and $t$ indicates the frames.
Since the collected images and the demixed fingerprints are positive, this
problem corresponds exactly to the family of NMF problems. The NMF framework
has been employed in demixing scenarios, both in structural imaging [24] and
functional imaging [26, 27, 28]. In our case, the estimated rank $\rho$
approximately corresponds to the number of independent emitters $P$ and it can
be estimated from the data by minimizing the root mean squared residual of NMF
as a function of the rank (see Supplementary 1).

###  Fingerprint-based reconstruction

After the demixing step, the fingerprints are retrieved. Due to the ME,
emitters close to each other will produce highly correlated fingerprints, with
a spatial shift that is directly determined by their relative position [19].
By exploring these correlations, a position map of the emitters can be
recovered, thus yielding an image of the object. Several approaches can be
used to explore the correlations and calculate the shifts between
fingerprints. Usually, this process is performed by doing a cross-correlation
between the fingerprints and locating the position of the maximum [24].
However, here we introduce a novel approach based on deconvolution, that we
denote as Fingerprint-based Reconstruction (FBR). Compared to the cross-
correlation procedure, we found that this approach allows to suppress noise
and strongly improves the quality of the reconstruction.

No matter if two emitters are within the same ME patch or not, one can perform
the pairwise deconvolution of the $i$-th emitter by the $k$-th emitter, which
can be written as:

$\underset{o_{i,k}}{\arg\min}\frac{\mu}{2}||w_{i}-o_{i,k}\circledast
w_{k}||^{2}_{2}+||o_{i,k}||_{TV},$ (3)

where $\mu$ is a regularization parameter, $\circledast$ denotes the
convolution operator, $||\cdot||_{2}$ indicates the $L_{2}$ norm and
$||\cdot||_{TV}$ represents the Total Variation (TV) norm, and the two
fingerprints are $w_{i}$ (considered as the ”image”) and $w_{k}$ (considered
as the ”point spread function” (PSF)). When the two emitters lay within one ME
range, the pairwise deconvolution yields a uniform image with a narrow delta-
like peak, which is located at a distance from the center given by the
relative position of the two emitters
($\vec{r}_{i,k}=\vec{r}_{i}-\vec{r}_{k}$). If the two emitters are located
beyond the ME range, the deconvolution yields noise.

For a given emitter $k$, it is possible to obtain, $O_{k}$, the partial image
of the object in the vicinity of the emitter, by simply adding the result of
all the pairwise deconvolutions related to that emitter, $o_{i,k}$. (see
Fig.1b-d).

$O_{k}=\sum^{\rho}_{i=1}o_{i,k},$ (4)

Even if the ensemble of emitters expands well beyond the ME range, the full
spatial distribution can be recovered if the different isoplanatic patches are
”connected” by emitters (see Fig.1d). For example, if emitters $i$ and $k$ are
beyond the ME range but emitter $j$ is between them, we can always calculate
the shift between them as $\vec{r}_{i,k}=\vec{r}_{i,j}+\vec{r}_{j,k}$. The
global reconstruction $O^{Global}$ can be obtained by composing all the
partial images, $O_{k}$, into one image, taking into account their relative
positions with respect to the first emitter, $\vec{r}_{k,1}$:

$O^{Global}=\sum^{\rho}_{k=1}O_{k}(\vec{r}-\vec{r}_{k,1}),$ (5)

Experimentally, we prove that our technique can be used to recover very sparse
objects by using 2D distributions of beads with a diameter of
$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$. As shown in Fig.2, our approach can
reconstruct the objects without constraints from the ME range. Non-sparse,
continuous objects are common in scattering biological samples, and this often
poses a difficult challenge for non-invasive imaging through scattering media
approaches. To demonstrate that our technique can also work with non-sparse
and continuous objects, we use fluorescence-stained pollen grains, whose
reconstructed images are shown in Fig.3. The size of fluorescence-stained
pollen grains in Fig.3 is smaller than the ME range and the process of
reconstruction used for them is the same as what is presented in Fig.1.

Figure 2:  Experimental results of imaging through a scattering medium with
fluorescent beads. (a)-(b) Fluorescent images of beads recorded without
scattering medium. (c)-(d) Reconstruction of the object using NMF+FBR
approach. The estimated rank of NMF is $\rho=26$ for (c) and $\rho=16$ for
(d). In both cases, $t=5120$ fluorescent speckle patterns are captured. The
exposure time of (c) and (d) is set to $15\text{\,}\mathrm{ms}$ and
$20\text{\,}\mathrm{ms}$, respectively. Dashed circles indicate the optical
memory effect range.

In both Fig.2 and Fig.3 we report the reconstruction of hidden beads
distributions and pollen seeds. However, the more complex the objects are, the
larger the required number of independent illuminations. Indeed, given more
speckle patterns, NMF is able to provide more reliable fingerprints, thus a
more reliable reconstruction (See supplementary Fig.10). We note that a few
hundred illuminations are sufficient to recover the object reasonably well in
our case. Importantly, our technique is not limited by the number of
independent illuminations that we can generate with the rotating diffuser, as
it is possible to produce a very large number of independent illumination with
different rotating diffusers by tunning their scattering angle [28]. As an
alternative to a rotating diffuser, we also propose a version of the setup
using a SLM, which allows reproducible pattern projections without practical
limitations on the number of patterns (see Supplementary 4).

Figure 3:  Experimental results of imaging through scattering media with
continuous objects. (a)-(b) Fluorescent images of different pollen seed
structures recorded without scattering medium. (c)-(d) Reconstruction of the
object with the NMF+FBR approach. The estimated rank for the NMF is $\rho=68$
for (c) and $\rho=85$ for (d), respectively. In both cases, $t=5120$
fluorescent speckle patterns are recorded with an exposure time of
$10\text{\,}\mathrm{ms}$. Dashed circle indicates the optical memory effect
range.

## Discussion

Here, we report on successfully recovering the hidden objects by exploiting
the correlation between the fingerprints. However, we believe, based on the
spectral ME or the 3D ME, that the technique could be used to recover multi
spectral or 3D objects. This is a significant advantage as it does not require
the scattering medium to be static. Another important future direction will be
to explore the ME of dynamic scattering media to recover the hidden object
through it. Interestingly, our technique, contrarily to the techniques based
on the optical transmission matrix, does not require focusing light on the
object, and relies only upon the random video frames generated from the random
illumination. Another important point is the simplicity of the technique,
which does not require calibration of the PSFs [29] of the imaging system, and
can be implemented without the need of expensive SLMs. Thus, our technique can
be easily implemented in various scattering media and imaging scenarios. It is
necessary to stress that our approach reconstructs both non-sparse and
continuous objects over previous autocorrelation approaches [18, 19].

In conclusion, we have shown a non-invasive technique to computationally
retrieve images of objects hidden behind a scattering medium from low-contrast
fluorescent speckles using random illumination. We have demonstrated that our
approach works with both sparse objects, even beyond the ME range, and
continuous objects. Importantly, the proposed approach neither relies upon
ballistic light nor uses wavefront shaping, and it is adaptable to various
scattering media and objects. Our technique is flexible, robust, and opens a
promising avenue towards deep fluorescence imaging in highly scattering media.
Finally, it can be generalized to a wide range of incoherent contrast
mechanisms and illumination schemes.

## METHODS

### Experimental setup

A continuous-wave laser ($\lambda=$532\text{\,}\mathrm{nm}$$, Coherent
Sapphire) is expanded and illuminates the rotating holographic diffuser
(Edmund, DG10). Then the modulated light is delivered onto the fluorescent
sample through a $200\text{\,}\mathrm{mm}$ lens (LA1708-A, Thorlabs) and
objective (Zeiss W “Plan-Apochromat” $\times 20$, NA $1.0$). After excitation,
the fluorescence is scattered by the medium and collected with a
$150\text{\,}\mathrm{mm}$ tube lens (L, AC254-150-A, Thorlabs), which is
employed to produce an image onto the detector, a sCMOS camera (Hamamatsu ORCA
Flash). Two dichroic filters (short pass $532\text{\,}\mathrm{nm}$, Thorlabs
and $533\text{\,}\mathrm{nm}$ notch MF525-39, Thorlabs) are used to block any
signal that does not come from the fluorescence emission. The fluorescent
objects, which are made of orange beads ($540\text{/}560\text{\,}\mathrm{nm}$,
Invitrogen FluoSpheres, size $1.0\text{\,}\mathrm{\SIUnitSymbolMicro m}$) or
pollen seeds (Carolina, Mixed Pollen Grains Slide, $w.m.$), are placed below
the scattering medium. The distance between the scattering medium and the
fluorescent objects is $0.2\text{\,}\mathrm{mm}$. A transmission pathway that
consists of a microscope objective (Olympus ”MPlan N” $\times 20$, NA $0.4$),
a $150\text{\,}\mathrm{mm}$ tube lens (L, AC$254-150-$A, Thorlabs), and CCD
camera (Allied Vision, Manta), is used as a passive control only. This control
part is used to correctly select the position of fluorescent object with a
white light source (Moritex, MHAB 150W) and it also allows us to align the
experimental setup. For the scattering medium, we either use a single
holographic diffuser (Newport,
10DKIT-C1,$10\text{\,}\mathrm{\SIUnitSymbolDegree}$) or a combination of two
(Newport, 10DKIT-C1,$10\text{\,}\mathrm{\SIUnitSymbolDegree}$ and
10DKIT-C1,$1\text{\,}\mathrm{\SIUnitSymbolDegree}$) in order to get different
memory effect ranges.

The exposure time has been set from $10\text{\,}\mathrm{ms}$ to
$20\text{\,}\mathrm{ms}$, depending on the scattering medium and the
fluorescent sample. Once captured, the speckle images, which contain few tens
of speckle grains, are cropped from the raw images. Then, a high-pass Fourier
Gaussian filter is employed to remove the background from the cropped images
and the processed data set is analyzed with the NMF algorithm to obtain the
fingerprint of each emitter. The experimental setup is shown in Supplementary
5.

In Fig.2a-b, the size of each cropped image is $70\times 72$ pixels, the
number of patterns $t$ is 5120, the scattering medium is a holographic
diffuser (Newport, 10DKIT-C1, $10$), and the exposure time is
$15\text{\,}\mathrm{ms}$. In Fig.2c-d, the size of each cropped image is
$64\times 66$ pixels, the number of patterns $t$ is 5120, the scattering
medium are two holographic diffusers (Newport, 10DKIT-C1, $10$ \+ $1$), and
the exposure time is $20\text{\,}\mathrm{ms}$. In Fig.3, the size of each
cropped image is $70\times 74$ pixels, the number of patterns $t$ is 5120, the
holographic diffuser (Newport, 10DKIT-C1, $10$) is used as the scattering
medium, and the exposure time is $10\text{\,}\mathrm{ms}$.

### NMF+fingerprint-based reconstruction algorithm

For the NMF, knowing the rank of the system is necessary. The rank $\rho$ is
estimated by looking at the root mean square residual $||I_{fluo}-WH||_{Fro}$
as a function of rank $\rho$ (detailed in Supplementary 1), where
$||\cdot||_{Fro}$ indicates the Frobenius norm, and minimizing it. For the
NMF, a random initialization is employed. The retrieved fingerprints are used
as the input data of the FBR.

### Acknowledgements

This research has been funded by the FET-Open (Dynamic-863203), European
Research Council ERC Consolidator Grant (SMARTIES-724473), China Scholarship
Council (CSC) (201906960055).

## Supplementary information

##  I Estimating the rank of the system

In order to operate the NMF algorithm, the rank, $\rho$, is a unique parameter
which needs to be pre-determined. Several methods can be used to estimate this
rank, such as looking at clusters using a k-means algorithm [26], or by
estimating the ’elbow’ of the loss function [24, 30]. In our case, we use a
simple procedure based on the study of the evolution of the root mean square
residual (RMSR) of the NMF for different ranks to estimate the true rank of
the system. Once we acquire a dataset, it is possible to demix it by setting
any desired rank $\rho$. For each $\rho$, we can estimate the quality of the
NMF solution by the RMSR: $||I_{fluo}-WH||_{Fro}$. As we change the rank, the
RMSR varies, reaching a minimum when the estimated rank equals the true rank
of the system. Also, when overestimating the rank, it is possible to have a
big variance on the result of the NMF for multiple realizations. Both these
behaviors can be used to delimit the true value of the rank of the system. As
an example of the procedure, we show both a simulation with different numbers
of beads and the same procedure with experimental datasets in Fig.4. As can be
seen in Fig.4a, the minimum value of RMSR is $\rho=5,10$, and $15$,
respectively, which is in excellent agreement with the ground truth number of
sources ($5,10$, and $15$). The rank estimation procedure is performed 12
times with different random initializations in order to see the variance of
the results.

Figure 4:  NMF rank estimation procedure. (a) Simulation results with $P$
fluorescent sources. The triangle symbol marks the minimum mean value of the
RMSR, the solid line stands for the mean value of the RMSR over 12 times with
the random initialization value of NMF, and the asterisk is the error bar that
indicates its standard deviation in both positive and negative direction. The
minimum mean RMSR value is achieved when the estimated rank $\rho$ equals the
true rank $P$ of the system. (b) Experimental results with a fluorescent
object which contains $P=11$ beads. (c) Estimating the rank of a fluorescent
object shown in (b), the solid line stands for the mean value of the RMSR over
12 repetitions, the triangle symbol marks the minimum mean value of the RMSR,
and the error bars indicate the standard deviation of the RMSR.

Furthermore, we investigate the effect of incorrectly estimating the ranks
$\rho$ of the system on the quality of our reconstruction. We experimentally
choose a fluorescent object which contains 10 beads. Then, we recover the
object using our approach with different rank estimations, ranging from 6 to
15. The structural similarity index metric (SSIM) [31] is used to quantify the
performance of reconstruction. As shown in Fig.5, our technique is capable of
achieving reconstruction with different ranks even if the proposed rank is not
fully accurate. However, underestimating the rank may lead to losing some
beads in the global reconstruction. On the other hand, overestimating the rank
on the NMF will produce mixed patterns that hinder the capability to obtain
the relative position between different emitters.

Figure 5:  The plot of SSIM coefficient between different reconstructions and
reference image. The reconstruction corresponding to $\rho=10$ is treated as a
reference image (the number of beads in the sample is $10$).

##  II Optical memory effect

With the purpose of proving that our technique can be applied in imaging
beyond the ME, different diffusers with different scattering properties are
tested. We estimate the ME range of each diffuser by exploring the correlation
between the speckle patterns generated while laterally moving an emitter (a
small fluorescent bead). As the source moves, the correlation between speckle
patterns gradually decreases, reaching a minimum correlation when the distance
is longer than the ME range. Measuring this distance allows estimating the ME
range for each diffuser. As reported in Fig.6, the ME range is approximately
$50\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and
$20\text{\,}\mathrm{\SIUnitSymbolMicro m}$ corresponding to diffuser $\\#1$
and diffuser $\\#2$, respectively. Our non-invasive technique not only can
retrieve objects within ME range as shown in Fig.2a of the manuscript, which
is possible to be recovered via autocorrelation approaches, but also can
reconstruct an extended object, which could not be retrieved with
autocorrelation approaches. As shown in Fig.2b of the manuscript, the size of
the reconstructed extended object is around
$$58\text{\,}\mathrm{\SIUnitSymbolMicro
m}$\times$22\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ with the diffuser $\\#2$.

Figure 6:  Plot of maximum value of cross-correlations between different
patterns corresponding to different fluorescent bead positions for diffuser
$\\#1$ (a) and $\\#2$ (b), where the colorcode indicates the different data
sets.

##  III Fingerprint-based reconstruction

As described in the main text, it is possible to reconstruct the whole object
as long as the ME patches have some overlaps. To present the details of FBR,
we experimentally choose a fluorescent object which contains 3 fluorescent
beads and display the various results from the pairwise deconvolution,
$o_{i,k}$, and the different partial images, $O_{k}$. As shown in Fig.7, the
$O_{k}$ of the emitter $k$ can be recovered by choosing $w_{k}$ as the PSF. By
looking at the maximum value of those $o_{i,k}$, the shift $\vec{r}_{i,k}$
between fingerprints $w_{i}$ and $w_{k}$ can be retrieved, as shown in Fig.1c
of the manuscript. In practice, the fingerprints coming from two emitters
which are beyond the ME range will not provide the relative position
information of their emitters (as they will be totally uncorrelated). Thus, it
is necessary to infer whether the fingerprints $w_{i}$ and $w_{k}$ are within
or beyond the ME range. In our method, we study
$\alpha=\frac{max\\{o_{i,k}\\}}{max\\{o_{k,k}\\}}$ as a function of relative
distance, where $max\\{o_{i,k}\\}$ stands for maximum value of $o_{i,k}$. A
given threshold, $\alpha_{tres}$ of $\alpha$ is introduced to evaluate it. For
example, if $\alpha$ is greater than $\alpha_{tres}$, $w_{i}$ and $w_{k}$
belong to the same ME range. Otherwise, they belong to different ME ranges,
and their relative position cannot be retrieved. As shown in Fig.8, we
experimentally investigate the maximum value of various results of pairwise
deconvolutions as a function of distance between emitters. We set
$\alpha_{tres}$ to $0.01$. The full workflow of technique is depicted in
Alg.1.

Figure 7:  Detail of FBR. (a) Ground truth taken without diffuser. (b) The
detail of pairwise deconvolution. The estimated rank $\rho$ of this data set
is 3 and the exposure time is $10\text{\,}\mathrm{ms}$. Figure 8:  (a)-(b)
Plot of $\alpha=\frac{max\\{o_{i,k}\\}}{max\\{o_{k,k}\\}}$ as a function of
relative distance between different fluorescent emitters for diffuser $\\#1$
and diffuser $\\#2$, where the colorcode represents the different data.

Input: Series of camera images, $I_{fluo}(r,t)$.

Output: Image of the object, $O^{Global}$.

Estimate the rank ($\rho$) of the system from $I_{fluo}(r,t)$ (see
Supplementary I).

Retrieve the spatial fingerprints ($w_{i}$) by using NMF with the estimated
rank.

for  _$k=1,...,\rho$_ do

Perform the pairwise deconvolution between $w_{k}$ and all the other
fingerprints ($w_{i\neq k}$).

Retrieve the relative position between emitter $k$ and its neighbours inside
the ME range.

Calculate the partial image of the object in the vicinity of the emitter
($O_{k}$) by adding the result of all the pairwise deconvolutions related to
that emitter ($o_{i,k}$).

end for

Merge all the partial images ($O_{k}$) into the final reconstruction
($O^{Global}$) using the relative position between emitters.

Algorithm 1 Image retrieval procedure

##  IV Producing random illumination using a SLM

Our technique retains the potential of employing a SLM for producing random
illumination. The SLM can produce a large number of independent speckle
illuminations. Experimentally, we replace the rotating diffuser with the SLM
and perform experiments on fluorescent point-like objects and continuous
volumetric objects by generating random illuminations with the SLM. The
reconstruction is presented in Fig.9. To prove that it is possible to recover
a more reliable image with more patterns, we show the reconstruction with the
different number of patterns in Fig.10. Note that the results as shown in
Fig.9 and Fig.10 are just tests to recover with the illumination patterns
produced by SLM that allows generating a higher diversity of patterns than our
current rotating diffuser, but it is not a limitation of our technique, as we
could use different rotating diffusers to get more independent illumination
patterns.

Figure 9:  Experimental results of imaging through scattering media by
producing random illumination with a SLM. (a,c) are ground truths of
fluorescent beads object and fluorescent continuous volumetric object. (c,d)
are reconstructions corresponding to (a,b), respectively. The estimated rank
$\rho$ of (c) and (d) is $14$ and $53$, respectively. $t=10240$ for (c) and
$t=5120$ for (d) fluorescent speckle patterns are recorded. The dataset (c) is
recorded with an exposure of $50\text{\,}\mathrm{ms}$ and the exposure time of
(d) is $20\text{\,}\mathrm{ms}$. Dashed circle indicates the optical memory
effect range. Figure 10:  Performance of reconstruction with different number
of speckle patterns by producing random illumination with a SLM. (a) Ground
truth of fluorescent continuous volumetric object taken without scattering
medium. (c)-(d) Reconstructions with different number of illumination
patterns. The rank of the NMF is estimated at $\rho=27$ for (b), at $\rho=31$
for (c), at $\rho=33$ for (d), at $\rho=36$ for (e), at $\rho=41$ for (f), and
at $\rho=53$ for (g). The dataset is recorded with an exposure of
$20\text{\,}\mathrm{ms}$.

##  V Experimental setup

Figure 11:  Experimental setup. The expanded $532\text{\,}\mathrm{nm}$ laser
beam illuminates the rotating diffuser and the modulated light is imaged on
the back focal plane of objective 1 (OB1). The object is placed in the focal
plane of OB1. Camera1 is located in the imaging plane of the fluorescent
imaging system that is employed to capture fluorescent speckle. The passive
controlling part is made of objective 2 (OB2), tube lens 2 (TL2), and camera
2. DM: dichroic mirror, SF: spectral filter, Scat.: scattering medium.

## References

  * [1] Zhao, M., Beauregard, D. A., Loizou, L., Davletov, B. & Brindle, K. M.  Non-invasive detection of apoptosis using magnetic resonance imaging and a targeted contrast agent.  _Nature Medicine_ 7, 1241–1244 (2001).
  * [2] Artzi, N. _et al._ In vivo and in vitro tracking of erosion in biodegradable materials using non-invasive fluorescence imaging.  _Nature Materials_ 10, 890–890 (2011).
  * [3] Kozloff, K. M. _et al._ Non-invasive optical detection of cathepsin K-mediated fluorescence reveals osteoclast activity in vitro and in vivo.  _Bone_ 44, 190–198 (2009).
  * [4] Goodman, J. W.  Some fundamental properties of speckle.  _Journal of the Optical Society of America_ 66, 1145 (1976).
  * [5] Bender, N., Yılmaz, H., Bromberg, Y. & Cao, H.  Customizing speckle intensity statistics.  _Optica_ 5, 595 (2018).
  * [6] Abramson, N.  Light-in-flight recording by holography.  _Optics Letters_ 3, 121 (1978).
  * [7] Huang, D. _et al._ Optical coherence tomography.  _Science_ 254, 1178–1181 (1991).
  * [8] Mosk, A. P., Lagendijk, A., Lerosey, G. & Fink, M.  Controlling waves in space and time for imaging and focusing in complex media.  _Nature Photonics_ 6, 283–292 (2012).
  * [9] Rotter, S. & Gigan, S.  Light fields in complex media: Mesoscopic scattering meets wave control.  _Reviews of Modern Physics_ 89, 015005 (2017).
  * [10] Vellekoop, I. M. & Mosk, A. P.  Focusing coherent light through opaque strongly scattering media.  _Optics Letters_ 32, 2309 (2007).
  * [11] Horstmeyer, R., Ruan, H. & Yang, C.  Guidestar-assisted wavefront-shaping methods for focusing light into biological tissue.  _Nature Photonics_ 9, 563–571 (2015).
  * [12] Katz, O., Ramaz, F., Gigan, S. & Fink, M.  Controlling light in complex media beyond the acoustic diffraction-limit using the acousto-optic transmission matrix.  _Nature Communications_ 10, 1–10 (2019).
  * [13] Popoff, S. M. _et al._ Measuring the transmission matrix in optics: An approach to the study and control of light propagation in disordered media.  _Physical Review Letters_ 104, 1–4 (2010).
  * [14] Hofer, M. & Brasselet, S.  Manipulating the transmission matrix of scattering media for nonlinear imaging beyond the memory effect.  _Optics Letters_ 44, 2137 (2019).
  * [15] Freund, I., Rosenbluh, M. & Feng, S.  Memory effects in propagation of optical waves through disordered media.  _Physical Review Letters_ 61, 2328–2331 (1988).
  * [16] Yılmaz, H. _et al._ Angular memory effect of transmission eigenchannels.  _Physical Review Letters_ 123, 203901 (2019).
  * [17] Osnabrugge, G., Horstmeyer, R., Papadopoulos, I. N., Judkewitz, B. & Vellekoop, I. M.  Generalized optical memory effect.  _Optica_ 4, 886\.
  * [18] Bertolotti, J. _et al._ Non-invasive imaging through opaque scattering layers.  _Nature_ 491, 232–234 (2012).
  * [19] Katz, O., Heidmann, P., Fink, M. & Gigan, S.  Non-invasive single-shot imaging through scattering layers and around corners via speckle correlations.  _Nature Photonics_ 8, 784–790 (2014).
  * [20] Ruan, H., Liu, Y., Xu, J., Huang, Y. & Yang, C.  Fluorescence imaging through dynamic scattering media with speckle-encoded ultrasound-modulated light correlation.  _Nature Photonics_ (2020).
  * [21] Lichtman, J. W. & Conchello, J. A.  Fluorescence microscopy.  _Nature Methods_ 2, 910–919 (2005).
  * [22] Mangeat, T. _et al._ Super-resolved live-cell imaging using random illumination microscopy.  _Cell Reports Methods_ 1, 100009 (2021).
  * [23] Hofer, M., Soeller, C., Brasselet, S. & Bertolotti, J.  Wide field fluorescence epi-microscopy behind a scattering medium enabled by speckle correlations.  _Optics Express_ 26, 9866 (2018).
  * [24] Boniface, A., Dong, J. & Gigan, S.  Non-invasive focusing and imaging in scattering media with a fluorescence-based transmission matrix.  _Nature Communications_ 11, 6154 (2020).
  * [25] Berry, M. W., Browne, M., Langville, A. N., Pauca, V. P. & Plemmons, R. J.  Algorithms and applications for approximate nonnegative matrix factorization.  _Computational Statistics and Data Analysis_ 52, 155–173 (2007).
  * [26] Moretti, C. & Gigan, S.  Readout of fluorescence functional signals through highly scattering tissue.  _Nature Photonics_ 14, 361–364 (2020).
  * [27] Pégard, N. C. _et al._ Compressive light-field microscopy for 3D neural activity recording.  _Optica_ 3, 517 (2016).
  * [28] Pnevmatikakis, E. A. _et al._ Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data.  _Neuron_ 89, 285 (2016).
  * [29] Antipa, N. _et al._ DiffuserCam: lensless single-exposure 3D imaging.  _Optica_ 5, 1 (2018).
  * [30] Hutchins, L. N., Murphy, S. M., Singh, P. & Graber, J. H.  Position-dependent motif characterization using non-negative matrix factorization.  _Bioinformatics_ 24, 2684–2690 (2008).
  * [31] Wang, Z., Bovik, A., Sheikh, H. & Simoncelli, E.  Image Quality Assessment: From Error Visibility to Structural Similarity.  _IEEE Transactions on Image Processing_ 13, 600–612 (2004).