## INTRODUCTION

For decades, optical microscopy has been a vital tool in biomedical research to observe live specimens with a submicrometer resolution and with minimal invasiveness. Yet, imaging conditions required for these exquisite performances are rarely gathered. For instance, both the resolution and the contrast drop as the imaging depth increases inside a biological tissue. This observation is a consequence of the spatial variations of the specimen’s refractive index that distort the wavefront of both the incoming and outgoing light. When these variations exhibit low spatial frequencies, we use the term aberrations while scattering describes the effect of the higher spatial variations. Both these effects limit the use of conventional microscopy to shallow depths or to semitransparent specimens. Imaging deeper requires to simultaneously compensate for these detrimental phenomena.

To mitigate the aberrations induced by the specimen, the concept of adaptive optics (AO) has been adapted to microscopy from astronomy where it was developed decades ago (*1*, *2*). Astronomers faced the same impediment as fluctuations in the atmosphere severely distort the wavefront of the light coming from stars and prevent to obtain a diffraction-limited stellar image. Astronomers then proposed to measure these distortions using a wavefront sensor and to counterbalance it with a dynamic programmable element such as deformable mirrors. Following this concept and the development of deformable mirrors with increasing number of elements, AO already demonstrated its benefits in various imaging techniques such as digital holography (*3*, *4*), confocal microscopy (*5*, *6*), two-photon microscopy (*7*–*10*), or optical coherence tomography (OCT) (*11*, *12*). Unfortunately, AO methods usually require a guide star or are based on an image sharpness metric. In addition, they are limited to a small region called the isoplanatic patch (IP), the area over which the aberrations can be considered as spatially invariant. Therefore, there is a need to extend the field of view of AO methods by tackling the case of multiple IPs. This issue is particularly decisive for deep imaging where IP size becomes extremely tiny: <10 μm beyond a depth of 1 mm (*13*). Note that multiconjugate AO can deal with multiple IPs, but this is at the price of a much more complex optical setup (*14*–*16*).

On the other hand, multiple scattering was long thought to be too complex to be compensated. For deep imaging, a gating mechanism is generally used to reject the multiply scattered photons and capture only the ballistic light. This gating can be spatial (*17*) as in confocal microscopy or temporal (*18*) as in OCT, but they are still depth limited as they rely on the exponentially attenuated ballistic light. In a pioneering experiment, Vellekoop and Mosk (*19*) demonstrated in 2007 the possibility to restore a diffraction-limited spot through a scattering medium by properly shaping the incoming light. Subsequently, a matrix approach of light propagation through complex media was developed (*20*). Relying on the measurement of the Green’s functions between each pixel of a spatial light modulator (SLM) and of a charge-coupled device (CCD) camera across a scattering medium, the experimental access to the transmission matrix allows taking advantage of multiple scattering for optimal light focusing (*20*) and communication across a diffusive layer (*21*, *22*) or a multimode fiber (*23*, *24*). However, a transmission configuration is not adapted to noninvasive and/or in vivo imaging of biological media. An epidetection geometry should thus be considered (*25*). During the last few years, the reflection matrix R had been investigated to perform selective focusing/detection (*26*, *27*) or energy delivery (*28*, *29*) through strongly scattering media. With regard to the specific purpose of imaging, the matrix approach has been recently used to implement AO tools in postprocessing. The single scattering component of the reflected wave field through biological tissues has been enhanced in depth by compensating for high-order aberrations (*30*, *31*).

Here, we propose to go beyond a matrix approach of AO by introducing a novel operator: the so-called distortion matrix D. Unlike previous works that investigated R either in the focal plane (*27*) or the pupil plane (*26*, *30*, *31*), we here consider the medium response between those dual bases (*32*, *33*). Unlike R, the D-matrix does not consider the reflected wave field as a building block but its deviation from an ideal wavefront that would be obtained in the absence of aberrations and without multiple scattering. This operation may seem trivial but it markedly highlights the input/output correlations of the wave field. While the canonical reflection matrix exhibits a random feature in a turbid medium, the distortion matrix displays strong field-field correlations over each IP. Thanks to this new operator, some relevant results of information theory can thus be fruitfully applied to optical imaging. A singular value decomposition (SVD) of D allows a partition of the field of illumination (FOI) into orthogonal isoplanatic modes (IMs) and extracts the associated wavefront distortion in the pupil plane. The Shannon entropy ℋ of the singular values allows one to define the effective rank of the imaging problem. A combination of the ℋ first eigenstates yields an image of the focal plane with an excellent contrast and a diffraction-limited resolution as if the medium ahead was made perfectly transparent.

Several experiments with an increasing order of complexity are presented to demonstrate the benefits of the D-matrix for optical imaging in turbid media. For the sake of simplicity, the first experiment involves the imaging of a single IP through a thick layer of biological tissues. This configuration allows us to lay down the D-matrix concept and highlight the physics behind it. Then, a second proof-of-concept experiment considers a thin but strong aberrating layer introduced between the microscope objective (MO) and a resolution target. This imaging configuration involves a spatially varying aberration across the FOI (i.e., several IPs). Last, we describe an imaging experiment through a turbid nonhuman primate cornea that induces high-order aberrations (including forward multiple scattering) and a strong diffuse multiple scattering background. The D-matrix decomposes the imaging problem into a set of IMs whose degree of complexity increases with their rank (i.e., smaller spatial extent in the focal plane and higher phase distortion in the pupil plane). This last experiment demonstrates the ability of our matrix approach to discriminate between forward multiple scattering paths, which can be taken advantage of for imaging, and the diffuse background, which shall be removed from the final image.

## RESULTS

### Time-gated reflection matrix

The D-matrix concept first relies on the measurement of the time-gated reflection matrix R from the scattering sample. Until now, optical transmission/reflection matrices have always been investigated either in the **k**-space (plane-wave basis) (*20*, *30*) or in the real space (focused basis) (*27*). Here, the R-matrix will be defined between those dual bases. This choice is dictated by our will to go beyond the study of restricted isoplanatic fields of view and tackle space-variant aberrations. Waves produced by nearby points inside a complex medium can generate highly correlated, but tilted, random speckle patterns in the far field (*34*). In a focused basis, this corresponds to a spatially invariant point spread function (PSF) over an area called the IP. As we will see, only a dual-basis matrix can highlight these angular correlations that persist over a restricted spatial domain in the focal plane.

The experimental setup has already been described in a previous work (*27*) and is displayed in fig. S1. The experimental procedure is detailed in Materials and Methods. In a few words, the sample is illuminated through a MO by a set of focused waves (input focusing basis; see Fig. 1A). For each illumination, the amplitude and phase of the reflected wave field is recorded by phase-shifting interferometry on a CCD camera placed in the pupil plane (output pupil basis). A coherent time gating is also applied to select ballistic and snake photons while eliminating a (large) part of the diffuse photons. A set of time-gated reflection coefficients, *R*(**u**_{out}, **r**_{in}), is lastly measured between each virtual source in the focal plane identified by the vector **r**_{in} at the input and each point of the pupil plane identified by the vector **u**_{out} at the output. These coefficients form the reflection matrix R (see Fig. 1D).

The first imaging problem that we consider in this paper deals with an experiment through biological tissues (see Fig. 1A). A positive U.S. Air Force (USAF) 1951 resolution target placed behind an 800-μm-thick rat intestinal tissue is imaged through an immersion objective [40×, NA (numerical aperture), 0.8; Nikon]. The rat intestinal tissue displays a refractive index *n* ∼ 1.4, a scattering mean free path 𝓁* _{s}* of the order of 150 μm and an anisotropy factor

*g*≃ 0.9 (

*35*). The reflection matrix R is measured over a FOI Ω × Ω = 41 × 41 μm

^{2}with

*N*

_{in}= 729 input wavefronts, a spatial sampling δ

*r*

_{in}= 1.6 μm, and an input pupil aperture

. This reduced pupil diameter corresponds to the size of the illumination beam (see fig. S2). At the output, the wave field is recorded over a pupil size of

${\mathcal{D}}_{\text{out}}\times {\mathcal{D}}_{\text{out}}=4.5\times 4.5{\text{mm}}^{2}$ with *N*_{out} = 6084 pixels and a spatial sampling δ*u*_{out} = 68 μm. The corresponding field of view is 60 × 60 μm^{2}. This experimental configuration corresponds to an imaging condition for which time gating guarantees that the reflection matrix contains a fraction of ballistic or snake photons reflected by the resolution target (see fig. S3). However, aberrations are so intense that the full-field image of the resolution target is dominated by the diffuse multiple scattering background (see Fig. 1A4).

Figure 1B1 displays examples of reflected wave fields for several virtual sources **r**_{in}. Each wave field is stored along a column vector and forms the reflection matrix R = [*R*(**u**_{out}, **r**_{in})]. R exhibits a four-dimensional (4D) structure but is concatenated both at the input and output to be displayed in 2D (see fig. S4). The phase of R is displayed in Fig. 1C1. Its spatial and angular correlations in the focal and pupil planes are displayed in Fig. 1 (C2 and C3, respectively). As it could be conjectured from the column vectors displayed in Fig. 1B1, the matrix R only displays short-range correlations. This is quite unexpected as the object to be imaged is deterministic and contained in a single IP. To understand this seemingly randomness of R and reveal its hidden correlations, we now investigate its theoretical expression. The reflection matrix can be expressed as follows (see fig. S5)

(1)or, in terms of matrix coefficients

$$R({\mathbf{u}}_{\text{out}},{\mathbf{r}}_{\text{in}})=\int T({\mathbf{u}}_{\text{out}},\mathbf{r})\mathrm{\gamma}(\mathbf{r}){H}_{\text{in}}(\mathbf{r},{\mathbf{r}}_{\text{in}})d\mathbf{r}$$(2)

H_{in} = [*H*_{in}(**r**, **r**_{in})] is the input focusing matrix. Its columns are none other than the input focal spots centered around each focusing point **r**_{in} (see fig. S5). Under a single scattering assumption, Γ is a diagonal matrix whose elements γ(**r**) map the reflectivity of the object in the focal planes. This object is here assumed to cover the whole FOI. T is the transmission matrix between the focal and pupil planes (see fig. S5). Its elements *T*(**u**_{out}, **r**) describe the propagation of the wave field from a point **r** in the MO focal plane to a detector **u**_{out} in the output pupil plane. Theoretically, the correlation length *r _{P}* of the reflected wave field in the pupil plane scales as λ

*f*/Ω (see section S1), while its correlation length

*r*in the focal plane is dictated by the coherence length of the input focal spots, that is to say, the input diffraction limit,

_{F}, in a strong aberration regime (see section S2). This accounts for the spatial incoherence exhibited by R both at its input (Fig. 1C2) and output (Fig. 1C3), respectively. In the next section, we show how to reveal the hidden correlations in R to, subsequently, extract the transmission matrix T.

### Principle of the distortion matrix

The holy grail for imaging is indeed to have access to this transmission matrix T. Its inversion or pseudoinversion would actually allow to reconstruct a reliable 3D image of the scattering medium, thereby overcoming aberration and multiple scattering effects generated by the medium itself. However, in most imaging configurations, the true transmission matrix T is not accessible as it would require an invasive measurement. The common assumption in wave imaging is thus to consider a homogeneous medium model. The free space transmission matrix T_{0} should then be considered. Its elements *T*_{0}(**u**_{out}, **r**) are simply given by

(3)where *f* is the MO’s focal length and λ the central wavelength.

In this work, we will use T_{0} as a reference matrix. The columns of T_{0} are actually the reflected wave fields that would be obtained in an ideal case, i.e., without aberrations. In the Fourier space, the phase of the complex wave field, or wavefront, is particularly adequate to study the effect of aberrations. Figure 1B compares few examples of reflected wavefronts (columns of R; see Fig. 1B1) with the corresponding ideal wavefronts (columns of T_{0}; see Fig. 1B2). While the latter ones display plane-wave fringes whose orientation and spatial frequency are related to the position **r**_{in} of the input focusing point, the recorded wavefronts consist in a stack of this geometrical component with a distorted phase component induced by the biological tissues. The key idea of this paper is to isolate the latter contribution by subtracting the recorded wavefront by its ideal counterpart. Mathematically, this operation can be expressed as a Hadamard product between R and

(where * stands for phase conjugate)

$$\mathrm{D}=\mathrm{R}\circ {\mathrm{T}}_{0}^{*}$$(4)which, in terms of matrix coefficients, can be written as

$$D({\mathbf{u}}_{\text{out}},{\mathbf{r}}_{\text{in}})=R({\mathbf{u}}_{\text{out}},{\mathbf{r}}_{\text{in}})\times {T}_{0}^{*}({\mathbf{u}}_{\text{out}},{\mathbf{r}}_{\text{in}})$$.(5)

The matrix D = [*D*(**u**_{out}, **r**_{in})] is the so-called distortion matrix. Removing the geometrical component of the reflected wave field in the pupil plane as done in Eq. 4 amounts to a change of reference frame. While the original reflection matrix is recorded in the object’s frame (static object scanned by the input focusing beam; see Fig. 2A), the D-matrix is a reflection matrix in the frame of the input focusing beam (moving object illuminated by a static beam; see Fig. 2B). Physically, this corresponds to a descan of the reflected light as performed in confocal microscopy.

The D-matrix deduced from R is displayed in Fig. 1D1. Compared to R (Fig. 1C1), it exhibits long-range correlations both in the pupil (Fig. 1D3) and focal (Fig. 1D2) planes, respectively. On the one hand, by virtue of the van Cittert Zernike theorem (*36*), the coherence length *d _{P}* of the distorted wave field in the pupil plane scales as λ

*f*/δ

_{in}, with δ

_{in}being the spatial extension of the incoherent input focal spot ∣

*H*

_{in}∣

^{2}(see section S2). On the other hand, its correlation length

*d*in the focal plane corresponds to the size 𝓁

_{F}*of the IP (see section S2). This is illustrated by examples of distorted wave fields displayed in Fig. 1B3. While the original reflected wavefronts did not exhibit any similarity, the distorted component displays similar Fresnel rings whatever the focusing point*

_{c}**r**

_{in}. The D-matrix thus reveals input/output correlations of the wave field that were originally completely hidden in the original R-matrix (Fig. 1C).

### SVD of the distortion matrix

The next step is to extract and take advantage of those field-field correlations for imaging. To that aim, an SVD of the distortion matrix is performed. It consists in writing D as follows

$$\mathrm{D}={\mathrm{U}\mathrm{\Sigma}\mathrm{V}}^{\u2020}$$(6)or, in terms of matrix coefficients

$$D({\mathbf{u}}_{\text{out}},{\mathbf{r}}_{\text{in}})=\sum _{p=1}^{{N}_{\text{in}}}{\mathrm{\sigma}}_{p}{U}_{p}({\mathbf{u}}_{\text{out}}){V}_{p}^{*}({\mathbf{r}}_{\text{in}})$$(7)Σ is a diagonal matrix containing the real positive singular values σ* _{i}* in a decreasing order σ

_{1}> σ

_{2}> . . > σ

_{Nin}. U and V are unitary matrices whose columns,

**U**

_{p}= [

*U*(

_{p}**u**

_{out})] and

**V**

_{p}= [

*V*(

_{p}**r**

_{in})], correspond to the output and input singular vectors, respectively. The symbol † stands for transpose conjugate. Mathematically, the SVD extracts a signal subspace associated with the largest singular values and characterized by an important correlation between its lines and/or columns. In the D-matrix, these correlations are induced by the isoplanicity of the input PSF

*H*

_{in}. The single scattering and forward multiple scattering contributions are expected to lie along the signal subspace since they exhibit a spatial invariance over each IP (

*37*). On the contrary, the diffuse photons induced by the scattering layer ahead of the focal plane give rise to a fully incoherent wave field that will be equally distributed over all the eigenstates of D (

*38*). Hence, the pollution of the signal subspace by the multiple scattering noise scales as the inverse of the number of independent input focusing points mapping each IP,

. A large IP enables the SVD to drastically decrease the multiple-to-single scattering ratio.

To know which of the input or output correlations will dictate the SVD of D, relevant parameters are the numbers of independent speckle grains, *M _{D}* and

*N*, exhibited by D at its input and output, respectively. The correlation degree of the distorted wave field in each plane is actually inversely proportional to the corresponding number of independent speckle grains. In the focal plane,

_{D}*M*is given by the squared ratio between the FOI Ω and the coherence length

_{D}*d*of the distorted wave field in the focal plane

_{F}(8)*d _{F}* is the minimum between the isoplanatic length 𝓁

*and the characteristic fluctuation length 𝓁*

_{c}_{γ}of the object’s reflectivity (see section S2). In the pupil plane, the number

*N*of independent speckle grains scales as (see section S1)

_{D}(9)where δ^{0out} is the diffraction-limited resolution at the output (Eq. 13). The domination of input correlations implies the following condition

(10)

If 𝓁_{γ} > 𝓁* _{c}*, then the last equation can be translated as follows: The number

*M*= (Ω/𝓁

_{D}*)*

_{c}^{2}of IPs supported by the FOI should be smaller than the number

*N*of resolution cells that map each input focusing beam (Eq. 9). As we will see, this strong aberration condition is fulfilled in the experiments presented in this work.

_{D}When input correlations dominate, the effective rank of the signal subspace then corresponds to the number of independent spatial modes required to map the distorted wave field in the focal plane, i.e., the number *M _{D}* of IPs. As shown in section S3, the SVD decomposes the FOI onto a set of orthonormal IMs defined by the input singular vectors

**V**

_{p}. The corresponding output singular vectors

**U**

_{p}yield the associated aberration phase laws in the pupil plane. Their coherent combination can then lead to the retrieval of the transmission matrix T. In the next sections, we will check all these promising properties of D experimentally and see how we can take advantage of it for deep imaging.

### Imaging over a single IP

The reflection and distortion matrices corresponding to the imaging experiment through a thick layer of rat intestine are shown in Fig. 1 (C1 and D1, respectively). The long-range spatial correlations exhibited by D (Fig. 1D2) seem to indicate that the isoplanatic hypothesis is close to being fulfilled in this experiment. The SVD of D confirms this intuition by exhibiting a predominant eigenstate. The corresponding singular vectors **V**_{1} and **U**_{1} are displayed in Fig. 3. The modulus of **V**_{1} displays a contrasted image of the resolution target (Fig. 3B), but its resolution is limited by the low spatial sampling of the illumination scheme. The output singular vector **U**_{1} corresponds to the wavefront induced by a virtual coherent reflector of scattering distribution ∣*H*_{in}(**r** − **r**_{in})∣^{2}, hence located on the optical axis in the focal plane (see Fig. 2C). This virtual scatterer results from a coherent summation of the descanned input focal spots through the SVD process (see section S3). Its phase is made of Fresnel rings mainly induced by the irregular surface of the sample and its index mismatch with the surrounding fluid (Fig. 3D). Its finite support is explained by the finite size δ_{in} of the coherent reflector (Fig. 3C). To make this virtual scatterer point like and retrieve a diffraction-limited image (Fig. 2D), a normalized vector

should be considered, such that

${\tilde{U}}_{1}({\mathbf{u}}_{\text{out}})={U}_{1}({\mathbf{u}}_{\text{out}})/\mid {U}_{1}({\mathbf{u}}_{\text{out}})\mid $.

${\tilde{\mathbf{U}}}_{1}$can be used to build an estimator

$\widehat{\mathrm{T}}$of the transmission matrix between the pupil and focal planes, such that its coefficients read

$${\widehat{T}}_{p}({\mathbf{u}}_{\text{out}},{\mathbf{r}}_{\text{in}})={\tilde{U}}_{p}({\mathbf{u}}_{\text{out}}){T}_{0}({\mathbf{u}}_{\text{out}},{\mathbf{r}}_{\text{in}})$$(11)with *p* = 1 in the present case. This estimator can be used to project the recorded matrix R in the focal basis both at input and output, such that

(12)

The coefficients *R*_{1}(**r**_{out}, **r**_{in}) are the impulse responses between each input focusing point **r**_{in} and each output imaging point **r**_{out}. In other words, once reshaped in 2D, each column of R_{1} yields the PSF of the imaging system at the input focusing point **r**_{in}. The PSF for an input focusing point on the optical axis (**r**_{in} = 0) is displayed in Fig. 3F. For the sake of comparison, the corresponding initial focal spot is displayed in Fig. 3E. The latter one is extracted from the focused matrix R_{0} deduced from R using

. While the initial PSF exhibits a random speckle pattern (Fig. 3E), the PSF after correction shows a nearly diffraction-limited focal spot with almost all the energy concentrated in the vicinity of the incident focusing point (Fig. 3F). The apparent width of this PSF yields an estimation of the local output resolution δ_{out} at **r**_{in}. Here, δ_{out} goes from 20 μm on the raw data (Fig. 3E) to 1.2 μm after the matrix correction (Fig. 3F). This value should be compared to the diffraction-limited resolution

(13)with

${\text{NA}}_{\text{out}}={\mathcal{D}}_{\text{out}}/(2f)=0.45$being the output NA. The numerical application of this formula yields

${\mathrm{\delta}}_{\text{out}}^{0}\simeq 0.9\mathrm{\mu}\mathrm{m}$ in our experimental configuration. The mismatch between δ_{out} and

comes from the noisy aspect of **U**_{1} at large spatial frequencies (see Fig. 3D), which prevents from an efficient aberration compensation over the whole NA.

If the spatial sampling was equivalent at input and output, then a confocal image could be deduced from the diagonal elements (**r**_{in} = **r**_{out}) of R_{0} and R_{1} (*27*). Here, as a sparse illumination scheme has been used

, a full-field image is considered and obtained by summing R over its input elements

$${\mathcal{F}}_{p}({\mathrm{r}}_{\text{out}})={\displaystyle \sum _{{\mathbf{r}}_{\text{in}}}}\mid {R}_{p}({\mathbf{r}}_{\text{out}},{\mathbf{r}}_{\text{in}})\mid $$(14)with *p* = 0 or 1 here. The corresponding images ℱ_{0} and ℱ_{1} are displayed in Fig. 3 (G and H, respectively). While the patterns of the resolution target are hardly visible on the raw image, the D-matrix approach provides a highly contrasted image of the target. To quantify this gain in image quality, the Strehl ratio is a relevant parameter (*39*). It is defined as the ratio of the PSF peak intensity with and without aberration. Equivalently, it can also be defined in the pupil plane as the squared magnitude of the mean aberration phase factor. Its initial value

can thus be directly derived from the D-matrix coefficients

$${\mathcal{S}}_{0}={\mid \mathrm{\u3008}\text{exp}(j\text{arg}\{D({\mathbf{u}}_{\text{out}},{\mathbf{r}}_{\text{in}}){V}_{1}({\mathbf{r}}_{\text{in}})\})\mathrm{\u3009}\mid}^{2}$$(15)where the symbol 〈⋯〉 denotes an average over **u**_{out} and **r**_{in}. In the present case, the original Strehl ratio is S_{0} = 8 × 10^{−5}. This experiment corresponds to imaging conditions far from being in the range of operation of conventional AO and explains why the patterns of the resolution target are so hardly visible on the raw image (Fig. 3G). The Strehl ratio

after the **U**_{1} correction can be directly extracted from the SVD of D (Eq. 7)

(16)

The **D**-matrix correction leads to a Strehl ratio S_{1} = 3 × 10^{−3}. However, Eq. 16 gives the same weight to bright and dark areas of the resolution target in the focal and pupil planes. One possibility is to consider a weighted average instead of Eq. 16 by the object reflectivity ∣*V*_{1}(**r**_{in})∣^{2}. This weighted Strehl ratio

then reaches the value of 1.1 × 10^{−2}. Such a Strehl ratio value is relatively low, but it should be kept in mind that the distortion matrix is associated with a PSF in reflection that convolves the transmit and receives PSFs. Our measurement of the Strehl ratio is thus degraded by (i) the subsistence of input aberrations and (ii) the presence of a diffuse multiple scattering background that acts here as an additive noise. Note, however, that the gain in terms of Strehl ratio is absolute; this is the relevant quantity to assess the benefit of our matrix approach. This gain here is spectacular

and accounts for the satisfying image of the resolution target obtained after the D-matrix correction (see Fig. 3H). Figure S3 shows how this drastic improvement of the Strehl ratio allows us to push back the imaging depth limit from 450 μm to almost 1 mm.

This first experiment demonstrates the benefit of the D-matrix in the simple case of an FOI containing a single IP. In the next section, the case of multiple IPs is tackled.

### Imaging over multiple IPs

The first element of the group 6 in the resolution target is now imaged through an aberrating layer consisting in a rough plastic sheet placed *d* = 1 mm above the resolution target (USAF 1951) (see Fig. 4A). The reflection matrix R is measured over a FOI of 57 × 57 μm^{2} with *N*_{in} = 441 input wavefronts, a spatial sampling δ*r*_{in} = 2.85 μm, and an input pupil aperture

. At the output, the wave field is recorded over a pupil size of

${\mathcal{D}}_{\text{out}}\times {\mathcal{D}}_{\text{out}}=2\times 2{\text{mm}}^{2}$ with *N*_{out} = 12,321 pixels and a spatial sampling δu_{out} = 18 μm.

The full-field image F_{0} (Eq. 14) and an example of PSF (Eq. 12) are displayed in Fig. 4 (A and B, respectively). The PSF is strongly degraded with a characteristic focal spot dimension δ_{out} ∼ 45 μm. This PSF dimension allows an estimation of the coherence length 𝓁* _{c}* of the aberrating layer. Under a thin phase screen model (

*37*), the IP dimension 𝓁

*coincides with the coherence length of the aberration transmittance. It turns out that the PSF width is inversely proportional to 𝓁*

_{c}*in this experiment: δ*

_{c}_{out}∼ λ

*d*/𝓁

*. The IP size and the number of IPs supported by the FOI can be deduced from the PSF width δ: 𝓁*

_{c}*∼ 18 μm and*

_{c}*M*∼ (Ω/𝓁

_{D}*)*

_{c}^{2}∼ 10.

A D-matrix is deduced from R (Eq. 4). Its analysis leads to the following estimation of the initial Strehl ratio:

${\mathcal{S}}_{0}^{\prime}=1.6\times {10}^{-6}$(Eq. 15). This particularly strong aberration level accounts for the highly blurred aspect of the full-field image in Fig. 4A. This experimental situation is thus particularly extreme, even almost hopeless, for a successful imaging of the resolution target. Yet, the SVD of D will provide the solution.

Figure 4D displays the histogram of the normalized singular values

${\tilde{\mathrm{\sigma}}}_{p}={\mathrm{\sigma}}_{p}/{\sum}_{i=1}^{{N}_{\text{in}}}{\mathrm{\sigma}}_{i}$. If recorded data were not corrupted by experimental noise, then the matrix would be of effective rank *M _{D}*. We could use all the eigenstates of D associated with nonzero singular values to retrieve an image of the object. In Fig. 4D, only few singular values seem to emerge from the noise background. Hence, it is difficult to determine the number of eigenstates that we need to consider to properly reconstruct an image of the object. This issue can be circumvented by computing the Shannon entropy H of the singular values (

*40*,

*41*), such that

(17)Shannon entropy delivers the maximally noncommittal dataset at a given signal-to-noise ratio, that is to say, the most information with the least artifact. The Shannon entropy can be used as an indicator of how many eigenstates are needed to build an adequate image of the object without being affected by experimental noise.

The singular values of Fig. 4D have an entropy *H*≃ 8.4. Hence, only the eight first singular states shall be considered. Figure 4G displays the phase of the four first singular vectors **U**_{p} in the pupil plane. They display phase distortions whose typical coherence length *u _{c}* scales as

*f*𝓁

*/*

_{c}*d*∼ 100 μm. The phase conjugation of these singular vectors should compensate for the detrimental effect of the aberrating layer in different parts of the FOI. A set of focused reflection matrices R

_{p}can be deduced (Eq. 12). Figure 4F displays an example of corrected PSF extracted from a column of R

_{1}. Its comparison with the original PSF in Fig. 4C shows how the phase conjugation of

**U**

_{1}allows one to compensate for the aberrations at this incident focusing point. On the one hand, the PSF width is narrowed by a factor of 20, going from δ

_{out}∼ 45 μm to 2.25 μm. The latter value should be compared with the diffraction-limited resolution

(Eq. 13) in our experimental conditions. The Strehl ratio is increased by a factor of 2.2 ×10^{3} to reach the final value

(Eq. 16). Again, this value of the Strehl ratio is probably underestimated because of input aberrations and multiple scattering.

Confocal images can be deduced from the focused reflection matrices R_{p}

(18)where *l _{p}* is the aperture of the numerical confocal pinhole (

*27*). This finite aperture enables an average of the output image over neighbor incident focusing points to smooth out the sparse illuminations. Figure 4H displays the confocal images I

*for*

_{p}*l*= 2 μm. For a specular object such as a resolution target, the SVD has indeed the property of decomposing into a set of orthogonal IMs of spatial period 𝓁

_{p}*(see section S3). Their shape depends on the autocorrelation function of the aberrating phase screen. A general trend is that the spatial frequency content of the eigenvectors increases with their rank. If this function presents an exponential or sinc model, then the FOI will be spatially decomposed into sinusoidal wave functions (*

_{c}*42*) analogous to optical fiber modes or to prolate spheroidal wave functions (

*43*), respectively. Here, the autocorrelation function of the aberrating phase displays a Gaussian-like shape. The FOI is thus spatially mapped onto Hermite-Gaussian wave functions analogous to laser cavity modes (

*44*).

The normalized pupil singular vectors

${\tilde{\mathbf{U}}}_{p}$yield a set of orthogonal phase transmittances that map aberrations onto each IM. A coherent combination of these singular vectors should lead, in principle, to a satisfying estimator of the transmission matrix (see section S3)

$$\widehat{\mathbf{T}}={\displaystyle \sum _{p=1}^{\mathcal{H}({\tilde{\mathrm{\sigma}}}_{i})}}{\tilde{\mathbf{U}}}_{\mathrm{p}}\circ {\mathbf{T}}_{0}$$(19)In practice, a final image I of the resolution target can be obtained by summing the previous IMs I_{p}

(20)The final result is displayed in Fig. 4E. The comparison with the initial full-field image (Fig. 4B) illustrates the benefit of the D-matrix approach. Spatially varying aberrations are overcome, and a contrasted image of the resolution target is recovered over the whole FOI. This experiment demonstrates how the D-matrix enables a decomposition of the FOI into several IMs and a mapping of each of them onto orthonormal distorted phase laws. However, this demonstration has been restricted to the case of a 2D aberrating phase layer. In the next section, we consider the case of a cornea with deteriorated transparency as a 3D aberrating and scattering structure.

### Imaging through a hazy cornea

The experimental configuration is displayed in Fig. 5A. The number “3” of the group 5 in the resolution target is imaged through a 700-μm-thick edematous nonhuman primate cornea. The reflection matrix R is measured over an FOI of 52 × 52 μm^{2} by means of *N*_{in} = 625 input wavefronts, a spatial sampling δ*r*_{in} = 2.1 μm, and an input pupil aperture

. At the output, the wave field is recorded over an output pupil size

${\mathcal{D}}_{\text{out}}\times {\mathcal{D}}_{\text{out}}=2\times 2{\text{mm}}^{2}$ with *N*_{out} = 1296 pixels and a spatial sampling length δ*u*_{out} = 56 μm . Figure 5C displays the confocal image I_{0} deduced from R_{0} with *l _{p}* = 1 μm (Eq. 18). Multiple scattering and aberrations induced by the cornea induce a random speckle-reflected wave field that prevents from imaging the resolution target. On the contrary, as we will see, the D-matrix analysis allows us to nicely recover the pattern 3 of the resolution target (see Fig. 5D).

Figure 5C displays the spectrum of the singular values

${\tilde{\mathrm{\sigma}}}_{p}$. The first singular value emerges from the rest of the spectrum, but it is difficult to know until which rank the eigenstates can be considered as belonging to the signal subspace. As previously, the Shannon entropy of the singular values yields an unambiguous answer:

$\mathcal{H}({\tilde{\mathrm{\sigma}}}_{i})=10.7$. The 11th first singular states should thus be considered. Figure 5E displays the 1st, 6th, and 11th singular vectors **U*** _{p}* in the pupil plane. The complexity of the wavefront distortion, i.e., their spatial frequency content, increases with the rank of the corresponding singular values. The corresponding IMs I

*(Eq. 18) are displayed in Fig. 5F. While the first singular vector*

_{p}**U**

_{1}allows a wide-field correction of low-order aberration, the higher-rank singular vectors are associated with high-order aberrations that are effective over IMs of smaller dimension. In Fig. 5D, the whole spatial frequency spectrum of wavefront distortions is lastly compensated by smartly combining the confocal images I

*associated with each singular state from D’s signal subspace (Eq. 20). The comparison of the initial (Fig. 5C) and final (Fig. 5D) images is spectacular with a Strehl ratio gain*

_{p}. The comparison of I (Fig. 5D) and I_{1} (see the first inset of Fig. 5F) illustrates the benefit of a matrix approach of aberration correction compared to conventional AO, since the latter one would yield I_{1} in theory.

This decomposition of complex aberration phase laws over a set of IMs opens important perspectives for aberrometry. It actually goes well beyond state-of-the-art techniques that basically consist in a simple projection over a set of Zernike polynomials. Moreover, an estimator of the single-to-multiple scattering ratio (SMR) can be built on the relative energy between the signal and noise subspaces of D

$$\text{SMR}=\frac{{\displaystyle {\sum}_{p=1}^{\mathcal{H}({\mathrm{\sigma}}_{i})}}{\mathrm{\sigma}}_{p}^{2}}{{\displaystyle {\sum}_{i=\mathcal{H}({\mathrm{\sigma}}_{i})+1}^{{N}_{\mathit{\text{in}}}}}{\mathrm{\sigma}}_{p}^{2}}$$(21)The SMR can actually be a quantitative biomarker of the corneal opacification or a quantitative measure of corneal transparency (*45*). On the basis of a fit with a recent analytical study of the SMR (*38*), the cornea thickness *L* can be estimated in terms of scattering mean free path 𝓁* _{s}*:

*L*∼ 9𝓁

*(see fig. S3). As the corneal thickness is known (*

_{s}*L*= 700 μm), the scattering mean free path can be deduced: 𝓁

*∼ 80 μm. This value is in excellent agreement with recent ex vivo measurements of 𝓁*

_{s}*in pathological corneas with compromised transparency (*

_{s}*45*). The value of 9𝓁

*highlights the difficult experimental conditions under which the imaging of the resolution target has been successfully achieved.*

_{s}In conclusion, this last experiment shows the potential of a matrix approach for eye aberrometry and turbidimetry, such as for improved quality control of donor tissue assessment before corneal transplantation (*45*). Of course, this method is by no means limited to ophthalmic applications. It can be applied to the characterization of any kind of biological tissues provided that we have access to the associated reflection matrix.

## DISCUSSION

Here, we present a novel and noninvasive method for aberration compensation and diffraction-limited imaging at large optical depths. This approach relies on a new operator, the so-called distortion matrix, that connects a set of input focusing points with the distorted component of the reflected the wave field in the pupil plane. This operator connecting position and spatial frequency has some analogy with the Wigner distribution function (*46*). However, the Wigner transform applies to a single variable of a function, i.e., to a single vector in a discrete formalism. Here, our position-momentum analysis is performed between the input and output of a reflection matrix.

The concept of distortion matrix is to measure the backscattered waves in a descanned frame while scanning the sample with focused illuminations. This approach has some similarity with a previous AO approach (*10*) in its hardware configuration. The main difference is that, in this study, wavefronts are averaged by the Shack-Hartmann type of analysis and this AO approach thus relies on an isoplanatic condition. Here lies one of the strengths of our approach. While conventional methods estimate the aberrated wavefront for a single location or averaged over the whole FOI, we propose to study the spatial and angular correlations of the distortion operator through an SVD. In this manner, we demonstrate the efficient compensation of both low- and high-order aberrations over multiple IPs. Moreover, our approach relies on the Shannon entropy that provides an objective criterion to determine the number of IPs supported by the FOI. This is in contrast with recent works based on a far-field reflection matrix in which the FOI was arbitrarily divided into subareas where different corrections were applied (*47*, *48*).

Besides aberration correction, our approach leverages the correlations of the output wave field to get rid of the multiple scattering background. The latter contribution is actually spatially incoherent. It thus mainly lies along the noise subspace of the D-matrix. Thanks to these features, we were able to image through almost 10 scattering mean free paths of biological tissues, which is beyond the imaging depth of conventional OCT systems for these specimens (see fig. S3). Compared to the previously developed smart-OCT method that was able to detect few bright scatterers at large penetration depth (12𝓁* _{s}*) (

*27*), the D-matrix approach yields a direct image of the sample reflectivity at a diffraction-limited resolution. In addition, our approach enables to quantitatively estimate the amount of multiply scattered light. Combined with a conventional image, this parameter is of importance for characterization purposes.

The distortion operator thus opens a new route toward real-time deep optical imaging of biological tissues. In that respect, the experimental setup and procedure used in this paper are clearly perfectible. While postprocess operations take less than 1 min on a regular laptop, the main limitation in the current experimental configuration is the acquisition time. In particular, the scanning illumination scheme was not optimized because of the SLM speed. While the use of a galvanometric mirror or a high-speed deformable mirror would drastically reduce the acquisition time at the cost of a more complex setup, we counteracted this issue with a sparse illumination. However, this, in return, limited the available number of angular degrees of freedom at the input, which prevents us from an aberration correction of the incident wave field. By optimizing the experimental apparatus and acquisition scheme, large reflection matrices can be measured in a few seconds. For instance, Yoon *et al*. (*48*) recently demonstrated the acquisition of a 10,000-mode matrix in 15 s with the same degree of control for the incident and reflected waves. In that case, a simultaneous correction of aberrations at the input and output is absolutely possible under the distortion matrix approach by alternatively projecting the incident and reflected wave fields in the focal and pupil planes. In view of 3D imaging, our approach can also be coupled to computational AO techniques (*12*) to tackle depth-dependent aberrations and restore a diffraction-limited resolution in all directions. An alternative is to switch from a scanning to a full-field illumination scheme. A measurement of the coherent reflection matrix R can be performed under a spatially incoherent illumination (*49*, *50*). This full-field configuration would allow to record the reflection matrix over millimetric volumes in a moderate acquisition time.

Last, we used a negative resolution target as the sample to be imaged in this work. The reason is that this highly contrasted object was the ideal specimen to clearly highlight the issue of multiple isoplanatic areas. Beyond the proof-of-concept experiments presented in this article, a direct imaging of biological specimens over large penetration depth will be the next step. The assumption on which our method is based (Eq. 10) can easily be met in biological tissues since a strong aberration regime takes place beyond a few scattering mean free paths. Note also that, even when this condition is not fulfilled and far-field correlations dominate, the distortion matrix approach can still work but the FOI has to be beforehand subdivided into individual IPs (*47*, *51*). The ability of identifying multiple IPs will also be particularly promising to map the specimen-induced aberration and the SMR. Aside from aberrometry and/or turbidimetry, future in vivo implementations of our approach have implications beyond that of ocular media characterization, most notably for imaging through nontransparent ocular media (e.g., retinal imaging through a turbid cornea or through cataract opacities) (*52*).

In summary, we have introduced, in this work, a new operator, the so-called distortion matrix D, which reveals the hidden correlations of the reflected wave field. This matrix results from the mismatch between the phase of the recorded reflection matrix and those of a reference matrix that would be obtained in an ideal configuration. As shown in this paper, D gives access to the noninvasive transmission matrix between each sensor and each voxel of the FOI. Then, by solving the corresponding inverse problem, an image of a scattering sample can be obtained as if the medium ahead was made transparent. The D-matrix concept is very general. It can be extended to any kind of waves and experimental configurations for which a measurement of the amplitude and phase of the reflected wave field is possible under multiple illuminations (*53*–*56*). A recent work actually demonstrates the benefits of this concept for ultrasound imaging in a random scattering regime (*51*). This D-matrix concept thus opens a new route toward a global and noninvasive matrix approach of deep imaging in biological tissues.

**Acknowledgments: **We wish to thank L. Cobus, W. Lambert, P. Balondrade, and S. Meimon for fruitful discussions. **Funding:** The authors are grateful for the funding provided by Labex WIFI (Laboratory of Excellence within the French Program Investments for the Future; ANR-10-LABX-24 and ANR-10-IDEX-0001-02 PSL*). A.B. acknowledges financial support from the French “Direction Générale de l’Armement” (DGA). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement nos. 610110 and 819261, HELMHOLTZ* and REMINISCENCE projects, respectively). K.I. acknowledges financial support from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement no. 709104. **Author contributions:** A.A. initiated and supervised the project. A.B. built the experimental setup and performed the experiments. K.I., A.C.B., and M.F. initiated the ophthalmic application. K.I. provided corneal samples and guidance on the ophthalmic experiment. A.B., V.B., and A.A. analyzed the experiments. V.B. and A.A. performed the theoretical study. A.B. and A.A. prepared the manuscript. A.B., V.B., K.I., A.C.B., M.F., and A.A. discussed the results and contributed to finalizing the manuscript. **Competing interests:** A.A., M.F., A.C.B., A.B., and V.B. are inventors on a patent related to this work held by CNRS (no. WO2020016249, published January 2020). All authors declare that they have no other competing interests. **Data and materials availability:** All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.