Multiplexed wavefront sensing with a thin diffuser

In astronomy or biological imaging, refractive index inhomogeneities of e.g. atmosphere or tissues induce optical aberrations which degrade the desired information hidden behind the medium. A standard approach consists in measuring these aberrations with a wavefront sensor (e.g Shack-Hartmann) located in the pupil plane, and compensating them either digitally or by adaptive optics with a wavefront shaper. However, in its usual implementation this strategy can only extract aberrations within a single isoplanatic patch, i.e. a region where the aberrations remain correlated. This limitation severely reduces the effective field-of-view in which the correction can be performed. Here, we propose a new wavefront sensing method capable of measuring, in a single shot, various pupil aberrations corresponding to multiple isoplanatic patches. The method, based on a thin diffuser (i.e a random phase mask), exploits the dissimilarity between different speckle regions to multiplex several wavefronts incoming from various incidence angles. We present proof-of-concept experiments carried out in wide-field fluorescence microscopy. A digital deconvolution procedure in each isoplanatic patch yields accurate aberration correction within an extended field-of-view. This approach is of interest for adaptive optics applications as well as diffractive optical tomography.


INTRODUCTION
Image degradation caused by aberrations can be a critical issue when imaging objects of interest located behind a complex medium, e.g. a turbulent atmosphere or biological tissues, in which the distribution of the refractive index is both heterogeneous and dynamic.To correct aberrations and increase the image contrast and resolution, the most common technique is to use pupil adaptive optics (AO), in which the aberration is first characterized in the pupil plane of the optical system, and then corrected either digitally or by a wavefront shaper [1][2][3][4].
Aberration characterization methods can be classified into either direct or indirect wavefront sensing techniques.Indirect techniques estimate aberrations by optimizing, over time, a feedback image or signal with certain metrics (e.g.image contrast, sharpness or intensity [5][6][7], nonlinear signals [8][9][10] etc.).These techniques, also referred to as "sensorless", do not require any dedicated optical device to estimate aberrations, thus potentially reducing cost and hardware complexity.However, the convergence of the optimization process can be unstable and/or may lead to local minima.More importantly, these techniques require multiple measurements, which limits their application to static, or slowly varying, aberrating media, thus excluding most astronomical and ophthalmologic applications.
The other class of approaches, direct wavefront sensing, relies on a deterministic measurement of aberrations, typically using point sources called guide-stars (GS).The most common wavefront sensor (WFS) is indisputably the Shack Hartman [11].This WFS generally samples the pupil plane with a micro-lens array, where each micro-lens directly yields the local phase gradient of the wavefront in the corresponding subset of camera pixels, which defines a macropixel.A significant advantage of this method is that it only requires a single image acquisition to measure a wavefront, thus allowing higher temporal resolution.Direct wavefront sensing is preferred -and often crucial -whenever dynamic aberrations are involved, like those caused by the atmospheric turbulence in astronomy or the lachrymal film in ophthalmology [4,12].
However, aberrations only remain correlated within a limited region called the isoplanatic patch, whose typical size depends on the position and properties of the aberrating medium.A pupil aberration measurement typically assumes a 2D projection of the 3D refractive index of this medium along the direction of the GS.Wavefront measurements obtained with a single GS can therefore only be used to correct the aberration within a limited field-of-view (FoV) in the case of volumetric aberrating media.Aberration correction is then only valid in the angular vicinity of the GS, but decorrelates away from it.Extensive efforts have been made to address this critical issue using direct or indirect wavefront sensing methods.Conjugate AO can increase the isoplanatic patch size by conjugating the WFS (and the corrector) to the main aberrating layer [13][14][15][16].This approach is therefore highly efficient to widen the FoV in situations where aberrations originate from a single, well-defined planar layer.However, it requires a much larger number of sensing/shaping modes and it is not optimal in the most common scenarios where aberrations originate from a non-planar interface, from several layers, or from a volume.Ideally, a tomographic characterization of the aberrating medium is needed to retrieve high resolution over a large FoV or volume.This concept is at the heart of Multi-Conjugate AO, which uses several WFS (typ.>6) to probe aberrations along several GSs directions in order to estimate volumetric aberrations and correct them in several planes [17,18] .While increasingly used in astronomy, this approach is too complex, cumbersome and expensive to be routinely applied in microscopy and ophthalmology.Some approaches operating in the pupil plane rather than in a conjugated plane have also emerged.An interesting method sequentially applies various masks in the pupil plane to retrieve the local phase gradient in multiple isoplanatic patches [19][20][21].Alternatively, the beam or GS can be scanned to sequentially measure pupil aberrations in different locations [22].Both approaches however require a time-consuming scan which degrades temporal resolution.Shack-Hartmann WFS could in principle provide parallelized multi-GS, multi-angle measurements in a single shot but their angular dynamics is limited by aliasing problems and their ability to measure multiple, distributed GSs is strongly constrained (See Supplementary S1).
While some of these emerging methods increasingly allow aberration characterization along multiple angles, they can hardly address dynamic modifications.Conversely, fast methods are essentially limited to a single GS.In microscopy, the speed constraint is partially relieved because aberrations in tissues typically evolve over time scales of a few minutes.However, the development of fast, multiplexed wavefront sensing methods remains essential when imaging moving organisms, or deep inside tissues, where wavefronts are also affected by rapidly evolving, multiple scattering [23].It is also crucial for the volumetric imaging (e.g.light sheet, light field, multiplane imaging …) of large samples since aberrations need to be characterized for a high number of isoplanatic volumes, which is time-consuming with indirect or sequential methods [24][25][26][27].In most fields requiring aberration correction, a simultaneous parallel measurement of multiple pupil aberrations would provide a broad FoV.Such approaches could also put single shot tomographic microscopy within reach, while drastically simplifying the instrumentation.
In this context, recent advances in light manipulation through scattering media appear as decisive.The unique signature associated to random speckle patterns makes them particularly suited to the unambiguous retrieval of multiplexed information.Notably, the use of diffusers as encoding masks has recently shown its potential in areas including compressive ghost imaging [28], single-exposure 2D or 3D imaging [29][30][31], hyperspectral imaging [32][33][34] and 3D super-resolution microscopy [35].Furthermore, the analysis of speckle pattern distortions has recently emerged as an accurate way to perform high-resolution wavefront sensing [36][37][38][39].
This paper addresses the single-shot multiplexed measurement of multiple pupil aberrations, i.e. the parallel measurement of several wavefronts coming from different propagation directions.To solve the assignment ambiguity problem inherent to multiple wavefronts and periodic mask patterns (e.g.Shack-Hartmann), we propose to use a thin diffuser (i.e. a random phase mask) as Hartmann mask.The wavefronts related to several GSs are thus encoded into an incoherent sum of distorted speckle patterns.By exploiting both the dissimilarity between different speckle regions and the large memory effect of a thin diffuser [40,41], we show that each speckle pattern can be unambiguously ascribed to its GS, allowing to extract the associated phase gradient maps from speckle pattern distortions.We experimentally validate this concept in microscopy with the multiplexed acquisition of wavefronts originating from several fluorescent GSs.This demonstrates the efficacy of this single shot method for the characterization of pupil aberrations in multiple isoplanatic patches.Finally, we show in a proof-of-concept experiment that this multiplexed WFS can feed a deconvolution procedure, to digitally correct aberrations in several isoplanatic patches and obtain highresolution images in a large FoV.

PRINCIPLE
Figure 1(a) illustrates the general concept of multiplexed wavefront sensing in the context of fluorescence microscopy.Fluorescent GSs are imaged by a microscope objective, through an aberrating medium.The latter induces spatially varying pupil aberrations as the GSs are located in different isoplanatic patches.In the pupil plane, each GS will then generate an aberrating wavefront carrying a global tilt offset according to the GS lateral position.A thin-diffuser-based WFS, conjugated to the pupil plane, then allows measuring these aberrating wavefronts incoming with different propagation directions from a single image acquisition.
Before detailing the multiplexed wavefront sensing concept, let's first recall the basics of diffuser-based WFS in the case of a single wavefront acquisition.We recently demonstrated that a thin diffuser can be used as a Hartmann mask to perform broadband wavefront sensing with a high accuracy and resolution [39,42,43].The WFS is composed of a thin diffuser used as a Hartmann-type mask located at a short distance d from a camera sensor [see Fig. 1(b)].A "reference" speckle pattern () is first acquired by illuminating the WFS with a plane wave.Due to the large angular memory effect of thin diffusers [40,41], a tilt to the wavefront does not modify or decorrelate the speckle pattern but simply translates it.For a distorted wavefront, speckle grains will be locally shifted on the detector according to the local wavefront gradient, as shown in Fig. 1(b): a speckle grain at position r is locally displaced by a vector (), and the resulting speckle map () is distorted as compared to the reference () according to [36,43]: where () is the normalized intensity map.The phase gradient map  ⊥  can then be estimated by measuring the overall speckle grain displacement map u [39]: with  0 the wavenumber.
Equation 2 shows that thin-diffuser-based WFS and microlensarray-based Shack Hartmann WFS have similar behaviors, with speckle grains playing a role similar to that of focal spots in each subset of camera pixels, or macropixel (typ.7x7 camera pixels [39]), and the refractive diffuser surface behaving as random, aberrating microlenses [29].Interestingly, the use of a thin diffuser fundamentally provides a major advantage in the context of wavefront multiplexing.While the intensity pattern generated by micro-lens arrays in a Shack Hartmann is periodic, and thus prone to ambiguity between identical, unrecognizable spots, a diffuser generates a superposition of unique patterns which can be identified unambiguously.The ambiguity problem strongly limits the phase dynamic range and the number of wavefronts that can be multiplexed with a Shack Hartmann (See Supplementary S1).In contrast, random speckle grains provide a unique signature.Two uncorrelated, random speckle patterns are indeed statistically orthogonal relative to a zero-mean cross-correlation product.This inherent property of speckles alleviates the ambiguity issue and allows to retrieve multiple wavefronts locally encoded in the form of multiple, orthogonal speckle patterns [35].
To illustrate this concept, Fig. 1(c) and (d) respectively show a reference speckle pattern R and a multiplexed speckle pattern M encoding three wavefronts.Note that although the three speckle patterns associated to each wavefront are represented using three colors for clarity, the concept clearly applies to a single wavelength, monochrome detector, and to more than N=3 wavefronts.Since the N wavefronts originate from different incoherent point sources, the pattern M can be described as the sum of N reference speckle patterns R, each shifted and deformed by N different non-rigid transformations   (r): where j stands for the index of the wavefront among the N multiplexed ones (N=3 in Fig. 1).This superposition is illustrated for a chosen subset of (), the macropixel   () (here e.g.12x12 pixels) shown in  The reference speckle pattern R and the multiplexed speckle pattern M measured by the gray-level camera contain information on the three wavefronts.M is the incoherent sum of three distorted speckle patterns related to each wavefront.The pattern measured within a chosen region, the macropixel MA (gray square in d.), can be described as the sum of three patterns of the reference (green, blue and red boxes in c.) each shifted by uj, with j=1,2,3, the index of each wavefront.(e) Owing to the orthogonality between different speckle regions, the cross-correlation of   () with the full reference speckle pattern R reveals N=3 peaks.The peak position gives access to the local speckle grains displacements related to each wavefront, and so, to (f) the phase gradient maps of the wavefronts.(g) A 2D integration step finally allows independent wavefronts reconstructions.
() is the incoherent sum of three different subregions of the reference speckle R having undergone different lateral shifts,  1 ,  2 and  3 (See zoomed macropixel in Fig. 1(c) and (d)).Since these subregions correspond to different regions of the speckle pattern created by different parts of the diffuser Hartmann mask, they are statistically orthogonal.The cross-correlation   () ⋆  therefore allows to extract, without ambiguity, the speckle grain displacements vector map   () associated to each wavefront through an estimate (centroid) of the position of the peaks, and the normalized intensities   through the maximum values in each peak.This process is illustrated in Fig. 1(e), where three correlation peaks yielding  1 ,  2 and  3 are clearly visible.The local phase gradient at the pupil coordinate r can then be estimated using Eq. ( 2).Using a Digital Image Correlation (DIC) algorithm [44], the same process can be reproduced in all macropixels   (′) of the full multiplexed speckle pattern M in order to extract the phase gradient vector maps associated to each wavefront [see Fig. 1(f)].Finally, a 2D integration of these phase gradient maps can be used to retrieve the N wavefronts [see Fig. 1(g)].
It should be noted that this method implicitly requires a sparsity assumption about GSs density which is clearly visible in Fig 1(e): the position of the correlation peaks, related to the maximum phase gradients of wavefronts (delimited here by dashed circles), can only be estimated if the correlation peaks do not overlap (See Supplementary S2).However, this constraint is common to all WFS, and is not critical to the targeted application since wavefronts coming from GS located in different isoplanatic patches are generally well separated angularly.Furthermore, as discussed in the experimental section, the validity of this assumption can be predicted by calculating the "global" cross-correlation  ⋆ .
This approach can in principle be applied to any number of GS and wavefronts.However, the contrast of the multiplexed speckle pattern () resulting from the superposition of N patterns ( ∝ 1 √ ⁄ ) decreases for large N-value.This reduces the signal to noise ratio of the cross-correlation map   () ⋆ , the accuracy of the determination of speckle grain displacement   , and thus the accuracy of wavefront reconstruction.A first solution to maintain a high reconstruction fidelity while multiplexing a large number N of wavefronts would be to use larger macropixels MA (e.g.45x45 pixels, Supplementary S3).However, this method drastically degrades the spatial resolution on the rebuilt wavefront.To alleviate resolution degradation, we propose an iterative DIC process to converge towards each speckle pattern   [See Eq. ( 4)]: after a first DIC step (which provides a first estimation of the intensities   and displacement maps   ), the distorted speckle pattern   can be isolated from the multiplexed speckle pattern  by subtracting the contributions of all others GSs  ≠  : This process restores the contrast of speckle pattern   as compared to  and thus allows, through a second DIC step, to recover a more accurate estimation of the intensities   and displacement maps   .Noteworthy, since the subtracted patterns are prone to uncertainties, the contrast of the pattern   is restored but remains noisy before the algorithm converges.For this reason, this process can be repeated T times to increase the reconstruction accuracy and can be seen as a gradient descent algorithm which iteratively minimizes the quantity A detailed description of the corresponding algorithm, as well as numerical simulation results demonstrating its relevance for high-resolution wavefront sensing are provided in Supplementary S3.Briefly, these simulations show that the number T of iterations which are necessary to retrieve N wavefronts increases with N (typ.T=2 for N=5 wavefronts).Supplementary S3 also shows that more than N=16 wavefronts can be retrieved using this iterative DIC approach, while reducing the RMS error by an order of magnitude compared to direct DIC (See Fig. S6.).Importantly, it also shows that high-resolution is preserved by this method since small phase pixel size (typ.7x7 pixels) can be reached.When using a 4.2M pixel camera, this provides typ.16 multiplexed wavefronts, each with 85K phase and intensity pixels.

Description of the setup
To validate this concept experimentally, we built a wide-field fluorescence microscope based on a commercial inverted microscope (Olympus IX-71).1µm-diameter fluorescent beads randomly distributed on a glass slide are used as sample (Orange 540/560, Thermo Fisher).An aberrating layer (1° Holographic Diffusers, Edmund Optics) is positioned 150 µm away from the microbeads sample to induce spatially varying pupil aberrations [see Fig. 1(a)].In order to excite the fluorescence of multiple beads chosen within the FoV, a 532nm laser beam is shaped using a phaseonly spatial light modulator (SLM, X13138-01, Hamamatsu) conjugated to the back focal plane of the microscope objective (LUCPlanFLN, NA=0.45, x20, Olympus).The computer-generated hologram displayed on this SLM is calculated using a Gerchberg-Saxton algorithm [45] so as to illuminate the chosen beads and use them as GSs.The 532 nm excitation is spectrally filtered by a dichroic mirror (NFD01-532, Semrock) and a notch filter (λ=533 ± 2nm, Thorlabs) to collect the emitted fluorescence (See Supplementary S4).
The WFS is composed of a thin diffuser (1° scattering angle holographic diffuser, Edmund Optics) and a sCMOS camera (Zyla 5.5, Andor).The diffuser-camera distance is set to d=3 mm (see Ref. [39]), but placing the diffuser so close to the sensor was not mechanically possible.A 1x magnification relay lens (not shown in Fig. 1(a)) is therefore used to image the diffuser at a distance d from the camera.To properly measure wavefront distortions coming from various GSs, the multiplexed WFS is conjugated with the back pupil plane of the microscope objective.The reference speckle pattern is acquired in a first step using a simple collimated beam generated from a multimode fiber (core diameter 10µm, Thorlabs) and a long focal length lens (f'=400mm).This reference speckle pattern is shown in Fig. 2(a).

Multiplexed wavefront sensing validation
To demonstrate the possibilities of the method, we first excited simultaneously N=3 GSs located in different isoplanatic patches, i.e. separated by more than 120 µm in the sample plane (See Supplementary S4 for the characterization of the sample isoplanatic patch size).On the WFS, their fluorescence yields a superimposition of 3 speckle patterns, as shown in the bottom of Fig. 2(a).The crosscorrelation  ⋆  between the reference and multiplexed speckles is shown in Fig. 2(b), clearly revealing the number and location of the three excited GSs.Here, the position of each correlation peak indeed indicates the average propagation direction (or global tip/tilt: ⟨⟩  = ⟨⟩   ⁄ ) corresponding to each GS.The absence of overlap between peaks ensures that the sparsity hypothesis is valid, and that the wavefronts can be reconstructed independently.To this aim, the iterative DIC algorithm (T=3 iterations) is used to recover the speckle grains displacement maps related to each wavefront.The wavefronts obtained after integration of the phase gradient maps are shown in Fig. 2(c).
The validity of these multiplexed measurements was compared to individual, sequential measurements with single GSs, as quantitatively validated in Ref. [39].The multiplexed [Fig.2(c)] and sequential [Fig.2(d)] aberration measurements appear in excellent agreement.To quantify this comparison, we performed a Zernike decomposition on both acquisitions.In Fig. 2(e), we compare the first 25 lowest order modes.Note that tip/tilt ( 1 −1 and  1 1 ) and defocus ( 2 0 ) are omitted here for clarity because they respectively dominate other modes by more than one order of magnitude.Moreover, defocus strongly depends on the residual divergence of the reference plane wave and is therefore non-significant.As can be seen in Fig. 3(e), the agreement between sequential and multiplexed measurements is excellent (RMSE < 0.06λ), showing that several wavefronts generated by multiple GSs undergoing different aberrations can indeed be measured simultaneously in the pupil.In Supplementary S5, we present another experiment where N=5 GSs are multiplexed.We also discuss the experimental gain brought by the iterative DIC approach compared to the direct approach (RMSE reduction by a factor 1.5 to 6).

Correction beyond the isoplanatic patch
Having validated the ability to characterize multiple wavefronts simultaneously, we propose a proof-of-concept experiment showing that the measured wavefronts can be used to correct spatially-varying pupil aberrations.For this purpose, a second sCMOS camera (Panda 4.2, PCO) conjugated to the sample plane by a 4-f system images the fluorescent sample (See Supplementary S4).Figures 3(a When imaging an incoherent object described by its brightness (), the image () obtained in the presence of aberrations can be written as [46]: where  and  are spatial coordinates and (, ) describes the point-spread function (PSF) for a point source at a position s in the FoV.The image is therefore the incoherent sum of contributions affected by aberrations which vary spatially (or angularly).
In most approaches, however, the PSF is assumed to be stationary, i.e. spatially invariant.Under this approximation,  does not depend on the observation direction, and can therefore be simply written (, ) ≈ () all across the field of view.The aberrated image described by Eq. ( 5) then becomes a simple convolution product: If the aberrated wavefront (′) is measured in the pupil plane using a given GS, the stationary PSF () can be estimated using a Fourier transform: () = |{(′)}| 2 , where |{•}| 2 is the square modulus of the Fourier transform.As shown by Eq. ( 6), a simple deconvolution of () by (), using a Richardson-Lucy algorithm (Matlab image processing Toolbox, 15 iterations) then yields a corrected image ().This is rigorously accurate in the direction of the GS used to measure (′) , and the correction remains acceptable in its vicinity (within the isoplanatic patch), but it quickly degrades away from it, beyond a distance Rpatch.This is clearly visible in Fig. 3(c)), where this deconvolution process was applied using a distorted wavefront measured with a single GS (top right).
In this context, the simultaneous measurement of wavefronts coming from multiple GSs described above can provide a better estimation of the angle-dependent PSF (, ), particularly in cases where aberrations decorrelation over time demands simultaneous measurements across the FoV.In the proof-of-concept experiment shown in Fig. 3(d), the wavefronts coming from N=5 GSs, located in different isoplanatic patches, were simultaneously measured to estimate   (′)(j=1 to 5, corresponding wavefronts shown in Fig. 3(b) and deduce the associated PSF   () in each isoplanatic patch.The aberrated image can then be divided into N=5 regions delineated by indicative functions   () ∈ {0,1} which are equal to 1 inside the j th region and 0 elsewhere [46]: The image shown in Fig. 3(d) was corrected using this piecewise approximation, i.e. by performing a deconvolution in each isoplanatic patch with the associated PSF   .The improvement over the single, stationary WF correction  [Fig.3(c)] is significant in the entire FoV, both in terms of resolution and contrast, providing image improvements beyond the isoplanatic patch (See also Supplementary S6).Since the N measurements are performed simultaneously, the method is of particular interest for applications involving time-dependent aberrations.

CONCLUSION AND DISCUSSION
We proposed and demonstrated the use of a thin diffuser to simultaneously acquire multiple aberrated wavefronts by recording a single speckle pattern image.This direct, multi-angle wavefront sensing approach provides deterministic and quantitative measurements.It exploits the large memory effect of thin diffusers as well as the statistical orthogonality of speckle patterns to solve the aliasing problem intrinsic to periodic Hartmann masks.The proposed DIC-based algorithm can successfully reconstruct several (5 or more, see Supplementary S3) angularly-multiplexed wavefronts with high precision (typ.RMSE<) and high resolution (85K phase and intensity pixels), even for large angular distances between wavefronts.When conjugated to the pupil plane of an imaging system, this multi-angle WFS can thus sense, in a single shot, the aberrations of GSs located in different isoplanatic patches of the FoV.We illustrated the potential of this method for the digital correction of aberrations: several PSFs can be estimated simultaneously, and used to accurately deconvolve an aberrated image in multiple isoplanatic patches in order to recover a high-resolution image in the entire FoV.Here, the correction is performed considering a discrete set of GS and aberrated wavefronts, in the corresponding patches, but a better image correction could be obtained by interpolating between patches, and precisely estimating the PSF in each point of the FoV [46].
Multi-angle wavefront sensing has the potential to benefit a wide range of applications, in linear or non-linear microscopy.Coupled to light-sheet excitation, this instantaneous, large FoV digital AO method could be highly valuable to image freely moving animals (e.g swimming Zebrafish or C-Elegans) or large-volume samples [47].Additionally, single-shot aberrations measurement in multiple patches (using multi-spot excitation) has an interesting potential to accelerate deep tissue imaging.In such cases, the isoplanatic patch size can be reduced to a few micrometers only (and can even be reduced down to the size of the PSF in the extreme case of the diffusive regime), thus requiring either extremely long acquisition times or drastically reduced FoV.The high-resolution capability of this multi-angle approach (associated to recently proposed integration algorithms that enable speckle field reconstruction [48]) could be especially relevant to measure highly perturbed wavefronts.Ultimately, the method has the potential to allow the single-shot characterization of the transmission matrix [49].
When integrated into a full AO system (i.e.including hardware wavefront compensation), a single multi-angle WFS could provide instantaneous tomographic-like characterization of a large aberrating volume while significantly reducing the complexity of multi-conjugate AO systems.Given the complexity associated with using multiple compensators, an interesting compensation strategy entails using a single compensator in the pupil to correct the average aberrations in the FoV, followed by a deconvolution process in each patch [20].Another promising strategy involves the use of a multi-angle (or "multi-pupil") compensator in the pupil [7] to correct several (2D or 3D) isoplanatic patches in real-time.
Beyond microscopy, the main fields where this parallel approach should prove most valuable are arguably those in which AO has become indispensable: astronomy and ophthalmology.In both cases, aberrations vary rapidly and isoplanatic patches are relatively small.While we demonstrated multi-angle WFS using fluorescent GSs, it can indeed be applied to other types of GS [50,51] provided that they are mutually incoherent.
While diffuser-based WFS promise cost-efficiency for future applications, the associated speckle patterns have the drawback of spreading the energy over many pixels.This is non-ideal when dealing with low light levels, especially when a large number of wavefronts is multiplexed.The use of an optimized phase mask able to generate orthogonal patterns with energy concentrated on a few pixels, such as a random contour [52,53], offers an interesting opportunity to improve signal-to-noise ratios.
Besides a wide range of applications in adaptive optics, we envision that multi-angle wavefront sensing should open new possibilities in optical diffraction tomography where the 3D refractive index mapping of the sample is usually obtained by sequentially measuring several wavefronts under multiple illumination angles [54,55].Multi-angular WFS should allow multiplexing these measurements to greatly increase the temporal resolution, and even enable single-shot tomography, either in the visible or in the X-ray domain [56].

S1. MULTIPLEXED MEASUREMENT WITH A PERIODIC PATTERN: AMBIGUITY ISSUE
This section aims at illustrating the main advantages of speckle patterns, as compared to periodic patterns, for the angular multiplexing of incoherent wavefronts.To this aim, we first consider a periodic pattern such as those generated by the micro-lens arrays traditionally used in Shack Hartmann wavefront sensors (WFS).In this example, we modeled a 6x6 micro-lens array, each with a micro-pupil size p=150µm and a focal length f=4.5mm.Classically, the conditions avoiding ambiguous overlap with neighboring macropixels gives a maximum dynamic range of  2 ⁄ ≈ 0.0167 ≈ 1° for this WFS.Here, we consider that three wavefronts are multiplexed in a single model experiment [see Fig. S1(a)]: A first wavefront (red) propagates along the optical axis, while the other two wavefronts (green and blue) are tilted by an angle 〈〉 =1.6° larger than the maximum range 1° along the x and y axes.In addition to this tilt, each wavefront is formed by a linear combination of various Zernike modes with random amplitudes.The amount of tilt also contributes to inducing a translation of the pattern that exceeds the dynamic range of the macropixel, leading to a translation overlapping adjacent macropixels.In the multiplexed pattern of Fig. S1(b), we can obviously see in the macropixel marked with a dashed circle that only the red wavefront information is left, while the other two (blue and green) are translated by the tilt to regions ascribed to other macropixels.We can also observe on the zoomed insets in Fig. S1(b) that periodicity induces an ambiguity: the red and blue dots present in the same macropixel as the red one should actually be ascribed to other regions.This problem can be circumvented by designing the Shack-Hartmann WFS in such a way that the spots do not shift more than one lenslet radius, i.e. in such a way than 〈〉 <  2 ⁄ .This can be achieved by increasing respectively the size of the micro-lenses p or decreasing f, but this comes at the price of a highly reduced phase sensitivity or spatial resolution, respectively.
In contrast, the random nature of speckle pattern implies that a speckle pattern, encoding different tilted wavefronts, are orthogonal in each macropixel of the WFS, as illustrated in Fig. 1 of the main text.This property allows multiplexed wavefronts to be distinguished even for large angular distances between tilted wavefronts, and without imposing any constraints on the dynamic range of the WFS. Figure S1(c) shows a comparison of these two wavefronts sensing methods in the phase gradient space.Note however that diffuser-based wavefront sensing imposes a constraint on the maximum detectable angle: the propagation directions of the wavefronts need to be in the range of the diffuser memory effect to allow speckle pattern recognition by intercorrelation.In practice, however, the use of a thin diffuser with a near-infinite memory effect range (i.e. a surface diffuser) alleviates this constraint.

S2. WAVEFRONT REASSIGNMENT CRITERION
We demonstrated that a thin diffuser can be used to simultaneously measure multiple wavefronts with a single acquisition by exploiting the orthogonality of speckle patterns.For a thin-diffuser based WFS [1], a single wavefront can be unambiguously retrieved using a digital-imagecorrelation (DIC) algorithm: it corresponds to the unique maximum of the local cross-correlation between each macro pixel MA of the distorted speckle and the reference speckle R.This cross-correlation peak  ⋆ , however, is not Dirac-like, as it is broadened by both the speckle grain size and the wavefront gradient distribution.
When N wavefronts are measured with a thin diffuser exhibiting a large memory effect [2,3], the multiplexed speckle pattern is an intensity superposition of N replicas of the reference speckle patterns, each having undergone intensity and geometrical transformations caused by the distortion of each wavefront.Therefore, the cross-correlation  ⋆  contains N peaks instead of a single one.In the phase gradient space, these peaks are separated by a distance corresponding to the angular distance 〈〉 separating the wavefronts.Strongly distorted wavefronts, i.e. with strong local wavefront gradients, separated by a small angular distance 〈〉, can therefore lead to situations where peak overlap forbids an unambiguous separation.Note that in these conditions, similar peaks overlap also occurs with Shack-Hartman WFS.
In order to distinguish two wavefronts in a multiplexed speckle, their angular separation 〈〉 must be larger than the lateral width of the cross-correlation peaks (which is driven by both the speckle grain size and the phase gradient distribution, excluding the tilt 〈〉).Here, we numerically investigate this condition.
To model the 1° holographic diffuser (Edmund Optics) used as Hartmann mask, we first characterize its physical parameters by conjugating the diffuser surface to a commercial high-resolution WFS (PHASICS, SID4).Based on the acquired quantitative phase map, we extract a phase correlation width   =46.1µm and a phase standard deviation  ̅ =0.156µm.Using these parameters, we numerically generate a pseudorandom phase mask with realistic statistical properties, as explained and validated in [4].
In order to investigate the condition to properly distinguish wavefronts in a multiplexed measurement, we performed the numerical simulation shown in Fig. S2 here.WF1 is assumed to propagate along the optical axis, and WF2, along the angular direction 〈〉 which varies from 0.4° to 1.8° with steps of 0.05°.The resulting speckle grains displacements   = 〈〉 j +  , can be described as the sum of an average displacement 〈〉  = 〈〉  •  and local displacements  , corresponding to the local phase gradients (excluding the tilt, thus only related to the distortion induced by the Zernike polynomials).For both wavefronts, we apply an amplitude coefficient to the Zernike polynomials which results in a chosen maximum local displacement   (namely,  ,1  =  ,2  = 1 or 5 px).The propagation of the optical field between the diffuser and the camera sensor is calculated using Fresnel diffraction, for a diffuser-camera distance d=3mm and a pixel size of 6.5µm.The angular distance 〈〉 induces a lateral shift 〈〉 •  between the two speckle patterns, as shown in Fig. S2(a).The speckle grain size is measured in Fig. S2(c) as the FWHM of a Gaussian fit on the autocorrelation of the reference speckle pattern R, which yields a size of 3.4 px.Using direct DIC, we evaluate wavefront WF1.
To evaluate the quality of the wavefront reconstruction, we calculate the standard deviation of the difference between the reconstructed and the input wavefront for various values of the angular separation 〈〉.Fig. S2(b) shows this error calculated for  ,1  =  ,2  =1px and 5px (angular broadening applied to both WF1 and WF2).As 〈〉 increases, the measurement error converges to a minimum in both cases.The inflection point, however, depends on the value of    : a strong WF gradient, corresponding to a broad distribution (   = 5 px), requires larger separations (〈〉=1.3°) to be properly reconstructed.Fig. S2(d) shows the conceptual scheme in phase gradient space for the two cases of    .The red disk (1 px) and blue disk (5 px) respectively represent the phase gradient pattern of fixed WF1.The green disk represents WF2, globally tilted along the y axis with a value of 〈〉 • .These phase gradient disks are additionally broadened by a kernel determined by the speckle grain size (orange).For 〈〉 = 0.4°, which corresponds to 〈〉 •  ≈ 3.1 px, the peaks cannot be resolved since the speckle grain size is around 3.4 px.For 〈〉 = 0.7°, the two peaks corresponding to    = 1 px are well separated, but the    = 5 px peaks are not.Finally, for 〈〉 =1.3°, all peaks are well separated, thus allowing proper multiplexed measurements of the wavefronts.
Fig 1(c) and (d): the intensity in the macropixel

Fig. 1 .
Fig. 1.Principle of angularly multiplexed wavefront sensing.(a) Schematic description of the concept, here applied to fluorescence microscopy.N=3 fluorescent guide stars are imaged through an aberrating layer by a microscope objective.This leads to N=3 tilted and aberrated wavefronts in its pupil plane.A wavefront sensor (WFS) based on a thin diffuser allows to (g) multiplex the acquisition of these wavefronts.(b) Basic working principle of the thin diffuser-based WFS: A plane wave illuminating a 1° holographic diffuser generates a reference speckle pattern R on a camera sensor located at a distance d.For a distorted wavefront, measurement of the speckle grain displacements u give access to the phase gradient ( ⊥  ≃  0   ⁄ ).(c,d)The reference speckle pattern R and the multiplexed speckle pattern M measured by the gray-level camera contain information on the three wavefronts.M is the incoherent sum of three distorted speckle patterns related to each wavefront.The pattern measured within a chosen region, the macropixel MA (gray square in d.), can be described as the sum of three patterns of the reference (green, blue and red boxes in c.) each shifted by uj, with j=1,2,3, the index of each wavefront.(e) Owing to the orthogonality between different speckle regions, the cross-correlation of   () with the full reference speckle pattern R reveals N=3 peaks.The peak position gives access to the local speckle grains displacements related to each wavefront, and so, to (f) the phase gradient maps of the wavefronts.(g) A 2D integration step finally allows independent wavefronts reconstructions.
) and3(b)  show the images of the sample without and with an aberrating medium, respectively.

Fig. 2 .
Fig. 2. Experimental validation: single-shot measurements of three multiplexed wavefronts.(a) Acquired reference speckle pattern R (for a plane wave illumination) and multiplexed speckle pattern M. Inset: zoom showing a reduced contrast for the multiplexed case due to the incoherent summation of speckle patterns.(b) Autocorrelation of the reference speckle pattern R, and cross-correlation of the multiplexed speckle pattern M with the reference R. The three peaks in the cross-correlation reveal the number of multiplexed wavefronts and their global tip/tilt ‫ۄۃ‬  which can be determined by measuring the peaks shift ‫ۄۃ(‬  = ‫ۄۃ‬   ) from the center.(c) Local cross correlation allows to reconstruct the three multiplexed wavefronts.Here, the tip/tilt of each wavefront has been subtracted.(d) Comparison with the "classical" sequential method and (e) Zernike decomposition (first 25 modes without piston, tip/tilt) for the three wavefronts.The excellent agreement between both measurements validates the multiplexing method.

Fig. 3 .
Fig. 3. Aberrations deconvolution using multiplexed WFS (proof-of-concept experiment).(a) Non-aberrated image of fluorescent beads used as ground truth.(b) Image acquired in presence of an aberrations medium (c) Aberration correction using a single wavefront measurement (from the GS indicated with a solid circle).The aberration is corrected properly within the isoplanatic patch (indicated by the dashed circle), while the correction fails in the other region.(d) Aberration correction with five wavefronts acquired in a single shot.Each of the multiplexed wavefront is used to correct a certain region within the corresponding isoplanatic patch and a proper image stitch enables to correct the aberrations over a larger FoV.A zoomed is also provided for the 4 cases.

Fig. S1 .
Fig. S1.Ambiguity problem of multiplexed measurement with periodic pattern.(a) Schematic of simultaneously sensing three wavefronts with a micro-lens array.The large tilt angle between the wavefronts induces translation, leading the local focus to shift across macropixels.(b) A numerical simulation is performed to illustrate the ambiguity issue of the multiplexed measurement with micro-lens array.The zoomed macropixel contains the wavefront information from the adjacent pixels because of the translation, while it is not distinguishable to extract the corresponding local phase gradient.(c) Representation of the ambiguity problem in the phase gradient space for the Shack Hartmann WFS and comparison with the diffuser based WFS.
Fig. S2.Numerical simulation: criterion to reassign wavefronts in a multiplexed speckle Pattern M. (a) Two wavefronts corresponding to pure Zernike Polynomials with an average angular distance 〈〉 are multiplexed using a thin diffuser located at a distance d from the camera.The resulting speckle grains displacements   = 〈〉  +  , can be described as the sum of an average displacement 〈〉  = 〈〉  •  and local displacements  , corresponding to the local phase gradients (excluding the tilt, thus only related to the distortion induced by the Zernike polynomials).The two wavefronts are reconstructed from the multiplexed speckle pattern M and compared to the wavefronts used as input.(b) Evolution of the reconstruction error (standard deviation) plotted against the angular distance 〈〉.For both wavefronts, we apply an amplitude coefficient to the Zernike polynomials which results in a chosen maximum local displacement    =  ,1  =  ,2  = 1 ( ) or 5px (blue triangle).The inflection points indicated by arrows correspond to the minimum angular distances 〈〉 allowing wavefront isolation.(c) the reference speckle R and its autocorrelation, from which the speckle grain size is measured.(d) Conceptual description of the cross-correlation peaks in the presence of various phase gradients (radii of the green, blue and red discs) and speckle-induced broadening (orange discs).Broad peaks separated by a small angle 〈〉 (top) overlap, but can be separated (bottom) for larger angular distances.