Distortion matrix concept for deep optical imaging in scattering media
Amaury Badon, Victor Barolle, Kristina Irsch, A. Claude Boccara, Mathias Fink, Alexandre Aubry
DDistortion matrix concept for deep optical imaging inscattering media
Amaury Badon , Victor Barolle , Kristina Irsch , , Albert C. Boccara ,Mathias Fink , and Alexandre Aubry , ∗ Institut Langevin, ESPCI Paris, PSL University, CNRS, 1 rue Jussieu, 75005 Paris, France Vision Institute / Quinze-Vingts National Eye HospitalSorbonne University, CNRS UMR 7210, INSERM U 068, 17 rue Moreau, 75012 Paris, France The Wilmer Eye Institute, The Johns Hopkins University School of Medicine, Baltimore, MD, USA ∗ To whom correspondence should be addressed; E-mail: [email protected].
Summary Sentence : The distortion matrix overcomes aberrations and multiple scattering,thereby enabling ultra-deep and wide-field optical imaging.
Abstract
In optical imaging, light propagation is affected by the inhomogeneities of the medium.Sample-induced aberrations and multiple scattering can strongly degrade the image reso-lution and contrast. Based on a dynamic correction of the incident and/or reflected wave-fronts, adaptive optics has been employed to compensate for those aberrations. However, itonly applies to spatially-invariant aberrations or to thin aberrating layers. Here, we proposea global and non-invasive approach based on the distortion matrix concept. This matrixbasically connects any focusing point of the image with the distorted part of its wave-frontin reflection. A singular value decomposition of the distortion matrix allows to correctfor high-order aberrations and forward multiple scattering over multiple isoplanatic modes.Proof-of-concept experiments are performed through biological tissues including a turbidcornea. We demonstrate a Strehl ratio enhancement up to 2500 and recover a diffraction-limited resolution until a depth of ten scattering mean free paths. a r X i v : . [ phy s i c s . op ti c s ] J u l NTRODUCTION
For decades, optical microscopy has been a vital tool in biomedical research to observe livespecimens with a sub-micron resolution and with minimal invasiveness. Yet, imaging con-ditions required for such exquisite performances are rarely gathered. For instance, both theresolution and the contrast drop as the imaging depth increases inside a biological tissue. Thisobservation is a consequence of the spatial variations of the specimen’s refractive index that dis-tort the wave-front of both the incoming and outgoing light. When these variations exhibit lowspatial frequencies we use the term aberrations while scattering describes the effect of the higherspatial variations. Both these effects limit the use of conventional microscopy to shallow depthsor to semi-transparent specimens. Imaging deeper requires to simultaneously compensate forthese 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 ). In-deed, astronomers faced the same impediment as fluctuations in the atmosphere severely distortthe wave-front of the light coming from stars and prevent to obtain a diffraction-limited stel-lar image. Astronomers then proposed to measure these distortions using a wave-front sensorand to counterbalance it with a dynamic programmable element such as deformable mirrors.Following this concept and the development of deformable mirrors with increasing number ofelements, AO already demonstrated its benefits in various imaging techniques such as digitalholography (
3, 4 ), confocal microscopy (
5, 6 ), two-photon microscopy ( ), or optical coher-ence tomography (OCT) (
11, 12 ). Unfortunately AO methods usually require a guide star orare based on an image sharpness metric. Additionally, they are limited to a small region calledthe 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 the2ase of multiple IPs. This issue is particularly decisive for deep imaging where IP size becomesextremely tiny: < µ m beyond a depth of 1 mm ( ). Note that multi-conjugate AO can dealwith multiple IPs but this is at the price of a much more complex optical set up ( ).On the other hand, multiple scattering was long thought to be too complex to be compen-sated. For deep imaging, a gating mechanism is generally used to reject the multiply-scatteredphotons and capture only the ballistic light. This gating can be spatial ( ) as in confocalmicroscopy or temporal ( ) as in optical coherence tomography, but they are still depth lim-ited as they rely on the exponentially attenuated ballistic light. In a pioneering experiment,Vellekoop and Mosk demonstrated in 2007 the possibility to restore a diffraction-limited spotthrough a scattering medium by properly shaping the incoming light ( ). Subsequently, a ma-trix approach of light propagation through complex media was developed ( ). Relying on themeasurement 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 ac-cess to the transmission matrix allows taking advantage of multiple scattering for optimal lightfocusing ( ) and communication across a diffusive layer (
21, 22 ) or a multimode fiber (
23, 24 ).However, a transmission configuration is not adapted to non-invasive and/or in-vivo imaging ofbiological media. An epi-detection geometry should thus be considered ( ). During the lastfew years, the reflection matrix R had been investigated to perform selective focusing/ detec-tion (
26, 27 ) or energy delivery (
28, 29 ) through strongly scattering media. With regards to thespecific purpose of imaging, the matrix approach has been recently used to implement AO toolsin post-processing. The single scattering component of the reflected wave-field through biolog-ical tissues has been enhanced in depth by compensating for high-order aberrations (
30, 31 ).In this paper, we propose to go beyond a matrix approach of AO by introducing a novel oper-ator: the so-called distortion matrix D . Unlike previous works that investigated R either in thefocal plane ( ) or the pupil plane (
26, 30, 31 ), we here consider the medium response between3hose dual bases (
32, 33 ). Unlike R , the D -matrix does not consider the reflected wave-field asa building block but its deviation from an ideal wave-front that would be obtained in absence ofaberrations and without multiple scattering. This operation may seem trivial but it dramaticallyhighlights the input/output correlations of the wave-field. While the canonical reflection matrixexhibits a random feature in a turbid medium, the distortion matrix displays strong field-fieldcorrelations over each IP. Thanks to this new operator, some relevant results of information the-ory 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 extract the associated wave-front distortion in the pupil plane. The Shannon entropy H ofthe singular values allows one to define the effective rank of the imaging problem. A combina-tion of the H first eigenstates yields an image of the focal plane with an excellent contrast anda diffraction-limited resolution as if the medium ahead was made perfectly transparent.Several experiments with an increasing order of complexity are presented to demonstratethe benefits of the D -matrix for optical imaging in turbid media. For sake of simplicity, thefirst 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 be-hind it. Then, a second proof-of-concept experiment considers a thin but strong aberratinglayer introduced between the microscope objective and a resolution target. This imaging con-figuration involves a spatially-varying aberration across the FOI ( i.e several IPs). At last, wedescribe an imaging experiment through a turbid non-human primate cornea that induces high-order aberrations (including forward multiple scattering) and a strong diffuse multiple scatteringbackground. The D -matrix decomposes the imaging problem into a set of IMs whose degreeof complexity increases with their rank ( i.e. smaller spatial extent in the focal plane and higherphase distortion in the pupil plane). This last experiment demonstrates the ability of our matrixapproach to discriminate between forward multiple scattering paths, that can be taken advantage4f for imaging, and the diffuse background, that 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 alwaysbeen investigated either in the k -space (plane-wave basis) (
20, 30 ) or in the real space (focusedbasis) ( ). Here the R -matrix will be defined between those dual bases. This choice is dictatedby our will to go beyond the study of restricted isoplanatic fields of view and tackle space-variant aberrations. Indeed, waves produced by nearby points inside a complex medium cangenerate highly correlated, but tilted, random speckle patterns in the far-field ( ). In a focusedbasis, this corresponds to a spatially invariant point spread function (PSF) over an area calledthe isoplanatic patch. As we will see, only a dual-basis matrix can highlight these angularcorrelations that persist over a restricted spatial domain in the focal plane.The experimental set-up has already been described in a previous work ( ) and is displayedin Fig. S1. The experimental procedure is detailed in the Methods section. In a few words, thesample is illuminated through a microscope objective (MO) by a set of focused waves (inputfocusing basis) (see Fig. 1A). For each illumination, the amplitude and phase of the reflectedwave-field is recorded by phase-shifting interferometry on a CCD camera placed in the pupilplane (output pupil basis). A coherent time gating is also applied in order to select ballisticand snake photons while eliminating a (large) part of the diffuse photons. A set of time gatedreflection coefficients, R ( u out , r in ) , is finally measured between each virtual source in the focalplane identified by the vectors r in at the input and each point of the pupil plane u out at the output.These coefficients form the reflection matrix R (see Fig. 1D).The first imaging problem we consider in this paper deals with an experiment through bio-5 eflection Geometric Distortion = = ++= + Input pupil
Illumination 1 Illumination 2 A. MO Focal plane r r Illumination 3Illumination 3Illumination 2Illumination 1
Image planeOutput pupil
C. B. D. x out y o u t u i n = + u out D ( r in , u out ) r
100 dB-30dB π-π ⊗ r in u out ⊗ ⊗ u out r in ⊗
15 µmNA=0.5
A2A1 A3 A4B1 B2 B3C1C2C3 D1 D2D3
15 µm
15 µm
15 µmNA=0.5 d P d F r P r F Figure 1:
Principle of the distortion matrix approach. ( A ) A resolution target (USAF 1951)is positioned underneath a 800- µ m-thick sample of rat intestine ( A1 ). In scanning microscopy,raster scanning in the focal plane is obtained using a set of plane wave illuminations in theinput pupil ( A2 ). In presence of sample-induced aberrations, the detected intensity will exhibita much larger extent compared to the ideal PSF ( A3 ). The resulting full-field image displaysa low contrast and a reduced resolution ( A4 ). ( B ), In the output pupil plane, the phase of thereflected wave-field ( B1 ) can be split into a diffraction ( B2 ) and a distortion ( B3 ) term. ( C , D )The reflected distorted wave-fields can be stored along column vectors to form the reflectionand distortion matrices, R and D , respectively. The phase of R and D is displayed in ( C1 )and ( D1 ), respectively. The auto-correlations of the complex reflected/distorted wave-fields arecomputed in the focal ( C2 / D2 , see Section S2) and in the pupil ( C3 / D3 , see Section S1) planes,both in dB. All the data shown here are extracted from the rat intestine imaging experiment.Photo Credit: Amaury Badon, CNRS. 6ogical tissues (see Fig. 1A). A positive U.S. Air Force (USAF) 1951 resolution target placedbehind an 800- µ m-thick of 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 ∼ . , a scattering mean free path (cid:96) s of the order of 150 µ m and an anisotropy factor g (cid:39) . ( ). The reflection matrix R is measured over a FOI Ω × Ω = 41 × µ m with N in =
729 input wave-fronts, a spatial sampling δr in = 1 . µ m and an input pupil aperture D in × D in = 1 . × . mm . This reduced pupil diameter corresponds to the size of the il-lumination beam (see Fig. S2). At the output, the wave-field is recorded over a pupil size of D out × D out = 4 . × . mm with N out = 6084 pixels and a spatial sampling δu out = 68 µ m.The corresponding field-of-view is × µ m . This experimental configuration correspondsto an imaging condition for which time gating guarantees that the reflection matrix contains afraction 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 thediffuse multiple scattering background (see Fig. 1A4).Figure 1B1 displays examples of reflected wave-fields for several virtual sources r in . Eachwave-field is stored along a column vector and forms the reflection matrix R = [ R ( u out , r in )] . R exhibits a 4D-structure but is concatenated both at the input and output to be displayed in2D (see Fig. S4). The phase of R is displayed in Fig. 1C1. Its spatial and angular correlationsin the focal and pupil planes are displayed in Figs. 1C2 and C3, respectively. As it couldbe conjectured from the column vectors displayed in Fig. 1B1, the matrix R only displaysshort-range correlations. This is quite surprising as the object to be imaged is deterministicand contained in a single IP. To understand this seemingly randomness of R and reveal itshidden correlations, we now investigate its theoretical expression. The reflection matrix can beexpressed as follows (see Fig. S5): R = T × Γ × H in (1)7r, in terms of matrix coefficients, R ( u out , r in ) = (cid:90) T ( u out , r ) γ ( r ) H in ( r , r in ) d r (2) H in = [ H in ( r , r in )] is the input focusing matrix. Its columns are none other than the inputfocal spots centered around each focusing point r in (see Fig. S5). Under a single scatteringassumption, Γ is a diagonal matrix whose elements γ ( r ) map the reflectivity of the object inthe focal planes. This object is here assumed to cover the whole FOI. T is the transmissionmatrix between the focal and pupil planes (see Fig. S5). Its elements T ( u out , r ) describe thepropagation of the wave-field from a point r in the MO focal plane to a detector u out in theoutput pupil plane. Theoretically, the correlation length r P of the reflected wave-field in thepupil plane scales as λf / Ω (see Section S1) while its correlation length r F in the focal planeis dictated by the coherence length of the input focal spots, that is to say the input diffractionlimit, δ ∼ λf / D in , in a strong aberration regime (see Section S2). This accounts for the spatialincoherence exhibited by R both at its input (Fig. 1C2) and output (Fig. 1C3), respectively. Inthe next section, we show how to reveal the hidden correlations in R in order 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 inversionor pseudo-inversion would actually allow to reconstruct a reliable 3D image of the scatter-ing medium, thereby overcoming aberration and multiple scattering effects generated by themedium itself. However, in most imaging configurations, the true transmission matrix T isnot accessible as it would require an invasive measurement. The common assumption in waveimaging, is thus to consider an homogeneous medium model. The free space transmission8atrix T should then be considered. Its elements T ( u out , r ) are simply given by T ( u out , r in ) = 1 jλf exp (cid:20) j πλf u out . r (cid:21) (3)where f is the MO’s focal length and λ the central wavelength.In this work, we will use T as a reference matrix. The columns of T are actually thereflected wave-fields that would be obtained in an ideal case, i.e without aberrations. In theFourier space, the phase of the complex wave-field, or wave-front, is particularly adequate tostudy the effect of aberrations. Figure 1B compares few examples of reflected wave-fronts(columns of R , see Fig. 1B1) with the corresponding ideal wave-fronts (columns of T , seeFig. 1B2). While the latter ones display plane wave fringes whose orientation and spatial fre-quency is related to the position r in of the input focusing point, the recorded wave-fronts consistin a stack of this geometrical component with a distorted phase component induced by the bio-logical tissues. The key idea of this paper is to isolate the latter contribution by subtracting therecorded wave-front by its ideal counterpart. Mathematically, this operation can be expressedas a Hadamard product between R and T ∗ (where ∗ stand for phase conjugate), D = R ◦ T ∗ (4)which, in term of matrix coefficients, can be written as D ( u out , r in ) = R ( u out , r in ) × T ∗ ( u out , r in ) (5)The matrix D = [ D ( u out , r in )] is the so-called distortion matrix. Removing the geometricalcomponent of the reflected wave-field in the pupil plane as done in Eq.4 amounts to a changeof reference frame. While the original reflection matrix is recorded in the object’s frame (staticobject scanned by the input focusing beam, see Fig. 2A), the D -matrix is a reflection ma-trix in the frame of the input focusing beam (moving object illuminated by a static beam, see9 irtual sample scanningbeam scanning CCD planeMOfocal plane de-scan coherent sum
A B C
R-matrix SVD of DD-matrix normalization
Normalization D virtual coherent reflector virtual point-like reflector ^ V p U p U p V p ^ Figure 2:
Extracting the aberration transmittance from the distortion matrix D . ( A ) Therecording of the R -matrix consists in scanning the objects with a moving input focusing beam.( B ) The removal of the geometric component in each reflected wave-front [Eq.4] amounts torecenter each incident focal spot at the origin. The D -matrix is equivalent to the reflectionmatrix for a moving object. ( C ) The SVD of D leads to a coherent sum of the distorted wave-fronts in the pupil plane. A coherent reflector is virtually synthesized in the focal plane andthe corresponding wave-front emerges along the output singular vector U . The correspondingimage of the object is provided by the first input singular vector V but its resolution is dictatedby the width δ in of the input focusing beam. ( C ) A normalization of U in the pupil planemakes the virtual scatterer point-like. The corresponding input singular vector ˆV yields adiffraction-limited image of the object in the focal plane.Fig. 2B). Physically, this corresponds to a descan of the reflected light as performed in confocalmicroscopy.The D -matrix deduced from R is displayed in Fig. 1D1. Compared to R (Fig. 1C1), itexhibits long-range correlations both in the pupil (Fig. 1D3) and focal (Fig. 1D2) planes, re-spectively. On the one hand, by virtue of the van Cittert Zernike theorem ( ), the coherencelength d P of the distorted wave-field in the pupil plane scales as λf /δ in , with δ in being the spa-tial extension of the incoherent input focal spot | H in | (see Section S2). On the other hand,its correlation length d F in the focal plane corresponds to the size (cid:96) c of the isoplanatic patch(see Section S2). This is illustrated by examples of distorted wave-fields displayed in Fig. 1B3.While the original reflected wave-fronts did not display any similarity, the distorted component10isplays similar Fresnel rings whatever the focusing point r in . The D -matrix thus reveals in-put/output correlations of the wave-field that were originally completely hidden in the original R -matrix (Fig. 1C). Singular value decomposition of the distortion matrix
The next step is to extract and take advantage of those field-field correlations for imaging. Tothat aim, a singular value decomposition (SVD) of the distortion matrix is performed. It consistsin writing D as follows D = UΣV † (6)or, in terms of matrix coefficients, D ( u out , r in ) = N in (cid:88) p =1 σ p U p ( u out ) V ∗ p ( r in ) . (7) Σ is a diagonal matrix containing the real positive singular values σ i in a decreasing order σ > σ > .. > σ N in . 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 associatedwith the largest singular values and characterized by an important correlation between its linesand/or columns. In the D -matrix, these correlations are induced by the isoplanicity of theinput PSF H in . The single scattering and forward multiple scattering contributions are expectedto lie along the signal subspace since they exhibit a spatial invariance over each isoplanaticpatch ( ). On the contrary, the diffuse photons induced by the scattering layer ahead of thefocal plane give rise to a fully incoherent wave-field that will be equally distributed over all theeigenstates of D ( ). Hence, the pollution of the signal subspace by the multiple scatteringnoise scales as the inverse of the number of independent input focusing points mapping eachisoplanatic patch, ( (cid:96) c /δ ) . A large isoplanatic patch enables the SVD to drastically decrease11he multiple-to-single scattering ratio.To know which of the input or output correlations will dictate the SVD of D , relevantparameters are the numbers of independent speckle grains, M D and N D , exhibited by D at itsinput and output, respectively. The correlation degree of the distorted wave-field in each planeis actually inversely proportional to the corresponding number of independent speckle grains.In the focal plane, M D is given by the squared ratio between the FOI Ω and the coherence length d F of the distorted wave-field in the focal plane M D = (Ω /d F ) . (8) d F is the minimum between the isoplanatic length (cid:96) c and the characteristic fluctuation length (cid:96) γ of the object’s reflectivity (see Section S2). In the pupil plane, the number N D of independentspeckle grains scales as (see Section S1) N D = (cid:0) δ in /δ (cid:1) , (9)where δ is the diffraction-limit resolution at the output (Eq. 13). The domination of inputcorrelations implies the following condition: M D < N D . (10)If (cid:96) γ > (cid:96) c , the last equation can be translated as follows: The number M D = (Ω /(cid:96) c ) of IPssupported by the FOI should be smaller than the number N D of resolution cells that map eachinput focusing beam (Eq. 9). As we will see, this strong aberration condition is fulfilled in theexperiments presented in this work.When input correlations dominate, the effective rank of the signal subspace then corre-sponds to the number of independent spatial modes required to map the distorted wave-field inthe focal plane, i.e the number M D of IPs. As shown in Section S3, the SVD decomposes theFOI onto a set of orthonormal IMs defined by the input singular vectors V p . The corresponding12utput singular vectors U p yield the associated aberration phase laws in the pupil plane. Theircoherent 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, andsee how we can take advantage of it for deep imaging. π EF GH
PSF Full-field -π B
15 µm b e f o r e c o rr e c t i o n a f t e r c o rr e c t i o n F F R (r out ,0)R (r out ,0) C A arg{U }U D NA=0.1 V | || | Figure 3:
Imaging through a thick layer of rat intestinal tissue . ( A ) Experimental configura-tion. ( B , C ) Modulus of the first input singular vector V of D in the focal plane. ( D ) Modulusand phase of the first output singular vector U in the pupil plane. ( E ) Example of PSF deducedfrom the central column ( r in = ) of the raw focused matrix R . ( F ) Corresponding correctedPSF deduced from the central column of the focused matrix R (Eq. 12). ( G , H ) Comparisonof the full-field images F and F (Eq. 14) before and after aberration correction. Photo Credit:Amaury Badon, CNRS. Imaging over a single isoplanatic patch
The reflection and distortion matrices corresponding to the imaging experiment through a thicklayer of rat intestine are shown in Figs. 1C1 and D1, respectively. The long-range spatial cor-relations exhibited by D (Fig.1D2) seem to indicate that the isoplanatic hypothesis is close tobeing fulfilled in this experiment. The SVD of D confirms this intuition by exhibiting a predom-inant eigenstate. The corresponding singular vectors V and U are displayed in Fig. 3. Themodulus of V displays a contrasted image of the resolution target (Fig. 3B) but its resolution13s limited by the low spatial sampling of the illumination scheme. The output singular vector U corresponds to the wave-front induced by a virtual coherent reflector of scattering distri-bution | H in ( r − r in ) | , hence located on the optical axis in the focal plane (see Fig. 2C). Thisvirtual scatterer results from a coherent summation of the de-scanned input focal spots throughthe SVD process (see Section S3). Its phase is made of Fresnel rings mainly induced by theirregular 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 makethis virtual scatterer point-like and retrieve a diffraction-limited image (Fig. 2D), a normalizedvector ˜U should be considered, such that ˜ U ( u out ) = U ( u out ) / | U ( u out ) | . ˜U can be used tobuild an estimator ˆT of the transmission matrix between the pupil and focal planes, such thatits coefficients read ˆ T p ( u out , r in ) = ˜ U p ( u out ) T ( u out , r in ) (11)with p = 1 in the present case. This estimator can be used to project the recorded matrix R inthe focal basis both at input and output, such that R p = ˆT † p R (12)The coefficients R ( 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 yields the PSF of the imaging system at the input focusing point r in . The PSF for an inputfocusing point on the optical axis ( r in = ) is displayed in Fig. 3F. For sake of comparison,the corresponding initial focal spot is displayed in Fig. 3E. The latter one is extracted from thefocused matrix R deduced from R using T : R = T † R . While the initial PSF exhibits arandom speckle pattern (Fig. 3E), the PSF after correction shows a nearly diffraction-limitedfocal 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 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 δ = λ out , (13)with NA out = D out / (2 f ) = 0 . being the output numerical aperture. The numerical applicationof this formula yields δ (cid:39) . µ m in our experimental configuration. The mismatch between δ out and δ comes from the noisy aspect of U at large spatial frequencies (see Fig. 3D), whichprevents from an efficient aberration compensation over the whole numerical aperture.If the spatial sampling was equivalent at input and output, a confocal image could be de-duced from the diagonal elements ( r in = r out ) of R and R ( ). Here, as a sparse illuminationscheme has been employed ( δr in > δ ), a full-field image is considered and obtained by sum-ming R over its input elements: F p ( r out ) = (cid:88) r in | R p ( r out , r in ) | (14)with p = 0 or here. The corresponding images F and F are displayed in Figs. 3G 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 gainin image quality, the Strehl ratio is a relevant parameter ( ). It is defined as the ratio of thePSF peak intensity with and without aberration. Equivalently, it can also be defined in the pupilplane as the squared magnitude of the mean aberration phase factor. Its initial value S can thusbe directly derived from the D -matrix coefficients: S = |(cid:104) exp ( j arg { D ( u out , r in ) V ( r in ) } ) (cid:105)| (15)where the symbol (cid:104)· · · (cid:105) denotes an average over u out and r in . In the present case, the originalStrehl ratio is S = 8 × − . This experiment corresponds to imaging conditions far from15eing in the range of operation of conventional AO and explains why the patterns of the reso-lution target are so hardly visible on the raw image (Fig. 3G). The Strehl ratio S after the U correction can be directly extracted from the SVD of D (Eq. 7): S = |(cid:104) exp ( j arg { U ∗ ( u out ) D ( u out , r in ) V ( r in ) } ) (cid:105)| (16)The D -matrix correction leads to a Strehl ratio S = 3 × − . However, Eq. 16 gives thesame weight to bright and dark areas of the resolution target in the focal and pupil planes. Onepossibility is to consider a weighted average instead of Eq.16 by the object reflectivity | V ( r in ) | .This weighted Strehl ratio S (cid:48) then reaches the value of . × − . Such a Strehl ratio valueis relatively low but it should be kept in mind that the distortion matrix is associated with aPSF in reflection that convolves the transmit and receive PSFs. Our measurement of the Strehlratio is thus degraded by: ( i ) the subsistence of input aberrations; ( ii ) the presence of a diffusemultiple scattering background that acts here as an additive noise. Note, however, that the gainin terms of Strehl ratio is absolute; this is the relevant quantity to assess the benefit of our matrixapproach. This gain here is spectacular ( S (cid:48) / S (cid:48) ∼ D -matrix correction (see Fig. 3H). Figure S3 shows howthis drastic improvement of the Strehl ratio allows us to push back the imaging depth limit from450 µ m to almost 1 mm.This first experiment demonstrates the benefit of the D -matrix in the simple case of a FOIcontaining a single IP. In the next section, the case of multiple IPs is tackled. Imaging over multiple isoplanatic patches
The first element of the group 6 in the resolution target is now imaged through an aberratinglayer consisting in a rough plastic sheet placed d = R is measured over a FOI of × µ m with N in =441 input wave-fronts, a spatial sampling δr in = 2 . µ m and an input pupil aperture16
10 20 30 40 5002468 10 -3 Image a. F u ll f i e l d PSF M a t r i x b. c.e.d. f. s i n g u l a r v a l u e s π 10 -π U P h a s e m a s k C o n f o c a l i m a g e π U U U U
20 µm20 µm g.h. N o r m a li z e d NA=0.1
Figure 4:
Matrix imaging over multiple isoplanatic patches . ( A ) Schematic of the exper-iment. A resolution target (USAF 1951) is positioned at a distance d = 1 mm underneath arough plastic film (see inset). ( B ) Original full-field image F (Eq.14). ( C ) Example of PSF de-duced from a column of the raw focused matrix R . ( D ) Plot of the normalized singular values ˜ σ i of D . The red circles correspond to the eight first singular values (signal subspace), while thenoisy singular values are displayed in blue. ( E ) Matrix image constructed from the eight firsteigenstates of D (Eq. 20). ( F ) Example of PSF deduced from a column of the corrected focusedmatrix R . ( G ) Phase of the four first singular vectors U p . ( H ) Confocal images deduced fromthe focused reflection matrices R p [Eq.18]. Photo Credit: Amaury Badon, CNRS.17 in × D in = 1 . × . mm . At the output, the wave-field is recorded over a pupil size of D out × D out = 2 × mm with N out = 12321 pixels and a spatial sampling δu out =18 µ m.The full-field image F (Eq. 14) and an example of PSF (Eq. 12) are displayed in Figs. 4Aand B, respectively. The PSF is strongly degraded with a characteristic focal spot dimension δ out ∼ µ m. This PSF dimension allows an estimation of the coherence length (cid:96) c of theaberrating layer. Indeed, under a thin phase screen model ( ), the IP dimension (cid:96) c coincideswith the coherence length of the aberration transmittance. It turns out that the PSF width isinversely proportional to (cid:96) c in this experiment: δ out ∼ λd/(cid:96) c . The IP size and the number of IPssupported by the FOI can be deduced from the PSF width δ : (cid:96) c ∼ µ m and M D ∼ (Ω /(cid:96) c ) ∼ . A D -matrix is deduced from R (Eq. 4). Its analysis leads to the following estimation of theinitial Strehl ratio: S (cid:48) = 1 . × − (Eq. 15). This particularly strong aberration level accountsfor the highly blurred aspect of the full-field image in Fig. 4A. This experimental situationis thus particularly extreme, even almost hopeless, for a successful imaging of the resolutiontarget. Yet the SVD of D will provide the solution.Fig. 4D displays the histogram of the normalized singular values ˜ σ i = σ i / (cid:80) N in j =1 σ j . Ifrecorded data was not corrupted by experimental noise, the matrix would be of effective rank M D . We could use all the eigenstates of D associated with non-zero singular values to retrievean image of the object. In Fig. 4D, only few singular values seem to emerge from the noisebackground. Hence, it is difficult to determine the number of eigenstates we need to consider toproperly reconstruct an image of the object. This issue can be circumvented by computing theShannon entropy H of the singular values (
40, 41 ), such that H (˜ σ i ) = − N in (cid:88) i =1 ˜ σ i log (˜ σ i ) . (17)Shannon entropy delivers the maximally-noncommittal data set at a given signal-to-noise ratio,18hat is to say, the most information with the least artifact. The Shannon entropy can be usedas an indicator of how many eigenstates are needed to build an adequate image of the objectwithout being affected by experimental noise.The singular values of Fig. 4D have an entropy H (cid:39) . . Hence, only the eight first singularstates shall be considered. Fig. 4G displays the phase of the four first singular vectors U i inthe pupil plane. They display phase distortions whose typical coherence length u c scales as f (cid:96) c /d ∼ µ m. The phase conjugation of these singular vectors should compensate for thedetrimental effect of the aberrating layer in different parts of the FOI. A set of focused reflectionmatrices R p can be deduced (Eq.12). Fig. 4f displays an example of corrected PSF extractedfrom a column of R . Its comparison with the original PSF in Fig. 4C shows how the phaseconjugation of U allows one to compensate for the aberrations at this incident focusing point.On the one hand, the PSF width is narrowed by a factor 20, going from δ out ∼ µ m to . µ m. The latter value should be compared with the diffraction-limited resolution δ ∼ µ m(Eq. 13) in our experimental conditions. The Strehl ratio is increased by a factor 2.2 × toreach the final value S (cid:48) = 3 . × − (Eq. 16). Again, this value of the Strehl ratio is probablyunderestimated because of input aberrations and multiple scattering.Confocal images can be deduced from the focused reflection matrices R p : I p ( r out ) = (cid:88) r in | R p ( r out , r in ) | e −|| r out − r in || / l p (18)where l p is the aperture of the numerical confocal pinhole ( ). This finite aperture enables anaverage of the output image over neighbour incident focusing points in order to smooth out thesparse illuminations. Fig. 4H displays the confocal images I p for l p = 2 µ m. For a specularobject such as a resolution target, the SVD has indeed the property of decomposing into a setof orthogonal IMs of spatial period (cid:96) c (see Section S3). Their shape depends on the auto-correlation function of the aberrating phase screen. A general trend is that the spatial frequency19ontent of the eigenvectors increases with their rank. If this function presents an exponential orsinc model, the FOI will be spatially decomposed into sinusoidal wave functions ( ) analogousto optical fiber modes or to prolate spheroidal wave functions ( ), respectively. Here, theautocorrelation function of the aberrating phase displays a Gaussian-like shape. The FOI is thusspatially mapped onto Hermite-Gaussian wave functions analogous to laser cavity modes ( ).The normalized pupil singular vectors ˜U p yield a set of orthogonal phase transmittances thatmap aberrations onto each isoplanatic mode. A coherent combination of these singular vectorsshould lead, in principle, to a satisfying estimator of the transmission matrix (see Section S3) ˆT = H (˜ σ i ) (cid:88) p =1 ˜U p ◦ T . (19)In practice, a final image I of the resolution target can be obtained by summing the previousIMs I p : I ( r out ) = H (˜ σ i ) (cid:88) p =1 I p ( r out ) . (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 anda contrasted image of the resolution target is recovered over the whole FOI. This experimentdemonstrates how the D -matrix enables a decomposition of the FOI into several IMs and amapping of each of them onto orthonormal distorted phase laws. However, this demonstrationhas been restricted to the case of a 2D aberrating phase layer. In the next section, we consider thecase of a cornea with deteriorated transparency as a three-dimensional aberrating and scatteringstructure. Imaging through a hazy cornea
The experimental configuration is displayed in Fig. 5A. The number “3” of the group 5 in theresolution target is imaged through a 700- µ m-thick edematous non-human primate cornea. The20
50 100012345 U U I confocal matrix edematouscornea U I AB CD EF I s i n g u l a r v a l u e s N o r m a li z e d NA=0.110 µm
Figure 5:
Imaging through corneal tissue with deteriorated transparency . ( A ) Schematicof the experiment. A resolution target (USAF 1951) is positioned below an edematous non-human primate cornea (see inset). ( B ) Plot of the normalized singular values ˜ σ i of D . The redcircles correspond to the eleven first singular values (signal subspace), while the noisy singularvalues are displayed in blue. ( C ) Original confocal image deduced from the focused reflectionmatrix R (Eq. 18). ( D ) Final matrix image constructed from the eleven first eigenstates of D (Eq. 20). ( E ) Real parts of U , U and U . ( F ) Corresponding confocal images deduced fromthe focused reflection matrices R p (Eq. 18). 21eflection matrix R is measured over a FOI × µ m by means of N in =625 input wave-fronts, a spatial sampling δr in = 2 . µ m and an input pupil aperture D in × D in = 1 × mm . Atthe output, the wave-field is recorded over an output pupil size D out × D out = 2 × mm with N out = 1296 pixels and a spatial sampling length δu out =56 µ m. Fig. 5C displays the confocalimage I deduced from R with l p = 1 µ m (Eq. 18). Multiple scattering and aberrationsinduced by the cornea induce a random speckle reflected wave-field that prevents from imagingthe resolution target. On the contrary, as we will see, the D -matrix analysis allows us to nicelyrecover the pattern “3” of the resolution target (see Fig. 5D).Fig. 5C displays the spectrum of the singular values ˜ σ i . The first singular value emergesfrom the rest of the spectrum but it is difficult to know until which rank the eigenstates canbe considered as belonging to the signal subspace. As previously, the Shannon entropy of thesingular values yields an unambiguous answer: H ( σ i ) = 10 . . The 11 th first singular statesshould thus be considered. Fig. 5E displays the 1 st , 6 th and 11 th singular vectors U i in the pupilplane. The complexity of the wave-front distortion, i.e their spatial frequency content, increaseswith the rank of the corresponding singular values. The corresponding IMs I p (Eq. 18) aredisplayed in Fig. 5F. While the first singular vector U allows a wide-field correction of low-order aberration, the higher rank singular vectors are associated with high-order aberrations thatare effective over IMs of smaller dimension. In Fig. 5D, the whole spatial frequency spectrumof wave-front distortions is finally compensated by smartly combining the confocal images I p associated with each singular state from D ’s signal subspace (Eq.20). The comparison of theinitial (Fig. 5C) and final (Fig. 5D) images is spectacular with a Strehl ratio gain S (cid:48) /S (cid:48) = 230 .The comparison of I (Fig. 5D) and I (see the first inset of Fig. 5F) illustrates the benefit ofa matrix approach of aberration correction compared to conventional AO, since the latter onewould yield I in theory.This decomposition of complex aberration phase laws over a set of IMs opens important22erspectives for aberrometry. It actually goes well beyond state-of-the-art techniques that basi-cally consist in a simple projection over a set of Zernike polynomials. Moreover, an estimatorof the single-to-multiple scattering (SMR) ratio can be built on the relative energy between thesignal and noise sub-spaces of D : SMR = (cid:80) H ( σ i ) i =1 σ i (cid:80) N in i = H ( σ i )+1 σ i . (21)The SMR can actually be a quantitative bio-marker of the corneal opacification or a quantitativemeasure of corneal transparency ( ). Based on a fit with a recent analytical study of theSMR ( ), the cornea thickness L can be estimated in terms of scattering mean free path (cid:96) s : L ∼ (cid:96) s (see Fig. S3). As the corneal thickness is known ( L = 700 µ m), the scattering meanfree path can be deduced: (cid:96) s ∼ µ m. Interestingly, this value is in excellent agreement withrecent ex-vivo measurements of (cid:96) s in pathological corneas with compromised transparency ( ).The value of 9 (cid:96) s highlights the difficult experimental conditions under which the imaging of theresolution target has been successfully achieved.In conclusion, this last experiment shows the potential of a matrix approach for eye aber-rometry and turbidimetry, such as for improved quality control of donor tissue assessment priorto corneal transplantation ( ). Of course, this method is by no means limited to ophtalmicapplications. It can be applied to the characterization of any kind of biological tissues providedthat we have access to the associated reflection matrix. DISCUSSION
In this article, we present a novel and non-invasive method for aberration compensation anddiffraction-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 distortedcomponent of the reflected the wave-field in the pupil plane. This operator connecting position23nd spatial frequency has some analogy with the Wigner distribution function ( ). However,the Wigner transform applies to a single variable of a function, i.e to a single vector in a discreteformalism. Here, our position-momentum analysis is performed between the input and outputof a reflection matrix.The concept of distortion matrix is to measure the back-scattered waves in a de-scannedframe while scanning the sample with focused illuminations. This approach has some similar-ity with a previous AO approach ( ) in its hardware configuration. The main difference is that,in this study, wave-fronts are averaged by the Shack-Hartmann type of analysis and this AOapproach 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 averagedover the whole FOI, we propose to study the spatial and angular correlations of the distortionoperator through an SVD. In this manner, we demonstrate the efficient compensation of bothlow- and high-order aberrations over multiple IPs. Moreover, our approach relies on the Shan-non entropy that provides an objective criterion to determine the number of IPs supported bythe FOI. This is in contrast with recent works based on a far-field reflection matrix in which theFOI was arbitrarily divided into sub-areas 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 spatiallyincoherent. It thus mainly lies along the noise subspace of the D -matrix. Thanks to thesefeatures, we were able to image through almost 10 scattering mean free paths of biologicaltissues, which is beyond the imaging depth of conventional OCT systems for such specimens(see Fig. S3). Compared to the previously developed smart-OCT method that was able to detectfew bright scatterers at large penetration depth (12 (cid:96) s ) ( ), the D -matrix approach yields adirect image of the sample reflectivity at a diffraction-limited resolution. Additionally, ourapproach enables to quantitatively estimate the amount of multiply-scattered light. Combined24ith a conventional image, this parameter is of importance for characterization purposes.The distortion operator thus opens a new route towards real-time deep optical imaging ofbiological tissues. In that respect, the experimental set up and procedure used in this paper areclearly perfectible. While post-process operations take less than one minute on a regular laptop,the main limitation in the current experimental configuration is the acquisition time. In particu-lar, the scanning illumination scheme was not optimized because of the SLM speed. While theuse of a galvanometric mirror or a high-speed deformable mirror would reduce drastically theacquisition time at the cost of a more complex setup, we counteracted this issue with a sparseillumination. However, this, in return, limited the available number of angular degrees of free-dom at the input, which prevents us from an aberration correction of the incident wave-field. Byoptimizing the experimental apparatus and acquisition scheme, large reflection matrices can bemeasured in a few seconds. For instance, Yoon et al. recently demonstrated the acquisition of a10000 modes matrix in 15 seconds with the same degree of control for the incident and reflectedwaves ( ). In that case, a simultaneous correction of aberrations at the input and output is ab-solutely possible under the distortion matrix approach by alternatively projecting the incidentand reflected wave-fields in the focal and pupil planes. In view of 3D imaging, our approachcan also be coupled to computational AO techniques ( ) in order to tackle depth-dependentaberrations and restore a diffraction-limited resolution in all directions. An alternative is toswitch from a scanning to a full-field illumination scheme. A measurement of the coherentreflection matrix R can be performed under a spatially incoherent illumination (
49, 50 ). Thisfull-field configuration would allow to record the reflection matrix over millimetric volumes ina moderate acquisition time.Finally, we used a negative resolution target as the sample to be imaged in this work. Thereason is that this highly contrasted object was the ideal specimen to clearly highlight the issueof multiple isoplanatic areas. Beyond the proof-of-concept experiments presented in this arti-25le, a direct imaging of biological specimens over large penetration depth will be the next step.Interestingly, the assumption on which our method is based (Eq. 10) can easily be met in biolog-ical 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, thedistortion matrix approach can still work but the FOI has to be beforehand sub-divided intoindividual IPs (
47, 51 ). The ability of identifying multiple IPs will also be particularly promis-ing to map the specimen-induced aberration and the single-to-multiple scattering ratio. Asidefrom aberrometry and/or turbidimetry, future in-vivo implementations of our approach haveimplications beyond that of ocular media characterization, most notably for imaging throughnon-transparent ocular media (e.g., retinal imaging through a turbid cornea or through cataractopacities) ( ).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 themismatch between the phase of the recorded reflection matrix and those of a reference matrixthat would be obtained in an ideal configuration. As shown in this paper, D gives access tothe non-invasive transmission matrix between each sensor and each voxel of the FOI. Then, bysolving the corresponding inverse problem, an image of a scattering sample can be obtained asif the medium ahead was made transparent. The D -matrix concept is very general. It can beextended to any kind of waves and experimental configurations for which a measurement of theamplitude and phase of the reflected wave-field is possible under multiple illuminations ( ).A recent work actually demonstrates the benefits of this concept for ultrasound imaging in arandom scattering regime ( ). This D -matrix concept thus opens a new route towards a globaland non-invasive matrix approach of deep imaging in biological tissues.26 ATERIALS AND METHODS
Experimental set up
The experimental configuration is identical to the one described in ( ) except for the MO thathad been replaced by a water immersion one. The following components were used in theexperimental setup (see Fig. S1): a femtosecond laser (Femtosecond Fusion 20-400, centralwavelength: 810 nm, bandwidth: 40 nm), an SLM (PLUTONIR-2, HOLOEYE), an objectivelens (40 × ; NA, 0.8; Nikon), and a CCD camera (Dalsa Pantera 1M60) with a dynamic range of60 dB. The incident light power in the back pupil plane of the MO was 10 mW in the experiment.Thus, the radiant flux was W/cm at the focal spot in free space. For each input wave-front, the complex-reflected wave field was extracted from four intensity measurements usingphase shifting interferometry. The acquisition time of the reflection matrix was approximately2 minutes. Image acquisition and data analysis
Both data acquisition and analysis were performed using Matlab custom-written codes. Thesecodes are available from the authors upon request.
SUPPLEMENTARY MATERIALS
Section S1. Correlations of the reflected and distorted wave-fields in the pupil plane.Section S2. Correlations of the reflected and distorted wave-fields in the focal plane.Section S3. Singular value decomposition of the distortion matrix.Fig. S1. Measuring the time-gated reflection matrix.Fig. S2. Conjugating the pupil, focal and imaging planes.Fig. S3. Predicting the single-to-multiple scattering ratio in biological tissues.27ig. S4. Building the reflection matrix R .Fig. S5. Modeling light propagation from the virtual source plane to the output pupil plane.Tab. S1. Glossary of the variables used in the study.Tab. S2. Glossary of the matrices used in the study.References ( ). REFERENCES AND NOTESReferences
1. H. W. Babcock, The possibility of compensating astronomical seeing.
Publications of theAstronomical Society of the Pacific , 229–236 (1953).2. R. Foy, A. Labeyrie, Feasibility of adaptive telescope with laser probe. Astron. Astrophys. , L29–L31 (1985).3. S. T. Thurman, J. R. Fienup, Correction of anisoplanatic phase errors in digital holography.
J. Opt. Soc. Am. A , 995–999 (2008).4. A. E. Tippie, J. R. Fienup, Multiple-plane anisoplanatic phase correction in a laboratorydigital holography experiment. Opt. Lett. , 3291–3293 (2010).5. M. J. Booth, M. A. Neil, R. Juˇskaitis, T. Wilson, Adaptive aberration correction in a confo-cal microscope. Proc. Natl. Acad. Sci. U. S. A. , 5788–5792 (2002).6. X. Tao, B. Fernandez, O. Azucena, M. Fu, D. Garcia, Y. Zuo, D. C. Chen, J. Kubby,Adaptive optics confocal microscopy using direct wavefront sensing. Opt. Lett. , 1062–1064 (2011). 28. D. D´ebarre, E. J. Botcherby, T. Watanabe, S. Srinivas, M. J. Booth, T. Wilson, Image-basedadaptive optics for two-photon microscopy. Opt. Lett. , 2495–2497 (2009).8. N. Ji, D. E. Milkie, E. Betzig, Adaptive optics via pupil segmentation for high-resolutionimaging in biological tissues. Nat. Methods , 141–147 (2010).9. I. N. Papadopoulos, J.-S. Jouhanneau, J. F. Poulet, B. Judkewitz, Scattering compensationby focus scanning holographic aberration probing (f-sharp). Nat. Photon. , 116–123(2017).10. M. Rueckel, J. A. Mack-Bucher, W. Denk, Adaptive wavefront correction in two-photonmicroscopy using coherence-gated wavefront sensing. Proc. Nat. Acad. Sci. U. S. A. ,17137–17142 (2006).11. B. Hermann, E. Fern´andez, A. Unterhuber, H. Sattmann, A. Fercher, W. Drexler, P. Prieto,P. Artal, Adaptive-optics ultrahigh-resolution optical coherence tomography.
Opt. Lett. ,2142–2144 (2004).12. S. G. Adie, B. W. Graf, A. Ahmad, P. S. Carney, S. A. Boppart, Computational adaptiveoptics for broadband optical interferometric tomography of biological tissue. Proc. Natl.Acad. Sci. U. S. A. , 7175–7180 (2012).13. B. Judkewitz, R. Horstmeyer, I. M. Vellekoop, I. N. Papadopoulos, C. Yang, Translationcorrelations in anisotropically scattering media.
Nat. Phys. , 684-689 (2015).14. F. J. Rigaut, B. L. Ellerbroek, R. Flicker, Adaptive Optical Systems Technology , P. L. Wiz-inowich, ed. (SPIE, 2000).15. Z. Kam, P. Kner, D. Agard, J. W. Sedat, Modelling the application of adaptive optics towide-field microscope live imaging.
J. Microscopy , 33–42 (2007).296. R. D. Simmonds, M. J. Booth, Modelling of multi-conjugate adaptive optics for spatiallyvariant aberrations in microscopy.
J. Opt. , 094010 (2013).17. M. Minsky, Confocal scanning microscope (1955).18. M. R. Hee, E. A. Swanson, J. G. Fujimoto, J. A. Izatt, Femtosecond transilluminationtomography in thick tissues. Opt. Lett. , 1107–1109 (1993).19. I. M. Vellekoop, A. P. Mosk, Focusing coherent light through opaque strongly scatteringmedia. Opt. Lett. , 2309–2311 (2007).20. S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, S. Gigan, Measuring theTransmission Matrix in Optics: An Approach to the Study and Control of Light Propagationin Disordered Media. Phys. Rev. Lett. , 100601 (2010).21. S. M. Popoff, G. Lerosey, M. Fink, A. C. Boccara, S. Gigan, Image transmission throughan opaque material.
Nat. Commun. , 1–5 (2010).22. M. Kim, Y. Choi, C. Yoon, W. Choi, J. Kim, Q.-H. Park, W. Choi, Maximal energy transportthrough disordered media with the implementation of transmission eigenchannels. NaturePhoton. , 583-587 (2012).23. T. Cizmar, K. Dholakia, Exploiting multimode waveguides for pure fibre- based imaging. Nat. Commun. (2012).24. I. N. Papadopoulos, S. Farahi, C. Moser, D. Psaltis, Focusing and scanning light througha multimode optical fiber using digital phase conjugation. Opt. Express , 10583-10590(2012).25. S. A. Alexandrov, T. R. Hillman, T. Gutzler, D. D. Sampson, Synthetic aperture fourierholographic optical microscopy. Phys. Rev. Lett. , 168102 (2006).306. S. M. Popoff, A. Aubry, G. Lerosey, M. Fink, A. C. Boccara, S. Gigan, Exploiting the time-reversal operator for adaptive optics, selective focusing, and scattering pattern analysis. Phys. Rev. Lett. , 263901 (2011).27. A. Badon, D. Li, G. Lerosey, A. C. Boccara, M. Fink, A. Aubry, Smart optical coherencetomography for ultra-deep imaging through highly scattering media.
Sci. Adv. , e1600370(2016).28. Y. Choi, T. R. Hillman, W. Choi, N. Lue, R. R. Dasari, P. T. C. So, W. Choi, Z. Yaqoob,Measurement of the time-resolved reflection matrix for enhancing light energy delivery intoa scattering medium. Phys. Rev. Lett. , 243901 (2013).29. S. Jeong, Y.-R. Lee, W. Choi, S. Kang, J. H. Hong, J.-S. Park, Y.-S. Lim, H.-G. Park,W. Choi, Focusing of light energy inside a scattering medium by controlling the time-gatedmultiple light scattering.
Nat. Photon. , 277–283 (2018).30. S. Kang, S. Jeong, W. Choi, H. Ko, T. D. Yang, J. H. Joo, J.-S. Lee, Y.-S. Lim, Q.-H.Park, W. Choi, Imaging deep within a scattering medium using collective accumulation ofsingle-scattered waves. Nature Photon. , 253–258 (2015).31. S. Kang, P. Kang, S. Jeong, Y. Kwon, T. D. Yang, J. H. Hong, M. Kim, K.-D. Song, J. H.Park, J. H. Lee, et al. , High-resolution adaptive optical imaging within thick scatteringmedia using closed-loop accumulation of single scattering. Nat. Commun. , 2157 (2017).32. J.-L. Robert, M. Fink, Green’s function estimation in speckle using the decomposition ofthe time reversal operator: Application to aberration correction in medical imaging. J.Acoust. Soc. Am. , 866-877 (2008).33. J.-L. Robert, M. Fink, The time-reversal operator with virtual transducers: Application tofar-field aberration correction.
J. Acoust. Soc. Am. , 3659-3668 (2008).314. G. Osnabrugge, R. Horstmeyer, I. N. Papadopoulos, B. Judkewitz, I. M. Vellekoop, Gener-alized optical memory effect.
Optica , 886–892 (2017).35. S. Jacques, Optical properties of biological tissues: a review. Phys. Med. Biol. , R37-R61(2013).36. M. Born, E. Wolf, Principles of optics: electromagnetic theory of propagation, interferenceand diffraction of light (CUP Archive, 1999).37. J. Mertz, H. Paudel, T. G. Bifano, Field of view advantage of conjugate adaptive optics inmicroscopy applications.
Appl. Opt. , 3498–3506 (2015).38. A. Badon, A. C. Boccara, G. Lerosey, M. Fink, A. Aubry, Multiple scattering limit inoptical microscopy. Opt. Express , 28914–28934 (2017).39. V. N. Mahajan, Strehl ratio for primary aberrations: some analytical results for circular andannular pupils. J. Opt. Soc. Am. , 1258–1266 (1982).40. L. L. Campbell, Minimum coefficient rate for stationary random processes. Inf. Control ,360–371 (1960).41. S. J. Roberts, W. Penny, L. Rezek, Temporal and spatial complexity measures for elec-troencephalogram based brain-computer interfacing. Med. Biol. Eng. Comput. , 93–98(1999).42. R. G. Ghanem, R. D. Spanos, Stochastic Finite Elements: A Spectral Approach (Springer-Verlag, New York, 1991).43. J.-L. Robert, M. Fink, The prolate spheroidal wave functions as invariants of the time re-versal operator for an extended scatterer in the fraunhofer approximation.
J. Acoust. Soc.Am. , 218-226 (2009). 324. A. Aubry, J. de Rosny, J.-G. Minonzio, C. Prada, M. Fink, Gaussian beams and legendrepolynomials as invariants of the time reversal operator for a large rigid cylinder.
J. Acoust.Soc. Am. , 2746-2754 (2006).45. R. Bocheux, P. Pernot, V. Borderie, K. Plamann, K. Irsch, Quantitative measures of cornealtransparency, derived from objective analysis of depth-resolved corneal images, demon-strated with full-field optical coherence tomographic microscopy.
PLOS ONE , 1-10(2019).46. M. Bastiaans, Wigner distribution function and its application to first-order optics. J. Opt.Soc. Am. , 1710–1716 (1979).47. M. Kim, Y. Jo, J. H. Hong, S. Kim, S. Yoon, K.-D. Song, S. Kang, B. Lee, G. H. Kim,H.-C. Park, et al. , Label-free neuroimaging in vivo using synchronous angular scanning mi-croscopy with single-scattering accumulation algorithm. Nature Commun. , 3152 (2019).48. S. Yoon, H. Lee, J. H. Hong, Y.-S. Lim, W. Choi, Laser scanning reflection-matrixmicroscopy for label-free in vivo imaging of a mouse brain through an intact skull. arXiv:1910.04681 (2019).49. A. Badon, G. Lerosey, A. C. Boccara, M. Fink, A. Aubry, Retrieving time-dependentgreen’s functions in optics with low-coherence interferometry. Phys. Rev. Lett. , 023901(2015).50. A. Badon, D. Li, G. Lerosey, A. C. Boccara, M. Fink, A. Aubry, Spatio-temporal imagingof light transport in highly scattering media under white light illumination.
Optica , 1160–1166 (2016). 331. W. Lambert, L. A. Cobus, T. Frappart, M. Fink, A. Aubry, Distortion matrix approachfor ultrasound imaging of random scattering media. Proc. Natl. Acad. Sci. U. S. A. ,14645-14656 (2020).52. A. Badon, V. Barolle, A. C. Boccara, M. Fink, K. Irsch, A. Aubry, Towards characteri-zation and compensation of loss of anterior segment transparency.
IEEE Transactions onEngineering in Medicine and Biology Conference (EMBC) (2019).53. V. Lauer, New approach to optical diffraction tomography yielding a vector equation ofdiffraction tomography and a novel tomographic microscope.
Journal of Microscopy ,165-176 (2002).54. S. Shahjahan, A. Aubry, F. Rupin, B. Chassignole, A. Derode, A random matrix approachto detect defects in a strongly scattering polycrystal: How the memory effect can helpovercome multiple scattering.
Appl. Phys. Lett. , 234105 (2014).55. T. Blondel, J. Chaput, A. Derode, M. Campillo, A. Aubry, Matrix approach of seismicimaging: Application to the erebus volcano, Antarctica.
J. Geophys. Res.: Solid Earth ,10936–10950 (2018).56. T. Zhang, K. Unger, G. Maire, P. C. Chaumet, A. Talneau, C. Godhavarti, H. Giovannini,K. Belkebir, A. Sentenac, Multi-wavelength multi-angle reflection tomography.
Opt. Ex-press , 26093–26105 (2018).57. J.-L. Robert, Evaluation of green’s functions in complex media by decomposition of thetime reversal operator: Application to medical imaging and aberration correction, Ph.D.thesis, Universite Paris 7 - Denis Diderot (2007).58. M. Priestley, Spectral analysis and time series (Academic Press, London, 1988).349. J. W. Goodman,
Statistical Optics (Wiley, New York, 2000).60. C. Prada, J.-L. Thomas, Experimental subwavelength localization of scatterers by decom-position of the time reversal operator interpreted as a covariance matrix.
J. Acoust. Soc.Am. , 235–243 (2003).61. C. Prada, M. Fink, Eigenmodes of the time reversal operator: A solution to selective focus-ing in multiple-target media.
Wave Motion , 151–163 (1994). Acknowledgements
The authors wish to thank Laura Cobus, William Lambert, Paul Balondrade and Serge Meimonfor fruitful discussions.
Funding:
The authors are grateful for the funding provided by LabexWIFI (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 fromthe French “Direction G´en´erale de l’Armement” (DGA). This project has received funding fromthe European Research Council (ERC) under the European Union’s Horizon 2020 research andinnovation programme (grant agreements nos. 610110 and 819261, HELMHOLTZ* and REM-INISCENCE projects, respectively). K.I. acknowledges financial support from the EuropeanUnion’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curiegrant agreement No. 709104.
Author Contributions
A.A. initiated and supervised the project. A.B. built the experimentalsetup and performed the experiments. K.I., A.C.B and M.F. initiated the ophthalmic applica-tion. 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. andA.A. prepared the manuscript. A.B., V.B., K.I., A.C.B., M.F., and A.A. discussed the resultsand contributed to finalizing the manuscript. 35 ompeting Interests:
A.A., M.F., A.C.B, A.B. and V.B. are inventors on a patent related tothis work held by CNRS (no. WO2020016249, published January 2020). All authors declarethat they have no other competing interests.
Data and materials availability:
All data needed to evaluate the conclusions in the paper arepresent in the paper and/or the Supplementary Materials. Additional data related to this papermay be requested from the authors. 36 upplementary MaterialsS1 Correlations of the reflected and distorted wave-fields inthe pupil plane
In this section, we derive the pupil correlations of the R - and D - matrices. Our aim is to providea theoretical proof of the experimental observation made in Fig.1C2 and 1D2. The distortedwave-fields exhibit correlations over a longer range than the reflected wave-fields in the pupilplane. For sake of analytical tractability but without loss of generality, we will assume in thissection: ( i ) a set of fully incoherent input focal spots ( i.e a strong aberration regime); ( ii ) a field-of-illumination (FOI) contained in a single IP. The main result is the following: While the pupilcorrelation length r P of R scales as the inverse of the FOI size ( r P ∼ λf / Ω ), the correlationlength d P of D is inversely proportional to the width δ in of the input PSF ( r P ∼ λf /δ in ). Theproofs of these two assertions are provided below. S1.1 Reflection matrix
To investigate the angular correlations of the reflected wave-field, the correlation matrix B R = N − in RR † should be considered. Using Eq. 2, its coefficients can be expressed as follows: B R ( u out , u (cid:48) out ) = N − in (cid:90) Ω d r (cid:90) Ω d r (cid:48) T ( u out , r ) T ∗ ( u (cid:48) out , r (cid:48) ) γ ( r ) γ ∗ ( r (cid:48) ) (cid:88) r in H in ( r , r in ) H ∗ in ( r (cid:48) , r in ) . (S1)In a strong aberration regime, the input focal spots can be considered as fully incoherent, (cid:104) H in ( r , r in ) H ∗ in ( r (cid:48) , r in ) (cid:105) = (cid:10) | H in ( r , r in ) | (cid:11) δ ( r − r (cid:48) ) . (S2)where δ is the Dirac distribution and the symbol (cid:104)· · · (cid:105) denotes an ensemble average. In astrong aberration regime, B R can be decomposed as the sum of a covariance matrix (cid:104) B R (cid:105) anda perturbation term δ B R : B R = (cid:104) B R (cid:105) + δ B R , (S3)1he correlation matrix B R (Eq. S3) should converge towards the covariance matrix (cid:104) B (cid:105) fora sufficiently large number M R ∼ (Ω /δ ) of independent speckle grains in the focal plane(Eq.S27). More precisely, the intensity of the perturbation term in Eq.S3, | δ B R | , scales as theinverse of M R ( ).Assuming the convergence of B R towards (cid:104) B R (cid:105) ( M R >> ), the correlation coefficients B R ( u out , u (cid:48) out ) (Eq. S1) can be expressed as follows: B R ( u out , u (cid:48) out ) = N − in (cid:90) d r T ( u out , r ) | γ ( r ) | T ∗ ( u (cid:48) out , r ) × (cid:88) r in (cid:10) | H in ( r , r in ) | (cid:11) , (S4)To go further, an isoplanatic configuration should be considered. On the one hand, this meansthat the input PSF is invariant by translation: H in ( r , r in ) = H in ( r − r in ) (S5)On the other hand, the output transmission matrix coefficients T ( u out , r ) can be decomposed asthe product of the transmittance ˆ H out ( u out ) of the aberrating layer and the free-space transmis-sion matrix coefficients T ( u out , r ) (Eq.3): T ( u out , r ) = ˆ H out ( u out ) T ( u out , r ) . (S6)Injecting these last equations and Eq. 3 into Eq. S4 leads to the following expression for : B R ( u out , u (cid:48) out ) = I in ˆ H out ( u out ) ˆ H ∗ out ( u (cid:48) out )ˆ γ ( u (cid:48) out − u out ) (S7)where I in = N − in (cid:88) r in (cid:10) | H in ( r − r in ) | (cid:11) is the mean input PSF intensity and ˆ γ ( u ) = (cid:90) d r | γ ( r ) | exp( − j π u . r /λf )
2s the 2D Fourier transform of the scattering distribution | γ ( r ) | in the focal plane. This quantity,which dictates the correlations displayed by R in the pupil plane, can be seen as an incoherentstructure factor of the object placed in the FOI. The corresponding coherence length r p scalesas r P ∼ λf / Ω , (S8)The number N R of independent speckle grains in the reflected wave-field is given by the squaredratio between the output pupil size D out = λf /δ and the pupil coherence length r P : N R ∼ (Ω /δ ) (S9) N R scales as the number of output resolution cells mapping the object.These theoretical predictions account for the incoherence of the reflected wave-field shownin Fig. 1C3. This figure plots the auto-correlation function B R (∆ u ) of the reflected wave-field inthe pupil plane. It is computed by averaging the correlation matrix coefficients B R ( u out , u (cid:48) out ) over couples ( u out , u (cid:48) out ) sharing the same relative position ∆ u = u out − u (cid:48) out . S1.2 Distortion matrix
As highlighted by Fig. 1C and demonstrated above, the reflection matrix displays a randomfeature at the output in the strong aberration regime. Now we will show how the realignment ofthe reflected wave-fronts in the pupil plane can reveal the angular correlations of the distortedcomponent.The distortion matrix D is defined as the Hadamard product between the reflection matrix R and the reference transmission matrix T ∗ (Eqs. 4-5). In the isoplanatic limit (Eqs. S5-S6)and using Eq.2, the D -matrix coefficients can be expressed as follows D ( u out , r in ) = ˆ H out ( u out ) (cid:90) d r T ( u out , r − r in ) γ ( r ) H in ( r − r in ) . (S10)3o investigate the angular correlations between distorted wave-fields, the spatial correlationmatrix B D = N − in DD † is investigated. Its coefficients can be expressed as follows: B D ( u out , u (cid:48) out ) = N − in ˆ H ( u out ) ˆ H ∗ ( u (cid:48) out ) (S11) × (cid:90) d r (cid:90) d r γ ( r ) γ ∗ ( r ) × (cid:88) r in H in ( r − r in ) T ( u out , r − r in ) H ∗ in ( r − r in ) T ∗ ( u (cid:48) out , r − r in ) As B R (Eq. S3), B D can be decomposed as the sum of a covariance matrix (cid:104) B D (cid:105) and a per-turbation term δ B D whose intensity decreases with the number M D ∼ (Ω /(cid:96) F ) of independentspeckle grains for the distorted wave-field from the focal plane (Eq. S35). For M D >> , B D converges towards (cid:104) B D (cid:105) , such that: B D ( u out , u (cid:48) out ) = N − in ˆ H ( u out ) ˆ H ∗ ( u (cid:48) out ) (S12) × (cid:90) d r (cid:90) d r γ ( r ) γ ∗ ( r ) × (cid:88) r in (cid:104) H in ( r − r in ) H ∗ in ( r − r in ) (cid:105) T ( u out , r − r in ) T ∗ ( u (cid:48) out , r − r in ) Assuming a strong aberration regime (Eq. S2), the expression of the correlation matrix coeffi-cients B D ( u out , u (cid:48) out ) can be simplified as follows B D ( u out , u (cid:48) out ) = I ˆ H ( u out ) ˆ H ∗ ( u (cid:48) out ) (cid:90) d r | γ ( r ) | (cid:88) r (cid:48) in T ( u out , r (cid:48) in ) T ∗ ( u (cid:48) out , r (cid:48) in ) γ D ( r (cid:48) in ) (S13)with r (cid:48) in = r − r in and γ D ( r (cid:48) in ) = (cid:10) | H in ( r (cid:48) in ) | (cid:11) , (S14)the intensity distribution of the virtual source synthesized in the focal plane at the input. UsingEqs. 3 and S6, Eq. S13 can be rewritten as B D ( u out , u (cid:48) out ) ∝ ˆ H ( u out ) ˆ H ∗ ( u (cid:48) out )ˆ γ D ( u (cid:48) out − u out ) (S15)4here ˆ γ D ( u ) = (cid:80) r γ D ( r ) exp( − j π u . r /λf ) is a discrete 2D Fourier transform of the scatter-ing distribution γ D ( r ) in the focal plane. The correlation length d p of the distorted wave-fieldin the pupil plane is thus inversely proportional to the spatial extension δ in of the input PSFintensity | H in | , such that d P ∼ λf /δ in . (S16)The number of independent speckle grains in the distorted wave-field is the squared ratio be-tween the output pupil size D out = λf /δ and the pupil coherence length d P : N D ∼ ( δ in /δ ) (S17) N D scales as the number of output resolution cells mapping the input PSF.As δ in is smaller than the FOI dimension Ω , d P / N D are larger/smaller than r P / N R (Eqs. S8-S9), respectively. This highlights the enhancement of the far-field correlations in D shownin Fig. 1D3. This figure plots the auto-correlation function B D (∆ u ) of the distorted wave-field in the pupil plane. B D (∆ u ) is computed by averaging the correlation matrix coefficients B D ( u out , u (cid:48) out ) over couples ( u out , u (cid:48) out ) of common relative position ∆ u = u out − u (cid:48) out . S2 Spatial correlations of the reflected and distorted wave-fields
In this section, we derive the input correlations of the matrices R and D . Our aim is to provide atheoretical proof of the experimental observation made in Fig.1C2 and 1D2. As seen previouslyin the pupil plane, the distorted wave-fields reveal spatial correlations in the focal plane thatwere originally hidden in the recorded wave-fields. Unlike the previous section, we derive ageneral expression for the input correlation matrices beyond the isoplanatic limit. The mainresult is the following: While the correlation length r F of the reflected wave-field in the focal5lane is restricted to the input diffraction limit resolution δ , the correlation length d F of D inthe focal plane corresponds to the isoplanatic length (cid:96) c . S2.1 Reflection matrix
To investigate the spatial correlations of the reflected wave-field, the correlation matrix C R = N − out R † R should this time be considered. Unlike in the previous section, the isoplanatic as-sumption is here not made. Using Eq. 2, the coefficients of C R can be expressed as follows: C R ( r in , r (cid:48) in ) = N − out (cid:90) d r (cid:90) d r (cid:48) γ ( r ) γ ∗ ( r (cid:48) ) H in ( r , r in ) H ∗ in ( r (cid:48) , r (cid:48) in ) (cid:88) u out T ( u out , r ) T ∗ ( u out , r (cid:48) ) (S18)As correlation matrices in the pupil plane, C R converges towards the covariance matrix (cid:104) C R (cid:105) for a large number N R ∼ (Ω /δ out ) of independent speckle grains for the reflected wave-fieldin the pupil plane (Eq. S27). For N R >> , the coefficients of C R are given by: C R ( r in , r (cid:48) in ) = N − out (cid:90) d r (cid:90) d r (cid:48) γ ( r ) γ ∗ ( r (cid:48) ) H in ( r , r in ) H ∗ in ( r (cid:48) , r (cid:48) in ) (cid:88) u out (cid:104) T ( u out , r ) T ∗ ( u out , r (cid:48) ) (cid:105) (S19)The mean correlation term (cid:104) T ( u out , r ) T ∗ ( u out , r (cid:48) ) (cid:105) can be developed by writing the transmis-sion matrix as a Hadamard product between the free-space transmission matrix T and an aber-ration matrix H out , such that T ( u out , r ) = ˆ H out ( u out , r ) T ( u out , r ) . It comes (cid:104) T ( u out , r ) T ∗ ( u out , r (cid:48) ) (cid:105) = (cid:68) ˆ H out ( u out , r ) ˆ H ∗ out ( u out , r (cid:48) ) (cid:69) T ( u out , r ) T ∗ ( u out , r (cid:48) ) (S20) = F ( r , r (cid:48) ) (cid:28)(cid:12)(cid:12)(cid:12) ˆ H out ( u out , r ) (cid:12)(cid:12)(cid:12) (cid:29) T ( u out , r ) T ∗ ( u out , r (cid:48) ) . (S21)The correlation function, F ( r , r (cid:48) ) = (cid:68) ˆ H out ( u out , r ) ˆ H ∗ out ( u out , r (cid:48) ) (cid:69) / (cid:28)(cid:12)(cid:12)(cid:12) ˆ H out ( u out , r ) (cid:12)(cid:12)(cid:12) (cid:29) , (S22)6escribes the spatial correlation of the aberration matrix ˆH out in the focal plane. Its supportis directly related to (cid:96) c , the IP size. For sake of simplicity but without lack of generality, weassume that the aberrating layer does not attenuate the wave-field: (cid:28)(cid:12)(cid:12)(cid:12) ˆ H out ( u out , r ) (cid:12)(cid:12)(cid:12) (cid:29) = 1 . (S23)Using Eq. S21, the sum over u out into Eq. S19 can then be rewritten as: N − out (cid:88) u out (cid:104) T ( u out , r ) T ∗ ( u out , r (cid:48) ) (cid:105) = F ( r , r (cid:48) ) (cid:88) u out T ( u out , r ) T ∗ ( u out , r (cid:48) ) (S24)Injecting the expression of the coefficients T ( u out , r in ) (Eq. 3), it finally comes N − out (cid:88) u out (cid:104) T ( u out , r ) T ∗ ( u out , r (cid:48) ) (cid:105) = F ( r , r (cid:48) ) (cid:88) u out exp (cid:18) i πλf u out . ( r − r (cid:48) ) (cid:19) = δ ( r − r (cid:48) ) (S25)The physical meaning of this last equation is that two virtual sources located at points r and r (cid:48) in the focal plane give rise to uncorrelated wave-fields in the pupil plane. Injecting this lastrelation into Eq. S19 leads to the following expression for C R ( r in , r (cid:48) in ) C R ( r in , r (cid:48) in ) = (cid:90) d r | γ ( r ) | H in ( r , r in ) H ∗ in ( r , r (cid:48) in ) (S26)To go further, a rough approximation is to assume an object of constant reflectivity in intensity: (cid:104)| γ ( r ) | (cid:105) = γ . The correlation length r F of the reflected wave-field then corresponds to thecoherence length of the input focal spots. In the strong aberration regime, r F thus scales as theinput diffraction limit δ . The number M R of independent speckle grains in the focal plane thencorrespond to the number of input resolution cells mapping the object: M R ∼ (Ω /δ ) (S27)These theoretical derivations account for the spatial incoherence exhibited by the reflectedwave-field in Fig. 1C2. This figure plots the auto-correlation function C R (∆ r ) of the reflectedwave-field in the focal plane. C R (∆ r ) is computed by averaging the correlation matrix coeffi-cients C R ( r in , r (cid:48) in ) over couples ( r in , r (cid:48) in ) of same relative position ∆ r = r in − r (cid:48) in .7 As highlighted by Fig. 1C and demonstrated above, the reflection matrix displays a randomfeature both at its output and input in the strong aberration regime. Now we will show how thede-scan of the input focal spots in the focal plane reveals the spatial correlations between wavedistortions.In the general case ( i.e beyond the isoplanatic limit), the D -matrix coefficients can be ex-pressed as follows D ( u out , r in ) = (cid:90) d r ˆ H out ( u out , r ) T ( u out , r − r in ) γ ( r ) H in ( r , r in ) (S28)To investigate the spatial correlations of the distorted wave-field in the pupil plane, the correla-tion matrix C D = N − out D † D should be considered. As the other correlation matrices, C D canbe decomposed as the sum of a covariance matrix (cid:104) C D (cid:105) and a perturbation term δ C D whose in-tensity is inversely proportional to the number, N D = ( δ in /δ ) , of independent pupil specklegrains in the distorted wave-field (Eq. S17).For N D >> , C D is shown to converge towards the covariance matrix (cid:104) C D (cid:105) . Its coeffi-cients can then be expressed as follows: C D ( r in , r (cid:48) in ) = N − out (cid:90) d r (cid:90) d r H in ( r , r in ) H ∗ in ( r , r (cid:48) in ) γ ( r ) γ ∗ ( r ) (S29) × (cid:88) u out (cid:104) ˆ H out ( u out , r ) ˆ H ∗ out ( u out , r ) (cid:105) T ( u out , r − r in ) T ∗ ( u out , r − r (cid:48) in ) Using Eqs. 3, 14 and S23, the sum over u out in Eq. S29 can be simplified as follows: N − out (cid:88) u out (cid:104) ˆ H out ( u out , r ) ˆ H ∗ out ( u out , r ) (cid:105) T ( u out , r − r in ) T ∗ ( u out , r − r (cid:48) in )= F ( r , r ) (cid:88) u out exp (cid:18) i πλf u out . ( r − r in − r + r (cid:48) in ) (cid:19) = F ( r , r ) δ ( r − r in − r + r (cid:48) in ) (S30)8f the statistical properties of the scattering medium are invariant by translation, then F ( r , r ) = F ( || r − r || ) . The spatial extension of the function F directly yields the isoplanatic length (cid:96) c .The injection of Eq. S30 into Eq. S29 yields C D ( r in , r (cid:48) in ) = F (∆ r ) (cid:90) d r γ ( r ) γ ∗ ( r − ∆ r ) H in ( r , r in ) H ∗ in ( r − ∆ r , r (cid:48) in ) . (S31)with ∆ r = r in − r (cid:48) in and ∆ r = | r in − r (cid:48) in | . The factor F (∆ r ) requires that the correlationcoefficients C D ( r in , r (cid:48) in ) cancel for points belonging to different IPs. The input PSFs can thusbe considered as locally invariant by translation, such that H in ( r − r in + r (cid:48) in , r (cid:48) in ) (cid:39) H in ( r − r in ) .Equation S31 simplifies into C D ( r in , r (cid:48) in ) ∝ F (∆ r ) (cid:90) d r γ ( r ) γ ∗ ( r − ∆ r ) | H in ( r , r in ) | , (S32)To go further, we can assume that the width of the input focusing beam δ in is larger than thecharacteristic fluctuation length (cid:96) γ of the sample reflectivity: C D ( r in , r (cid:48) in ) ∼ F (∆ r )( γ ∗ γ )(∆ r ) . (S33)where the symbol ∗ stands for the correlation product. Depending on the experimental con-ditions, the coherence length d F of the distorted wave-field can correspond to the correlationlength (cid:96) γ of the object’s reflectivity or the isoplanatic length (cid:96) c associated with the aberratinglayer d F = min { (cid:96) c , (cid:96) γ } (S34) d F is thus always larger than the coherence length r F ∼ δ of the incoherent reflected wave-field (Eq. S15). The number M D of independent focal speckle grains for the distorted wave-fieldis given by M D = (Ω /(cid:96) c ) (S35)If (cid:96) γ > (cid:96) c , this number M D coincides with the number (Ω /(cid:96) c ) of IPs contained by the object.9hese theoretical predictions confirm the experimental observations highlighted by Fig. 1.Spatial correlations are drastically enhanced between the input entries of D (Fig. 1D2) com-pared to R (Fig. 1C2). Figure 1D2 plots the auto-correlation function C D (∆ r ) of the distortedwave-field in the focal plane. This quantity is calculated by averaging the correlation matrix co-efficients C D ( r in , r (cid:48) in ) over couples ( r in , r (cid:48) in ) sharing the same relative position ∆ r = r in − r (cid:48) in .Now, we show how the long-range correlations exhibited by D can be leveraged for over-coming the aberrations and retrieving an image of the object with a resolution close to thediffraction limit. S3 Singular value decomposition of the distortion matrix
To take advantage of the correlations exhibited by the matrix D , its SVD (Eq. 6) is shown tobe an essential tool. It enables a decomposition of the FOI into IMs and an estimation of thetransmission matrix T between the CCD surface and the focal plane. To provide a theoreticalproof of this claim, the previous study of the correlation matrices B D and C D will be helpful.Their eigenvalue decomposition actually dictates the SVD of D . Correlations in the focal planeare shown to predominate in the experiments depicted in the accompanying paper, but also,more generally, in optical microscopy. Strikingly, an exchange of role is noticed between themedium’s reflectivity and the input PSF in the D -matrix compared to the original R -matrix.While the first singular vector of R yields the input PSF for a point-like reflector (
26, 60 ), thefirst singular vector of D directly yields the sample reflectivity for a point-like input focusingbeam in an isoplanatic configuration. Beyond this analogy made between R and D in thisasymptotic limit, a theoretical proof is then provided in the general case. We show how: ( i ) theSVD of D allows a decomposition of the FOI into a set of IMs V p ; ( ii ) a coherent combinationof the output eigenvectors U p can lead to an estimator of the transmission matrix T .10 The SVD of D (Eq. 6) can be directly deduced from the eigenvalue decompositions of itscorrelation matrices B D and C D . The latter ones can actually be written as follows B D = UΣ U † (S36)and C D = VΣ V † . (S37)or, in terms of matrix coefficients, B D ( u out , u (cid:48) out ) = N in (cid:88) p =1 σ p U p ( u out ) U ∗ p ( u (cid:48) out ) . (S38)and C D ( r in , r (cid:48) in ) = N in (cid:88) p =1 σ p V p ( r in ) V ∗ p ( r (cid:48) in ) . (S39)The eigenvalues of B D and C D are the square of the singular values σ p ; their eigenvectors, U p and V p , are the output and input singular vectors, respectively. The SVD of D is dictated eitherby the correlations between its lines or columns. To know which ones dominate over the other,the analytical expressions of the correlation matrices, B D and C D , should be investigated (seeEqs. S15 and S33).If the reflectivity of the object was fully random, i.e (cid:104) γ ( r ) ∗ γ ( r ) (cid:105) = δ ( r ) , the correlationmatrix C D (Eq. S33) would be diagonal. This means that the columns of D would be fullyuncorrelated. On the contrary, output correlations would subsist in D as they only depend onthe spatial extension of the input focal spot (Eq. S16). In this random speckle regime, the SVDof D is dominated by its correlations in the pupil plane and the analysis of D should rather berestricted to a FOI containing a single IP. This regime has been recently investigated in medicalultrasound imaging ( ) where scattering is often due to a random distribution of unresolvedscatterers. 11n optical microscopy, biological tissues induce a strong forward scattering: The involvedscatterers display a characteristic length (cid:96) γ larger than the wavelength. The auto-correlationof the sample reflectivity can span over several IPs especially at large depths. In this forwardscattering regime, correlations of the distorted wave-field in the focal plane may dominate overits far-field correlations.To know if this is the case, one can compare the number of independent speckle grains, N D and M D , in the pupil and focal planes, respectively. The correlation degree in each plane isactually inversely proportional to this number. Correlations in the focal plane will thus dominateif N D > M D . The latter condition is fulfilled in a strong aberration regime for which the numberof output resolution cells mapping each aberrated focal spot, N D = ( δ in /δ ) , is larger thanthe number of IPs mapping the object surface, M D = (Ω /(cid:96) c ) . This condition is checked inthe experiments of the accompanying paper. For instance, in the experiment depicted in Fig. 4, N D ∼ while M D ∼ . S3.2 Analogy with iterative time reversal
Now that the conditions for a domination of correlations in the focal plane have been derived,we now study the singular vectors of D . To that aim, an analogy with iterative time reversal isfirst explored to give a physical intuition of the SVD of D .If we compare the analytical expressions of the correlation matrices C R (Eq. S13) and C D (Eq. S31), we can notice an exchange of the role between the medium reflectivity γ and theinput PSF H in . While C R corresponds to a static object scanned by a moving illuminating beam(Fig.2A), C D corresponds to a static focused beam illuminating a moving object (Fig.2B). Inthe isoplanatic limit, the distortion matrix D (Eq. S31) is thus equivalent to a virtual reflectionmatrix associated with: ( i ) a coherent reflector of scattering distribution | H in ( r ) | (located on theoptical axis and at the focal plane); ( ii ) a virtual focusing beam associated with the PSF γ ( r in + ) . As shown by iterative time reversal experiments ( ), the reflection matrix is of rank 1 fora point-like scatterer, and its first input singular vector V shall directly yield the virtual inputPSF ( ). By analogy, for a point-like input focusing beam, the D -matrix shall be also of rank1 and its first input singular vector V shall directly provide the medium reflectivity γ ( r in ) .Interestingly, the SVD of D should therefore unscramble aberrations and sample reflectivity.However, this qualitative analysis has been made under strong hypotheses: the isoplanatic limitand a point-like input focusing beam. In the following, we make the problem more complexby first going beyond the isoplanatic limit and then by considering the finite size of the inputfocusing beams. S3.3 Isoplanatic modes
Let us first assume a point-like input focusing beam, H in ( r , r in ) = | H in ( r in , r in ) | δ ( r − r in ) ,beyond the isoplanatic limit. Equation S31 becomes C D ( r in , r (cid:48) in ) ∝ F (∆ r ) γ ( r in ) γ ∗ ( r (cid:48) in ) H in ( r in , r in ) H ∗ in ( r (cid:48) in , r (cid:48) in ) . (S40)A full-field intensity image F ( r in ) of the sample reflectivity can be retrieved by considering thediagonal of C D : F ( r in ) = C D ( r in , r in ) = | γ ( r in ) | | H in ( r in , r in ) | (S41) F ( r in ) can be a satisfying estimator of the sample reflectivity, | γ ( r in ) | . However, the inputfocusing beam intensity H in ( r in , r in ) pollutes the full-field image. The latter term can be detri-mental to imaging since it gives rise to a fluctuating contrast across the focal plane. Moreover,experimental noise and diffusive multiple scattering can still degrade the image. At last, we maywant to have access to the amplitude and phase of the reflectivity rather than only its squarenorm. For all these reasons, the singular value decomposition of D (Eq.6), or equivalently,the eigenvalue decomposition of C D (Eq.S39) is decisive. In the general case, the correlation13unction F (∆ r ) (Eq. S22) governs the eigenvalue decomposition of C D . The ratio betweenthe object surface Ω and the isoplanatic area (cid:96) c yields the effective rank M D = (Ω /(cid:96) c ) of C D . This rank scales as the number of IPs that fit in the object. The input eigenvectors V p can be derived by solving a second order Fredholm equation with Hermitian kernel ( ). Ananalytical solution can be found for certain analytical form of the correlation function F (∆ r ) (Eq. S22). For instance, a sinc kernel imply 3D prolate sphero¨ıdal eigenfunctions ( ); a Gaus-sian covariance function leads to Hermite-Gaussian eigenmodes ( ); exponential or triangularkernels yields cosine and sine eigenfunctions ( ). A general trend is that the spatial frequencycontent of the eigenvectors increases with their rank.The identification of Eqs. S39 and S40 leads to the following equality: M D (cid:88) p =1 σ p V p ( r in ) = H in ( r in , r in ) γ ( r in ) (S42)A coherent combination of the M D first eigenvectors V p can yield the amplitude and phaseof the reflectivity but the result is still polluted by the input illumination beam H in ( r in , r in ) . Inpractice, aberrations at the input can be corrected through the same process by exchanging inputand output, i.e by projecting the data in the pupil plane at the input and in the focal plane at theouput. In the experiments depicted in the accompanying paper, the sparse illumination schememakes the input basis incomplete and the spatial sampling insufficient. The image should thusbe built from the output to benefit from the excellent resolution with which the field is recordedby the CCD camera. To do so, Eq. S42 can be used to prove that the coherent combination ofoutput singular vectors U c = (cid:80) M D p =1 U p (Eq. 6) perfectly compensate for the output aberrationmatrix ˆH out . To that aim, let us apply the transpose conjugate U † c to the output of the matrix D (Eq.S28). It comes: (cid:90) d r (cid:88) u U ∗ c ( u ) H out ( u , r ) T ( u , r − r in ) γ ( r ) H in ( r , r in ) = H in ( r in , r in ) γ ( r in ) (cid:88) u U ∗ c ( u ) H out ( u , r ) T ( u , r − r in ) = δ ( r − r in ) (S43)which, under the matrix formalism, can be rewritten as ( U c ◦ T ) † T = I (S44)The matrix ˆT = ( U c ◦ T ) is an estimator of the transmission matrix T . The applicationof its transpose conjugate, ˆT † enables a perfect compensation for the aberrations contained inthe transmission matrix T . To obtain a diffraction-limited image of the object, the matrix ˆT † should be directly applied to the output of the matrix R (Eq.12 of the accompanying paper).This operation leads to a focused matrix R F whose coefficients can be expressed as R F ( r out , r in ) = γ ( r out ) H in ( r out , r in ) (S45)This matrix consists in an Hadamard product between the reflectivity of the focal plane at itsoutput and the input focusing matrix. In other words, aberrations are corrected at the output butsubsists at the input. Hence the resulting confocal image built from the diagonal of R F suffersfrom the same issue: I ( r out ) = R F ( r out , r out ) = γ ( r out ) H in ( r out , r out ) (S46)It is a relying estimator of the object’s reflectivity γ ( r out ) , but modulated by the amplitude andphase of the input illumination H in ( r out , r out ) . To reduce this detrimental effect on the imagecontrast, one can consider a full-field image integrated over all input focusing beams (see Eq. 14of the accompanying paper) or an adaptive confocal image integrating over a numerical pinhole(see Eq.18 of the accompanying paper). This integration over r in allows us to smooth themodulation of the image induced by the input focusing beams.15 All these theoretical developments have been made by considering a point-like input focusedbeam. This is, of course, not true in reality. The input focusing beam gives rise to a virtualcoherent reflector of finite size δ in . The issue we want to address is the impact of this size onthe SVD of D . Assuming incoherent input focusing beams (Eq. S2), Eq. S32 can be rewrittenas follows in the isoplanatic limit (Eqs. S5-S6): C D ( r in , r (cid:48) in ) ∝ (cid:18)(cid:90) d r γ ( r ) H in ( r − r in ) (cid:19) × (cid:18)(cid:90) d r (cid:48) γ ( r (cid:48) ) H in ( r (cid:48) − r (cid:48) in ) (cid:19) ∗ , (S47)By confronting this last equation with the eigenvalue decomposition of Eq. S39, it turns out thatthe distortion matrix D is of rank 1 and that its input singular vector V can be expressed as V ( r in ) = (cid:90) d r γ ( r ) H in ( r − r in ) = [ γ (cid:126) H in ]( r in ) . (S48)where the symbol (cid:126) stands for the convolution product. Albeit independent from output aber-rations, V is nevertheless a convolution product between the object’s reflectivity and the inputPSF H in (see Fig 2C of the accompanying paper). The output singular vector U can be deducedfrom V through the following matrix product: σ U = DV . (S49)Injecting Eq. S10 and Eq. S48 into this last equation yields the following expression for thecoefficients of U σ U ( u out ) = ˆ H out ( u out ) (cid:90) d r (cid:90) d r (cid:48) (cid:88) r in T ( u out , r − r in ) γ ( r ) γ ( r (cid:48) ) H in ( r − r in ) H in ( r (cid:48) − r in ) For a large number of resolution cells in the FOI, U will converge towards its ensemble aver-age, such that σ U ( u out ) = ˆ H out ( u out ) (cid:90) d r (cid:90) d r (cid:48) (cid:88) r in T ( u out , r − r in ) γ ( r ) γ ( r (cid:48) ) (cid:104) H in ( r − r in ) H in ( r (cid:48) − r in ) (cid:105)
16n a strong aberration regime (Eq. S2), the last equation can be rewritten as follows σ U ( u out ) = ˆ H out ( u out ) (cid:90) d r | γ ( r ) | (cid:88) r in T ( u out , r − r in ) (cid:10) | H in ( r − r in ) | (cid:11) Using the expression of the free-space transmission coefficients T ( u out , r in ) (Eq. 3), it finallyturns out that U ( u out ) ∝ ˆ H out ( u out ) (cid:104) ˆ H in ∗ ˆ H in (cid:105) ( u out ) . (S50)and σ ∝ (cid:90) d r | γ ( r ) | . (S51)While the singular value σ yields the object’s reflectivity integrated over the associated iso-planatic patch (here the FOI), the vector U corresponds to the aberration output transmittance ˆ H out modulated by the autocorrelation function of the aberration input transmittance ˆ H in (seeFig 2C of the accompanying paper). This last term tends to limit the angular aperture of thesingular vector U by the coherence angle of the input aberration ˆ H in . To circumvent that issue,the trick is to consider only the phase of the first singular vector U . Indeed, if we make therealistic hypothesis of a real and positive autocorrelation function ˆ H in ∗ ˆ H in , the normalizedvector ˜U is then given by ˜ U ( u out ) = exp ( j arg { U ( u out ) } ) = ˆ H ( u out ) (S52)A novel input vector ˜V can then be retrieved through the matrix product: ˜U † D = ˆV . (S53)Injecting the expression of ˜U (Eq. S52) and D (Eq. S10), the following expression can beretrieved for ˆV in the isoplanatic limit: ˆ V ( r in ) = H in ( r in , r in ) γ ( r in ) (S54)17f we compare this last equation with Eq. S48, the normalization of U allows us to virtuallyreduce the size of the input focusing beam (see Fig 2D of the accompanying paper). The matrix ˆT = ( ˜U ◦ T ) is then a satisfying estimator of the transmission matrix T in the isoplanaticlimit. The application of its transpose conjugate, ˆT † , allows a perfect compensation for theaberrations contained in the transmission matrix T . A diffraction-limited image of the objectcan be obtained by applying the matrix ˆT † to the output of matrix R (Eq. 12). S3.5 General case
In the general case ( i.e beyond the isoplanatic limit), the same method can be employed tovirtually reduce the size of the focal spot over each IM. The corresponding singular vectorsshould be normalized: ˜U p = exp ( j arg { U p } ) . The application of their transpose conjugateto the R -matrix should lead to an optimal aberration correction over each IM at the output.One open question is whether these output singular vectors can be combined coherently ornot, such that ˜U c = (cid:80) p ˜U p . In the present work, this coherent combination does not providebetter results than an incoherent summation of each IM image I p (Eq.20). This is because anincoherent sum of R F -matrix coefficients at the input is required to smooth the modulation ofthe image by H in (Eqs.14-18). 18 upplementary Figures Laser CCDSLM MirrorPZTP P P L L L L k o u t AB CD
Figure S1:
Measuring the time-gated reflection matrix . Experimental set up: P: polarizer,MO: microscope objective, BS: beam splitter, PBS: polarized beam splitter, SLM : spatial lightmodulator, PZT: piezo phase shifter, M: Mirror. A femtosecond laser beam (center wavelength:810 nm, bandwidth: 40 nm) is shaped by an SLM acting as a diffraction grating. A set of inci-dent plane waves is thus emitted from the SLM and focused at a different position in the focalplane of an immersion MO (NA=0.8). The backscattered wave-field is collected through thesame MO and interferes with a reference beam on a CCD camera. The latter one is conjugatedwith the back focal plane of the MO. The amplitude and phase of the wave-field is recordedby phase shifting interferometry ( ). The time of flight t is controlled by the length of theinterferometric arm and is matched with the position of the focal plane. For each input focusingpoint r in , a reflected wave-field is recorded in the pupil plane and stored along a column vectorin the matrix R ) = [ R ( u out , r in )] . 19 nput pupil output pupilobjective pupil D = 2 NA x f D in = 2 NA in x finput pupil aperture D out = 2 NA out x foutput pupil apertureobjective pupil D = 2 NA x f F O I diffraction limited input focal spot δ in =λ/2NA in output resolution cell δ out =λ/2NA out0 scanning stepδr in0 aberrated output focal spot δ out aberrated input focal spot δ in F O V Figure S2:
Conjugation relationships between the pupil, focal and imaging planes . Thedistortion matrix connects the focal plane of the MO with the output pupil plane. This figureillustrates the various parameters involved in the different planes of the system. Both inputand output pupil planes are ultimately limited by the MO edges. The input pupil D in is evenmore limited because the illumination beam underfills the objective pupil D . It results in areduced NA denoted NA in . In turn, the size of the input focal spot in the image plane is givenby δ = λ/ NA in if there is no aberration and δ in in the general case. In this focal plane, thefield-of-illumination depends on the scanning step (spatial sampling), denoted as δr in , and thenumber of measurements N in . In reflection, the output pupil D out is also smaller than the totalobjective pupil D due to the limited surface of the detector but larger than the input pupil D in .The resolution of the image is thus governed by the output resolution cell δ = λ/ NA in .20
00 400 600 800 1000depth (µm) s i ng l e - t o - m u l t i p l e sc a tt e r i ng r a t i o scattering mean free path l s (µm) s i ng l e t o m u l t i p l e sc a tt e r i ng r a t i o -1 -2 A
50 100 150 20010 -4 -2 edematous opaque cornea rat intestine tissues B o p t i c a l c o h e r e n c e t o m o g r a p h y m a t r i x a pp r o a c h imagingthreshold SMR=1 Figure S3:
Predicting the single-to-multiple scattering ratio in biological tissues. ( A ) SMRas a function of depth for the imaging experiment through the rat intestinal tissue (see Figs. 1and 3). The red curve (before aberration correction) is plotted for a Strehl ratio S = 3 × − ,while the blue curve (after matrix correction) corresponds to S = 1 . × − . The detectionthreshold yields an imaging depth limit of ∼ µ m for conventional OCT and 900 µ m forour matrix approach. ( B ) SMR as a function of the scattering mean free path (cid:96) s for the corneaimaging experiment (see Fig. 5). A SMR of 1 is obtained for a scattering mean free path (cid:96) s ∼ µ m. In both panels, the theoretical curves are built by considering the model described inRef. ( ) and the experimental parameters described in the paper. Note also that the y-axis is inlog-scale. 21 ( u out , r in ) u o u t r in x in y i n v out w o u t N e l e m e n t s N e l e m e n t s Figure S4:
Building the reflection matrix R . For each focused illumination at r in = ( x in , y in ) ,the reflected wave-field ψ r in ( v out , w out ) is recorded in the pupil plane by each pixel of the CCDcamera whose position is denoted by the vector u out = ( v out , w out ) . Each wave-field is concate-nated and stored along a column vector. This set of column vectors forms the reflection matrix R = [ R ( u out , r in )] , such that R ( u out , r in ) = ψ r in ( v out , w out ) .22 out u in CCD SLM T( u out , r ) H in ( r,r in ) r in γ( r) R( u out , r in ) ABC D
Figure S5:
Modeling light propagation from the virtual source plane to the output pupilplane.
The reflection matrix R contains the impulse responses R ( u out , r in ) between each virtualsource point r in and each CCD pixel u out in the output pupil plane. ( A ) The virtual source point r in is produced by each transverse mode shaped by the SLM in the input pupil plane. ( B ) Thepropagation between the virtual source plane and the focal plane of the MO can be modelledby the input focusing matrix H in = [ H in ( r , r in )] whose columns correspond to each input focalspot in the sample plane for each incident focusing point r in . ( C ) The return trip of the wavefrom the sample to the CCD camera is modeled by the transmission matrix T = [ T ( u out , r )] that connects each point r in the focal plane to each pixel u out of the CCD camera. ( D ) Finally,based on these propagation matrices and the sample reflectivity matrix Γ , the reflection matrix R can be simply expressed as the matrix product of these three matrices (Eq.1).23 upplementary tables variable definition λ wavelength n optical index g anisotropy factor (cid:96) s scattering mean free path L thickness of the scattering layer d distance between the aberrating layer and the focal plane f focal length of the microscope objective Ω size of the field-of-illumination r in / r out position vector in the input focusing / output pupil planes r position vector in the focal plane of the microscope objective u out position vector in the output pupil plane N in / N out number of input focusing beams / pixels in the output pupil plane D in / D out input / output pupil aperture NA in / NA out input / output numerical aperture δr in / δu out spatial sampling in the input focusing / output pupil planes δ in / δ out characteristic width of the input/output point spread functions δ / δ diffraction limit resolution at input/ouput r P / r F correlation length of the reflected wave-field in the output pupil / input focusing plane d P / d F correlation length of the distorted wave-field in the output pupil / input focusing plane (cid:96) γ charateristic correlation length of the sample’s reflectivity (cid:96) c characteristic size of an isoplanatic patch M D / N D number of independent speckle grains for D in the input focusing / output pupil planes M R / N R number of independent speckle grains for R in the input focusing / output pupil planes S / S p / S (cid:48) p Strehl ratios: initial /final / weighted values σ p / ˜ σ p singular values of D : raw / normalized H (˜ σ p ) Shannon entropy of singular valuesSMR single-to-multiple scattering ratioTable S1: