Imaging with Local Speckle Intensity Correlations: Theory and Practice

Recent advances in computational imaging have significantly expanded our ability to image through scattering layers such as biological tissues by exploiting the auto-correlation properties of captured speckle intensity patterns. However, most experimental demonstrations of this capability focus on the far-field imaging setting, where obscured light sources are very far from the scattering layer. By contrast, medical imaging applications such as fluorescent imaging operate in the near-field imaging setting, where sources are inside the scattering layer. We provide a theoretical and experimental study of the similarities and differences between the two settings, highlighting the increased challenges posed by the near-field setting. We then draw insights from this analysis to develop a new algorithm for imaging through scattering that is tailored to the near-field setting by taking advantage of unique properties of speckle patterns formed under this setting, such as their local support. We present a theoretical analysis of the advantages of our algorithm and perform real experiments in both far-field and near-field configurations, showing an order-of magnitude expansion in both the range and the density of the obscured patterns that can be recovered.


Imaging with Local Speckle Intensity Correlations: Theory and Practice
MARINA ALTERMAN and CHEN BAR, Department of Electrical Engineering, Technion IOANNIS GKIOULEKAS, Robotics Institute, Carnegie Mellon University ANAT LEVIN, Department of Electrical Engineering, Technion Fig. 1. Near-field imaging through scattering. A latent object, comprising mutually incoherent light sources, is seen through a scattering sample (e.g., fluorescent particles inside tissue). The image measured by the camera is degraded due to scattering. Previous approaches suggest that due to speckle statistics, the latent image can be recovered from the auto-correlation of the speckle image. Despite the potential of this idea, many previous experimental demonstrations considered sources a few centimeters behind the scattering layer rather than inside it, as would be desired in a real biological application. Here, we attempt to bring the algorithm from the far into the near field regime, taking advantage of a special characteristic of near field speckles: their local support. In the right part, both the classic full-frame auto-correlation and our approach can successfully recover the latent object when it is composed only of a small number of illuminating points (lower right). However, our approach can also recover a significantly denser object, while the classic algorithm fails (top right).

INTRODUCTION
Developing techniques for imaging through scattering layers, and in particular through layers of biological tissue, is a core challenge of modern imaging. The fundamental difficulty in achieving this objective is the fact that, when an incident wave propagates through such a layer, it interacts with its microstructure multiple times. For example, incident light arising from a single coherent source will, after such a scattering process, result in images showing strong speckle patterns, spread over multiple pixels.
Despite their noise-like appearance, these speckle images have strong statistical properties relating to both the incident wave and the scattering material that produced them. One such property, which is the focus of this article, is the memory effect (ME): This refers to the fact that speckle patterns are correlated and approximately shift-invariant with respect to small tilts in the illumination or viewing angles [Akkermans and Montambaux 2007;Baydoun et al. 2016;Dougherty et al. 1994;Osnabrugge et al. 2017;Freund and Eliyahu 1992;Berkovits and Feng 1994;Fried 1982;Feng et al. 1988]. This property underlies many recent techniques for producing clean images through scattering media [Katz et al. 2014;Bertolotti et al. 2012;Takasaki and Fleischer 2014;Edrei andScarcelli 2016a, 2016b;Hofer et al. 2018;Wu et al. 2017Wu et al. , 2020Wang et al. 2020]. In particular, Katz et al. [2014] have demonstrated that it is possible to recover a clean, scattering-free image of a sparse set of mutually incoherent latent illuminators, observed through a thick scattering layer. Their algorithm works by simply computing the auto-correlation of the observed speckle image and performing phase retrieval [Fienup 1982]. This remarkable imaging capability has strong potential for applications in medical imaging, e.g., for imaging fluorescent cells beneath tissue and for performing noninvasive blood flow analysis.
Unfortunately, realizing this potential in practice remains difficult. One reason for this lies in the fundamental limitations restricting the applicability of memory effect algorithms: For instance, it is well-documented that the strength of the memory effect decays fast as the displacement between the latent illuminators increases Schott et al. 2015]; and we show that it is also negatively affected as the density of the latent illuminators increases. Another reason is that nearly all previous demonstrations of imaging through scattering with the memory effect use illuminators and sensors that are in the far field of the scattering layer, i.e., placed at a large distance from it. (A notable exception is the work of Chang and Wetzstein [2018], which we discuss in detail later in the article.) This setting is reasonable for applications such as non-line-of-sight imaging [Smith et al. 2018;Freund 1990;Katz et al. 2012;Batarseh et al. 2018;Metzler et al. 2020;Viswanath et al. 2018;Boger-Lombard and Katz 2019]. However, it is unrealistic for tissue imaging applications, which typically require near-field imaging conditions. For example, in fluorescent imaging, the fluorescing particles are inside the scattering layer, rather than at a distance from it. Consequently, the current experimental protocol used in research papers for evaluating theory and algorithms about the memory effect is incompatible with how these theory and algorithms would be applied in medical imaging practice.
Our goal in this work is to draw attention to this incompatibility, show that its implications are significant, and propose ways to align research and practice. To this end, we begin with a detailed study of the memory effect in the near-field and far-field settings and highlight the differences between the two settings. In particular, we introduce a new theorem that allows us to draw direct analogies between the two settings and use physically accurate simulations [Bar et al. 2019[Bar et al. , 2020 to both validate and generalize conclusions drawn from this theorem. The findings of our study suggest that, in the near-field setting, memory effect techniques are only practical for scattering layers of modest optical depth, namely, scattering layers whose thickness is only a few mean free paths. In such layers, mid-order scattering is dominant, meaning that light undergoes a small number of scattering events. In tissue, this mid-order scattering regime corresponds to layers that are still well-beyond the maximum penetration depth of standard microscopes. Therefore, memory-effect techniques operating in this regime can be of great practical importance for medical imaging.
Based on this observation, we proceed to investigate how to improve imaging-through-scattering techniques in the mid-order scattering regime. Specifically, we document a property characteristic of speckle patterns arising due to mid-order scattering: The speckle pattern formed on a sensor due to a single latent illuminator is typically much smaller than the sensor. This local support property has not been studied in the past and, as we show, is key for enhancing the performance of imaging-through-scattering techniques in both the far-field and near-field settings.
In particular, we first derive an analytical expression for the signal-to-noise ratio that can be achieved when recovering scattering-free images using the memory effect. Our analysis suggests that there exists an optimal matched filter, corresponding to the local spatial support of speckle patterns, that maximizes this ratio. We then use this theoretical result to motivate and develop a new algorithm for using the memory effect to image through scattering, taking advantage of the local support property. Inspired by ptychography techniques [Rodenburg et al. 2007], our algorithm optimizes for the auto-correlation of overlapping local windows, instead of the full-frame auto-correlation of the entire sensor as in previous algorithms. Our algorithm can be used to improve imaging-through-scattering performance in both the near-field and far-field settings, so long as they operate in the mid-order scattering regime. We demonstrate this improved performance through experiments we perform using both nearfield and far-field imaging prototypes. Our experiments show that, compared with previous auto-correlation approaches [Katz et al. 2014;Chang and Wetzstein 2018], our algorithm results in an order-of-magnitude expansion of both the range and density of independent illuminators that can be recovered.
Implications and future outlook. Together, our theory, simulations, and experiments shed light on a fundamental limit on the performance of memory effect algorithms for imaging through scattering: namely, their ability to recover a clean image of obscured incoherent illuminators deteriorates as the density of the illuminators increases. This is due to the fact that speckle contrast decays quickly when summing speckle patterns from multiple illuminators, as shown, e.g., when comparing the sparse and dense input images in Figure 1. The density of illuminators is a fundamental limit that has not previously received much attention in the literature and is distinct from the better-studied fundamental limit imposed by the memory effect's finite range. That is, recovering a clean image can be unsuccessful in the presence of a large number of independent illuminators, even if all of them are within the memory effect's range.
Our article additionally comprehensively catalogues similarities and differences between the far-field and near-field variants of the imaging-through-scattering problem. In particular, our article demonstrates that the near-field variant of the problem is harder in two ways. First, ME correlation applies to much shorter displacements. Second, as we analyze in this article, due to restrictions on illuminators density, the latent patterns that can be recovered in the near-field setting are sparser and a lot more constrained in terms of their spatial layout than in the far-field setting. Despite these difficulties, our experimental results demonstrate that our algorithm can reconstruct latent near-field patterns of considerable size and density, both an order-of-magnitude larger than what was possible using previous algorithms. These latent patterns already have potential for applications in medical imaging, e.g., sparse blinking fluorescent sources used for STORM localization [Betzig et al. 2006] or sparse cell nuclei observed through microscopes. At the same time, the challenges and limitations analyzed in our article point towards the development of fully robust near-field memory effect algorithms as an important future research direction.
Imaging through scattering. Other techniques for imaging through scattering with coherent illumination use adaptive optics to focus at specific points inside the scattering sample [Rueckel et al. 2006;Vellekoop and Mosk 2007;Yaqoob et al. 2008;Katz et al. 2010;van Putten et al. 2011;Choi et al. 2011;Katz et al. 2012;Vellekoop et al. 2012;Lai et al. 2015;Boniface et al. 2019Boniface et al. , 2020. The main challenge for these techniques is the non-invasive recovery of the aberration correction pattern that the adaptive optics need to apply to achieve focusing. The memory effect can help alleviate this challenge by allowing to adapt a previously recovered pattern to focus at different nearby points [Osnabrugge et al. 2017].
Other approaches for imaging through scattering use incoherent illumination and rely on incoherent intensity models for scattering [Durduran et al. 2010]. Many of these techniques take advantage of additional information available in time-resolved measurements, captured using so-called transient imaging systems [Satat et al. 2015[Satat et al. , 2016[Satat et al. , 2017]. Noteworthy within this category are diffuse optical tomography techniques [Boas et al. 2001;Liu et al. 2020], which use diffusion theory to achieve larger depth penetration, at the cost of reduced resolution compared to coherent techniques.
Improving speckle correlation algorithms. Speckle autocorrelation algorithms for imaging through scattering have recently received increased attention, with several works focusing on improving depth penetration, angular extent, and overall robustness. For example,  proposed to decompose the auto-correlation as a superposition of multiple local auto-correlations, resulting in a threefold improvement in angular extent. Li et al. [2018a] extract spatially varying point spread functions from the speckle correlation using a sequence of illumination patterns. Other techniques [Liao et al. 2019;Chang and Wetzstein 2018] improve robustness by adding sparsity priors on the latent image. Complementary to these techniques are works [Li et al. 2018b;Guo et al. 2020] that use learning-based approaches to allow recovering illuminator patterns wider than the memory effect range. However, these come at the cost of reduced generality-only patterns similar to those available in constrained training datasets (e.g., handwritten digits) can be recovered.
Finally, related to our work are techniques that use ptychography [Rodenburg et al. 2007] to increase the memory effect range [Zhou et al. 2020;Gardner et al. 2019;Li et al. 2019aLi et al. , 2019b. These techniques take multiple images as input, each corresponding to illuminating a different area on the scattering sample. By contrast, our algorithm works using just a single image as input.

PROBLEM SETTING AND BACKGROUND
In this section, we formalize the imaging through scattering problem and clarify the distinction between near-field and farfield imaging conditions. We additionally provide background on speckle statistics, the memory effect, and its use for imaging through scattering.
Imaging geometry. We consider the setup in Figure 2(a). Without loss of generality, we assume that the optical axis of the system is aligned with the z axis. A scattering sample (e.g., tissue layer) of thickness L is positioned between depth planes z min = − L /2, and z max = L /2. We assume that the scattering sample has a width in the x, y dimensions that is much larger than the depth L.
The sample is illuminated by multiple co-planar sources located at depth z i . The light propagates through the scattering sample and generates a speckle pattern, measured by a 2D sensor at depth z v . This can be either a lensless sensor physically located at z v (Figure 2(a)), or an imaging system focused at z v (Figure 2(b)). We denote by i, v the 3D position of illumination and viewing points, and by i x,y , v x,y their x −y restriction to the z =z i and z =z v planes, respectively.
We restrict our discussion to the transmissive setting, where the illuminators and the sensor (or imaging lens) are placed at opposite sides of the scattering sample. Within this setting, selecting z i = z min corresponds to cases where the illuminators are immediately at the back of the sample (e.g., isotropic fluorescent sources inside the sample or confocal illumination focused at that depth). We refer to such cases as near-field configurations (Figure 2(c)). By contrast, selecting z i z min corresponds to cases where the illuminators are placed at a large distance from the sample, as is (a) A sample is illuminated by sources at distance z i behind it. Light propagates through the sample to generate a speckle pattern on a sensor plane at depth z v . (b) The same scene is imaged using a lens focused at the illuminator plane so z v = z i . (c) In near-field configurations, the light is located inside (or at the back face of) the sample rather than far behind it. (d) A typical speckle image obtained on the sensor, as the superposition of scattering from two sources. We mark the illumination and viewing points and the displacements Δ i , Δ v , τ between them. common in prior experimental realizations of imaging through scattering using the memory effect [Katz et al. 2014;Bertolotti et al. 2012;Edrei and Scarcelli 2016a;Hofer et al. 2018;Li et al. 2018b]. We refer to such cases as far-field configurations. At the extreme case, the illuminators can be located at (negative) infinity, which corresponds to illuminating the scene with directional plane waves. When not clear from context, we denote these directional far-field sources and sensors using vectors with a circumflex,î,v, corresponding to their (unit-norm) directions.
Speckle statistics. We denote by u i (v) the complex speckle field generated when light from source i propagates through the scattering sample and is observed at viewing point v. We denote by Consider a scattering sample illuminated by two mutually coherent sources at i 1 , i 2 , and measured at two sensor positions v 1 , v 2 . Then, we define the speckle covariance as: where * denotes complex conjugation, and expectation is taken with respect to multiple realizations of random media with the same statistical properties (e.g., multiple tissue layers of the same type and thickness). Similarly, we can define the intensity covariance 1 : (2) Using classical statistics, it is easy to show that for zero mean fields: We can now use these quantities to formally describe the memory effect (ME) property of speckle fields. We consider the speckle fields u i 1 , u i 2 generated by two nearby illuminators i 1 , i 2 , displaced relative to each other by a vector Δ i ≡ i 2 x,y − i 1 x,y . The ME refers to the fact that u i 1 , u i 2 will be correlated shifted versions of each other. That is, there exists a displacement vector Δ v in the view plane such that In Section 4, we prove that for any given illuminator displacement Δ i , we can compute an optimal view-plane displacement Δ v opt (Δ i ) that maximizes the above correlation.
We show examples of the ME property in Figure 5: Speckle patterns are similar when generated by nearby illuminators, but become different as the illuminator displacement increases. This points to the fact that the ME property holds only for small displacements |Δ i | between light sources, with the correlation decreasing as the displacement increases. To quantify how the ME correlation decays as a function of displacement, it is common to measure the sum of speckle correlations at all sensor pixels. For applications considering field correlations (e.g., adaptive optics for focusing through scattering [Judkewitz et al. 2014;Osnabrugge et al. 2017;Papadopoulos et al. 2016]), this corresponds to: where e ikθ opt (Δ i )v x,y is a phase correction, which is required to make C f (Δ i ) meaningful, since the correlation at individual pixels are complex with possibly different phases. We derive the optimal frequency θ opt (Δ i ) in Appendix A.1. For applications considering intensity correlations, such as the imaging-through-scattering task we study in this article, the ME is instead quantified using: In both Equations (5) and (6), the correlation is evaluated with the optimal view-plane displacement corresponding to the illuminator displacement. We assume that, for a wide homogeneous sample, this correlation depends only on the displacement Δ i , rather than on the exact spatial positions i 1 x,y , i 2 x,y of the illumination sources. We will call the ME range the maximum displacement Δ i for which the correlation C I (Δ i ) remains significant.
Using ME to image through scattering. We now briefly review the method of Katz et al. [2014] for imaging through scattering using the ME. We will be using this method as our baseline throughout the article. Consider a speckle image I generated when a scattering sample is illuminated simultaneously by K mutually incoherent sources i k on one side and imaged by a camera on the other side. The camera measures the incoherent summation of speckle intensities, where I i k (v) denotes the intensity image from the kth source. 2 We denote by S ≡ I (0,0) (v x,y ) the speckle image generated by a source i x,y = (0, 0). If all K sources are within the ME range, then  Figure 3. We note that, if we replace I and S bȳ μ (I ), μ (S ) being signal means, then the same correlation relation holds.
The imageĪ typically exhibits noise-like speckle, making it difficult to discern O. We can, however, consider its auto-correlation, As the intensity values inS are approximately zero-mean independent noise, its auto-correlation is approximately an impulse, S S ≈ δ . Thus, the auto-correlation ofĪ is approximately equal to the auto-correlation of O, Therefore, we can recover O fromĪ Ī using phase retrieval algorithms, e.g., the classical algorithm by Fienup [1982] or more robust strategies [Li et al. 2018b;Guo et al. 2020]. We refer to this procedure as the full-frame auto-correlation algorithm in the rest of the article. A serious shortcoming of the full-frame auto-correlation algorithm is that the range of illuminators it can recover is small, as the maximal displacement between the illuminators |i k 1 − i k 2 | needs to be within the ME range, which is typically very small. Additionally, our article highlights another shortcoming that has received less attention in the literature, namely, that recovery is only possible when the number of illuminators K contributing to Equation (7) is sufficiently small. Our goal in this article is to quantify and compare these constraints in the far-field and near-field settings; as well as to propose a new algorithm for imaging through scattering that can significantly relax these constraints.
We start by exploring some properties of speckle statistics, with a focus towards understanding the relationship between speckle statistics in the near-field and far-field settings. In particular, we explore the effect of the parameters z i and z v , shown in Figure 2, on speckle correlations. The depth z i of the illuminators is relative to the scattering sample controls whether we are in the near-field and far-field settings, and thus investigating how speckle correlations vary as a function of this parameter can help us understand the differences between the two settings. The depth z v of the viewing plane corresponds to experimental choices such as deciding whether to measure speckle using a bare sensor on the front face of the sample, versus using a lens to focus the sensor at a different plane. We will use the observations we make in this section to develop better imaging-through-scattering algorithms in subsequent sections.

Analytic Field Correlation Relationship
In this section, we focus on correlation of complex fields (Equation (5)); we will discuss the correlation of intensity images later. We derive a new technical result that allows converting between field correlations in the near-field and far-field settings. Equation (5), with illuminators and sensors placed at planes z j i , z j v , respectively. Then, correlations measured at different illuminator placements can be related through a displacement scaling as: where in both cases correlation is evaluated at the optimal view-plane displacement given by We provide the proof in Appendix A.1 relying on ideas in Bar et al. [2021]. We will be usinĝ to denote normalized displacements. For small angles, the normalized displacementΔ i is equal to a first-order approximation to the angle two illuminators displaced by Δ i at depth z i form with their midpoint on the plane z = 0 in the middle of the sample; and similarly forΔ v . We visualize these angles in Figure 4(a). Thus, we refer toΔ i ,Δ v as the angular displacements. Claim 1 has three important implications: First, it states that field correlation is a function of only the angular displacement between the illuminators, and not their actual distance. Second, it states that field correlation is invariant to the viewing plane z v , as long as the view-plane displacement is scaled as in Equation (12) to maintain a fixed angle with the scattering sample. Third, it states that we can convert field correlations between the far-field (large z i values) and near-field (small z i values) settings. This may be useful for translating knowledge about the far-field setting into the near-field setting. , highlighting the angular displacementΔ i corresponding to spatial displacements Δ i at different illumination planes. Plots (b)-(e) evaluate ME decay for three distances of the illumination source, and two focus settings. (b) Field correlation, plotted as a function of absolute displacement Δ i . Far sources has much wider ME extant. (c) When plotted as a function of angular displacement, all configurations lead to the exact same decay, regardless of the actual distance of illumination and view planes. That is, field correlation is a function of angular displacementΔ i rather than actual spatial displacements. (d) Intensity correlations as a function of spatial displacement are also much wider for far sources. (e) Intensity correlations as a function of angular displacement are not fully invariant to source distance. (f)-(g) Varying the view plane z v for a fixed illumination plane z i shows that the widest intensity correlation is obtained when the sensor focuses at the illuminator plane. (h) Setup for (i)-(j). (i) Near field intensity correlation as a function of absolute displacement, evaluated for a few material thickness. As the optical depth grows, ME correlation shrinks. For thick slices, the largest displacement at which correlation holds can be too small for any practical usage. (j) Same as (i) but for far sources. In this case, correlation at non zero displacements can be found even for thicker materials.
Claim 1 applies to the correlations of complex fields. However, imaging-through-scattering algorithms use correlations of speckle intensities. Unfortunately, we have not found a similar closed-form relationship for intensity correlations. However, the intuition that such correlations depend mostly on angular displacements still holds. In the following, we explore intensity correlations using numerical simulations as well as real tissue measurements.

Simulation-based Exploration
In this section, we use the physically accurate speckle rendering algorithms of Bar et al. (2019Bar et al. ( , 2020 to perform simulated experiments with two objectives in mind: First, we want to validate the predictions of Claim 1 for field correlations. Second, we want to explore properties of intensity correlations. In our simulations, we use scattering material parameters commonly used in medical imaging to model tissue [Igarashi et al. 2007]. In particular, we simulate a scattering sample with a Henyey-Greenstein phase function of anisotropy parameter д = 0.99, and mean free path MFP = 50μm at wavelength λ = 0.5μm. Except where noted otherwise, we set sample thickness to L = 200μm, resulting in an optical depth OD = 4. We consider three illumination plane settings: (i) z i = z min = −100μm, a near-field configuration where illuminators are placed at the back face of the sample; (ii) z i = −300μm, representing a small gap between the illuminator plane and the sample; (iii) z i = −1cm, which is large enough to correspond to a far-field configuration. The second case is representative of experimental near-field realizations where fluorescent patterns are not really attached to the sample, but rather are projected onto it using a-potentially misfocused-relay lens [Chang and Wetzstein 2018]; or where fluorescent particles need to remain separated from the sample by a thin cover glass. These situations motivate our experiments to understand the effect of the resulting gap on speckle correlations. For each z i setting, we consider two viewing plane settings: (i) z v = z i , corresponding to a sensor that is focused, through an imaging lens, at the illuminator plane; and (ii) z v = −z i , corresponding to a sensor placed at some distance from the front face of the sample. This case is representative of how far-field ME is often measured in practice, by placing a bare sensor far enough from the sample.
Dependence of correlation on angular displacement. In Figure 4(b)-(c), we simulate field correlation values C f as a function of either absolute displacement Δ i , or angular displace-mentΔ i , for the six different configurations of illuminator and view planes z i , z v described above. When parameterizing C f by Δ i , increasing the distance between the sample and illuminator plane increases the ME range. This is due to the fact that increasing this distance while keeping Δ i constant reduces the angular displacement of the sources. As suggested by Claim 1, parameterizing C f byΔ i makes all six configurations identical, showing that field correlations are invariant to the illuminator and view plane locations z i , z v .
In Figure 4(d)-(e), we repeat the above experiments, but this time simulating intensity correlation values C I , which are the actual input to the imaging-through-scattering algorithms we develop later. As in the field case, when we parameterize C I by absolute displacement Δ i , the correlation significantly increases as the distance of the illuminator plane z i increases. Parameterizing C I as a function of angular displacementΔ i brings the six different configurations closer to each other, suggesting that, similarly to field correlation, intensity correlation also depends more strongly on angular displacementΔ i than on absolute displacement Δ i . However, in contrast to the field correlation case, the six configurations are not equivalent, indicating that intensity correlations are not completely independent of absolute displacements. We explain this difference between field and intensity correlations in Appendix A.1; recall in particular that from Equations (5) and (6), Dependence of intensity correlation on view plane. Figure 4(g) demonstrates that, unlike field correlations C f , intensity correlations C I are sensitive to the selection of the viewing plane z v . In particular, we note that intensity correlation is maximized when the viewing plane coincides with the illuminator plane, z v = z i . This is an important practical observation. Often, such experiments use a bare sensor placed on the front face of the sample to measure speckle correlations. Our results suggest that we can increase the measured ME through a simple change in the imaging setup, namely, using a lens positioned in such a way that if there were no scattering sample, the sensor would image the illuminator plane z i . We refer to this imaging configuration as the focused configuration. We validate this observation below using real tissue measurements. ME ranges in practical near-field and far-field configurations. As ME range depends mostly on the angular displacementΔ i , the size of the latent pattern one can handle using ME techniques increases when this pattern is placed further away from the scattering sample. This observation explains why imaging-throughscattering experiments are easier to perform in the far-field than in the near-field. In particular, even though the intensity correlation decay as a function of angular displacement is similar for near-field and far-field configurations, in the near-field case, the corresponding maximal absolute displacement Δ i can become smaller than the wavelength.
To demonstrate this, in Figure 4(i), we simulate intensity correlation values in a near-field configuration, for scattering samples of progressively larger depths, with all other sample parameters remaining the same. This corresponds to increasing the optical depth of the sample, and thus increasing the average number of light scattering events; in turn, this results in a faster decay of intensity correlation as a function of absolute displacement Δ i . The evaluation shows indeed that, for thick samples, near-field ME vanishes for any realistic displacement. By contrast, when repeating the same simulations for a far-field configuration, as in Figure 4(j), we observe that the ME range remains non-negligible even for thick scattering samples. By moving the illuminator plane further away from the sample, we can scale the ME range to cover latent patterns of any size.
These simulations suggest that ME techniques in the near-field setting are only applicable for scattering samples of modest thickness, where mid-order scattering is dominant. Consequently, in the next section, we focus on exploring properties of speckle patterns formed under these conditions that can facilitate the development of imaging-through-scattering algorithms. We note here that the exact sample depth at which near-field ME is non-negligible will vary for different types of scattering materials. In particular, there are significant variations in the material parameters reported as representative of tissue in the literature [Cheong et al. 1990;Tuchin 2000;Igarashi et al. 2007]. We used one set of such parameters for our simulations, but the exact correlation values will be different for other tissue parameters. Therefore, our simulations, and in particular the near-field correlation plots in Figure 4(i), are not intended to precisely predict a tissue depth at which near-field ME vanishes, but rather to support our observation that near-field ME is only non-negligible under modest thicknesses corresponding to mid-order scattering. In practice, for chicken breast tissue, we detect near-field correlations for samples up to 200μm thick. Despite the modest thickness of the samples, images captured through them still contain considerable degradation that can benefit from speckle correlation techniques.

Qualitative Validation Using Real Measurements
Before concluding this section, we present results from real measurements of chicken breast tissue. Data was captured using a near-field experimental imaging setup, described in Section 8.1. These results lend support to the observations presented earlier in this section.
Empirical correlation decay. In Figure 5, we use near field speckle images to demonstrate how intensity correlation C I decays as a function of absolute displacement Δ i . We compute C I (Δ i ) by empirically correlating the captured speckle images. We perform measurements for different placements of the view and illuminator planes. We first compare the intensity correlation plots C I (Δ i ) measured when the sensor is focused at the illuminator plane, versus when it is focused at a different plane. The ME range is wider in the former case, as predicted by our simulations. We then compare the intensity correlations C I (Δ i ) measured for two placements of the illuminator plane: one where it is exactly at the back plane of the sample and another where it is at a distance of 200μm from the back plane. In both of these cases, the sensor is focused at the illuminator plane. In agreement with our earlier observations, our measurements show that placing the illuminators at a distance from the sample results in an increased ME range, when measured as a function of absolute displacement Δ i rather than angular dis-placementΔ i . That is, even a small distance of 200μm, which can occur in experiments due to misfocusing, can significantly impact the ME range.  5. Near-field ME. We present speckle images captured when placing an illuminator at three nearby positions; speckle patterns are shifted versions of each other, demonstrating the ME. The speckle spread is smallest when the illumination source is located exactly at the back plane of the sample and the objective distance is set to focuses on that plane. Focusing the objective on a closer plane (2nd row) results in wider speckles. Computing ME correlation empirically from the captured speckles (lower part), we see that ME correlation holds for larger displacements when the camera is properly focused. We also test the option of moving the light source 200μm further than the back layer while correctly focusing on the illumination plane. In this configuration, higher correlation is measured at wider displacements.
Summary of important observations. To conclude this section, we summarize three important observations we presented: (i) We showed that, in near-field settings, the ME range is non-negligible only for scattering samples of modest optical depth, with a thickness of only a few mean free paths. Such layers are dominated by mid-order scattering. (ii) We showed that it is important to use a focused configuration where the view plane coincides with the illuminator plane to maximize the ME ranges. (iii) We showed that care needs to be taken when placing the illuminators behind the scattering sample, as even a small gap between the two can artificially increase the ME range, even if such a gap is unrealistic for applications such as fluorescent imaging. In the rest of the article, Fig. 6. Local support property. Images of far-field illumination scattering through a chicken breast slice of thickness 250μm, for a focused and a bare sensor. Imaging with a focused sensor reduces the speckle support. We show speckle images for three different positions of the illumination source. The insets demonstrate the ME, namely, that the speckle patterns generated by different illuminators are shifted versions of each other. Lower panel: empirical speckle correlation C I (Δ, τ ) (as defined in Equation (17)) evaluated from the focused data at the second panel. The correlation is displayed as a function of the 2D displacement vector τ for three different choices of illuminator displacementsΔ 0 = (0, 0) • ,Δ 1 = (0.18, 0) • ,Δ 2 = (1.6, 0) • . Due to the modest sample thickness, speckle spread is local, and so is the correlation.
we focus on experimental settings using samples of modest thickness, properly focused sensor, and carefully placed illuminators, as informed by the above three observations.

THE LOCAL SUPPORT PROPERTY
In this section, we document and characterize a property of speckle patterns formed under conditions where mid-order scattering is dominant. In particular, in Figures 5 and 6, we show speckle patterns measured through a tissue layer of modest thickness, using our near-field and far-field imaging setups. We observe that the speckle patterns have local support, much smaller than the full extent of the sensor. Local supports are prevalent in biological tissue samples of modest thickness, as the phase functions characterizing these samples are strongly forward-scattering-their average cosine is typically д > 0.95. Given such phase functions, light entering the sample with directionî will, after undergoing a small number of scattering events, spread primarily towards outgoing directionsv ≈î. As we will show in the next section, the local support property of speckle patterns due to mid-order scattering is key for improving imaging-through-scattering algorithms based on the ME. Given its importance, we use this section to study this property in more detail.
Effect of focusing on speckle support. In Section 4, we showed that using a lens to focus the sensor at the same plane as the illuminators is important for increasing the ME range. Figures 5 and 6 demonstrate an additional advantage of using this focused configuration: In both the far-field and near-field settings, focusing decreases the speckle spread. Intuitively, the wider spread of nonfocused configurations can be explained by the fact that the support of the scattered field is convolved with a defocus blur kernel.
Far-field versus near-field speckle patterns. Speckle patterns formed under far-field and near-field settings both exhibit the local support property. However, far-field patterns, such as those shown in Figure 6, include many more speckle features compared to nearfield ones, such as those shown in Figure 5. We will see in subsequent sections that this difference makes imaging through scattering in the near-field setting more challenging than in the far-field setting.
Extended parameterization of intensity correlation. For the rest of this article, we restrict the discussion to focused configurations, where z v = z i . In this case, we get from Equation (12) that the displacement on the illuminator plane and the corresponding optimal displacement on the view plane are equal. Thus, we simplify notation using Δ ≡ Δ i = Δ v opt (Δ i ) for both displacements, leading to: We note that the definition of the intensity correlation C I (Δ) in Equation (6) treats all pixels v 1 x,y equally and does not consider the location of the pixel v 1 x,y relative to the illuminator location i 1 x,y . This is due to the fact that most prior literature on the ME focuses on settings where speckle patterns cover the entire sensor plane (e.g., cases where high-order scattering is dominant). Consequently, correlation does not vary significantly at different locations v 1 x,y on the sensor, and it is sufficient to analyze how correlation C I decays as a function of the displacement Δ alone. By contrast, in our setting, the local support property implies that speckle patterns generated by an illuminator at location i 1 x,y are concentrated at pixels in locations v 1 x,y adjacent to i 1 x,y . This suggests that correlation can vary at different locations v 1 x,y on the sensor: For example, as we move away from i 1 x,y , less light is measured, and we expect correlation to be reduced. To characterize this effect of the local support property, we will modify the definition of the intensity correlation C I so it takes as input the displacement between the illuminator and the pixel rather than only the displacement between the two sources. 3 To this end, we denote the 2D displacement between the illuminator and pixel locations as: We visualize both τ and Δ in Figure 2(d). For illuminator and pixel pairs satisfying Equation (14), it follows that: Then, we define the intensity correlation for illumination and pixel pairs satisfying both the illuminators displacement relation of Equation (14) and the illuminator-pixel displacement relation of Equation (16): where the intensity covariance C I (I i 1 (v 1 ), I i 2 (v 2 )) was defined in Equation (2). We note that we can relate this definition to the definition of intensity correlation as a function of Δ alone in Equation (6) through the equation C I (Δ) = τ C I (Δ, τ ).
To demonstrate empirically the importance of parameterizing the intensity correlation C I as a function of both Δ and τ , we use speckle images captured from chicken breast tissue samples with a far-field experimental imaging setup described in Section 8.1 below. In Figure 6, we show the intensity correlation C I (Δ, τ ) computed from the image measurements. We can observe that C I varies significantly as a function of both Δ and τ , and in particular that it quickly decays as the distance τ between the illuminator and pixel location increases. A schematic of the displacments Δ, τ is visualized in Figure 2(d). To our knowledge, the local support property and its effect on intensity correlations have not previously been used for imaging-through-scattering applications. In the next section, we provide a theoretical justification for using this property; then, in Section 7, we use it to develop an improved algorithm for imaging through scattering.

SIGNAL-TO-NOISE RATIO ANALYSIS
Previous studies of the full-frame speckle auto-correlation algorithm that we described in Section 3, for example by , have focused on how the limited ME range constrains the size of the latent illuminator pattern that can be recovered. In this section, we study a second constraint on the recoverable latent illuminator pattern that has received little attention in the literature (see limited discussion in the supplement of Katz et al. [2014]): the fact that reconstruction is usually successful only when the number of different illuminators K in Equation (7) is sufficiently small. When a large number of incoherent sources contribute to the measured intensity image, speckle contrast decays and correlation becomes noisier. For example, this difference in speckle contrast is noticeable when comparing the sparse and dense inputs of Figure 1. We show that, by taking advantage of the local support property we described in the previous section, we can significantly increase the signal-to-noise ratio (SNR) of the correlation, and consequently, the density of illuminators we can recover. We note that the reconstruction algorithm by Katz et al. [2014] has two parts: first computing speckle correlation, and then performing phase retrieval. The focus of our analysis is on the first part, the SNR at which correlation can be computed. Even though we expect that the performance of the phase retrieval part will also improve as the noise characteristics of its input improve, a detailed analysis of phase retrieval convergence is beyond the scope of this work.
We are given a speckle image I formed as in Equation (7) and want to examine whether illuminators i 1 , i 2 contributed to its formation. Denoting the illuminator displacement Δ = i 2 x,y − i 1 x,y as in Equation (14), we can multiply the zero-mean speckle imagē I with its shifted copy, then form a correlation estimate using weighted pixel averaging: We expect c emp to have a large value when an illuminator pair i 1 x,y , i 1 x,y + Δ exists and a value close to zero otherwise. When computing full-frame auto-correlation, as in Equation (10), the spatial weights w are uniform over the entire image I . However, if we know that the speckle patterns have local support, we can consider setting non-zero weights w only in a window around i 1 x,y , rather than in the entire image. We state a new technical result showing this can drastically improve SNR, and derive the optimal weighting strategy.
To formulate this result, we denote by P the number of sensor pixels and by F the number of speckle features in the image, where a feature refers to a diffraction-limited speckle spot. We have F ≤ P, where a gap F < P happens for two possible reasons: First, depending on the aperture, a diffraction limitted feature can spread over more than a single pixel; and second, even for single pixel features, the combined speckles from all illuminators may not cover the entire sensor. Additionally, we denote by K the number of illuminators. 4 Using these notations, we define the density of independent illuminators as Using the density definition, we state the following claim: Claim 2. The signal-to-noise ratio of the estimator of Equation (18) and is maximized by the matched filter w (Δ, τ ) ≡ C I (Δ, τ ), reaching We provide the proof in Appendix A.2. Algorithmically, using the matched filter requires averaging only within the local image window where we expect to have speckle from illuminator i 1 , and not within the entire sensor as in the full-frame auto-correlation algorithm. We provide algorithmic details in Section 7.
Implications. Claim 2 suggests that using the matched filter instead of uniform summation over the image can significantly improve SNR. To qualitatively characterize this improvement, we assume for simplicity that the size of a speckle feature is one pixel, the support of the speckle pattern due to one illuminator is N pixels, and all sensor pixels receive light from at least one illuminator, so F = P. Suppose also that C I (Δ, τ ) = 1 inside the support and 0 otherwise. From Claim 2, the matched and uniform filters achieve SNRs of SNR matched = 1 /(α 2 N ), SNR uniform = 1 /(α 2 P ).
4 In the supplement of Katz et al. [2014], density is defined as the area of high emission in the target, divided by the area of a diffraction limited spot. Fig. 7. SNR gain. We visualize the correlation c emp (i 1 x,y , Δ) over a 1D line (highlighted on O ). As this line includes three different illuminators, we expect to detect high correlations at three displacements. Correlations resulting from the matched filter (red curves) are less noisy than correlations from full-frame averaging (blue curves). We simulated observed images I due to the same illuminator arrangement O and three different speckle support sizes, visualized as insets at the top right corner of the corresponding observed images. As predicted by our theory, SNR improves for medium speckle support size, but decays for very small and very large support sizes. Therefore, using the matched filter versus full-frame averaging improves SNR by N /P. When the sensor size is a few megapixels and the speckle support of each illuminator is only 100 × 100 pixels, this translates into an SNR improvement of two orders of magnitude.
We can also use Claim 2 to understand what illuminator density α we can expect to reliably detect. Suppose that for good detection we seek an SNR larger than a threshold R. As before, assume speckle features are single-pixel wide and the full sensor is covered by speckles. The matched and uniform filters lead to different upper bounds on the recoverable density: Selecting, e.g., R = 100 as a threshold for reliable detection, and given P = F = 10 6 , we find that a uniform filter can reliably detect only one illuminator per 10 4 pixels. By contrast, if the speckle support from one illuminator includes N = 10 4 pixels (e.g., a 100 × 100 support), then the matched filter can reliably detect one illuminator per 10 3 pixels; for N = 10 2 , this becomes one illuminator per 10 2 pixels. We note that the above limits on SNR as a function of illuminator density hold even if all illuminators are within the ME range. Therefore, the constrained density is a fundamental limitation of the full-frame auto-correlation method for which there is limited discussion in the literature.
Visualizing the SNR gain. In Figure 7, we use a synthetic example to visualize the SNR gain achieved using the matched filter. We generate speckle images using the idealized formula I = S O, so all pixels are inside the ME range. We use three speckle patterns S of different support for the same latent image O. We then compute c emp (i 1 x,y , Δ) as in Equation (18), for the i 1 x,y point marked in Figure 7(a). For simplicity, we vary Δ only over one horizontal line marked in the figure. As the line contains only three illuminators, we ideally expect high correlation only for three translation values. We observe that, in agreement with Claim 2, uniform averaging produces significantly noisier correlations than the matched filter (compare the red and blue curves in Figure 7). In Figure 7(c), where we use a medium support N , the matched filter produces sharp correlation peaks at the correct displacements. In Figure 7(d), we increase the support N , and the matched-filter correlation becomes noisier; this agrees with Equation (22), which states that SNR matched decays as the support N increases.
We now consider Figure 7(b), where we use a small support N , and some sensor pixels do not receive light. The correlation becomes worse than that obtained with the medium support in Figure 7(c). To understand this, we note that when transitioning from the medium support in Figure 7(c) to the large support in Figure 7(d), the illuminator density α as defined in Equation (19) remains the same; as in both cases, the entire sensor is covered by speckle features so the feature number F remains the same. By contrast, in the case of small support in Figure 7(b), the density α is higher; this is because the number of sensor pixels covered by speckle, and thus the number of speckle features F , are both reduced. Consequently, from Equation (22), the SNR is also reduced. This example is important for understanding the near-field setting, where typically speckle support sizes are small and the speckles from all illuminators do not cover the entire sensor.

OPTIMIZING USING LOCAL SUPPORT
The previous section provides theoretical justification for using the matched filter, rather than uniform weights. In this section, we develop a new algorithm for imaging through scattering that explicitly takes into account the local support.
We begin by noting that, in many practical cases, we cannot measure the exact correlation C I (Δ, τ ), and thus cannot compute the exact matched filter of Equation (21). Instead, our algorithm will approximate it using two binary thresholds T τ ,T Δ , assuming that speckles from one illuminator are spread over pixels in a window of size T τ around it, and that ME correlation holds for displacements |Δ| < T Δ . The thresholds T τ ,T Δ are free parameters that we can fine-tune to improve reconstruction quality. We show in Appendix A.5 that performance is robust to their exact values. As we discuss below, our algorithm offers improved performance compared to the baseline full-frame auto-correlation algorithm in Fig. 8. Local window selection for optimization. We consider local subwindows w τ (light green and cyan frames) whose support is equivalent to the speckle support size. Each such window is correlated with a wider window w Δ (yellow and blue frames) around it, whose support is equivalent to the ME range. As speckle inside window w τ can arise from a source outside w τ , O w τ O w Δ may not match I w τ I w Δ . To overcome this, we use an extended non-binary sub-windoww 2τ = w τ w τ for O , whose support is indicated by dashed lines.
situations where T τ < T Δ , namely, when the support from one illuminator is lower than the ME range. For thick scattering slices, where high-order scattering is dominant, this relationship does not hold, and our approach reduces to the baseline full-frame autocorrelation algorithm of Equation (10).
Our algorithm searches for a latent image O such that the autocorrelation in its local windows will match the auto-correlation in the local windows of the input image I . We define w Δ and w τ to be binary windows with support T Δ ,T τ , respectively, andw 2τ = w τ w τ -note that, from its definition,w 2τ is non-binary. Then, we recover O by solving the optimization problem: denote windows of a given size cropped from the input and latent images, centered around the jth pixel.
Equation (24) uses windows of three different sizes, and we use Figure 8 to visualize their different roles: Each w τ j is a small window whose support is equivalent to the expected support size of the speckle pattern due to a single illuminator. w Δ j is a larger window around it, corresponding to the maximal displacement T Δ for which we expect to find correlation, as dictated by the ME range. If the windows w τ j , w Δ j are centered around pixel i x,y , then the Δ entry of the correlationĪ w τ Equation (18), where the matched filter is approximated by the binary window w τ .
We note, additionally, that the window cropped from O should be wider than that from I . This is because speckle at a certain pixel can arise from an illuminator within a window around it. For example, in Figure 8, no illuminator is located inside the cyan subwindow of O, but a part of the speckle pattern is contained within the corresponding cyan subwindow of I . As a result O w τ is easy to prove that this can be addressed using the larger, nonbinary windoww 2τ in the latent image, indicated in Figure 8 using dashed lines: In this case, Ow2τ The motivation for the cost of Equation (24) is that, even if two illuminators in the latent pattern O are at a distance larger than the ME range T Δ , they can be recovered if there exists a sequence of illuminators between them, where each two consecutive illuminators in the sequence are separated by a distance smaller than T Δ . For example, in Figure 8, the illuminators outside the yellow and cyan w Δ windows are recovered, thanks to the intermediate illuminators.
As pre-processing for our optimization procedure, we form an approximation for the zero mean speckle signal defined in Equation (8) by subtracting the local mean of each window: where G is a Gaussian blur filter. The optimization problem in Equation (24) is no longer a phase retrieval problem as in standard full-frame auto-correlation algorithms. We minimize it using the ADAM gradient-based optimizer [Kingma and Ba 2014]. Gradient evaluation is described in Appendix A.3 and reduces to a sequence of convolution operations that can be performed efficiently, e.g., using a GPU-based fast Fourier transform. For initialization, we set the latent image to random noise; we have observed empirically that the optimization is fairly insensitive to initialization. Finally, we note that even though we could place a window w j around every pixel of I , the empirical correlation is insensitive to small displacements of the central pixel j. Therefore, in practice, we consider windows only at strides T τ /2, which helps reduce computational complexity.
To conclude this section, we note that the optimization problem of Equation (24) is similar to ptychography algorithms [Rodenburg et al. 2007]. However, we emphasize that previous ptychographic approaches for extending the ME range recover the latent illuminators from multiple image measurements, captured by illuminating different areas on the scattering sample [Zhou et al. 2020;Gardner et al. 2019;Li et al. 2019aLi et al. , 2019bShekel and Katz Fig. 9. Hardware setup. Top row: schematic of far-field setup, demonstrating two illumination configurations used in experimental setup. The first translates a single point source (fibered laser). The second uses an LED source with masked area of emission. Second row: schematic of near-field setup. Lower panels visualize our hardware lab setup. 2020]. By contrast, our algorithm recovers the latent illuminators from a single shot.

EXPERIMENTS
We begin by evaluating our algorithm in the far-field setting, demonstrating that even in this setting it provides an order-ofmagnitude extension of both range and density of illuminators that can be recovered, compared to the full-frame auto-correlation algorithm. We then proceed to show experiments in near-field setting, demonstrating again significant improvement over previous approaches. We discuss the challenges of the near-field setting and show that they are in agreement with our theoretical analysis.

Experimental Setup
We built two hardware setups shown in Figure 9, implementing near-field and far-field imaging configurations.
For the near-field setup, we use a tube lens and an objective lens to focus a point source (the output of a single mode fiber connected to a 632nm laser) into a point source at the back side of a Fig. 10. Local versus global auto-correlation. The orientation of the auto-correlation evaluated in three different local windows of the image matches the orientation of the arc in the corresponding region of the latent image. By contrast, the auto-correlation of the full frame is much nosier and decays for large displacements due to limited ME. scattering sample. To image the sample, we use a camera placed at the opposite side of the sample, similarly equipped with a tube lens and an objective lens. As we discussed in Section 4, verifying that the point source is focused exactly on the back face of the sample is important for ensuring that our experiments are representative of realistic scenarios where illuminators are located inside the tissue rather than beyond it. We confirmed that using a second, control camera. Additionally, as discussed in Section 4, we use a focused configuration where the camera is also focused at the same plane as the point source, at the back face of the sample. We verified the camera focusing by scanning its objective lens along theẑ axis, capturing a focal stack, and selecting the position where speckle support size is smallest and ME range is largest. We used Nikon N20X-PF objectives with NA = 0.5 and ×20 magnification, and Thorlabs TTL200 tube lenses.
For the far-field setup, instead of placing the source at a large distance from the sample, we placed it at the Fourier plane of a lens, creating fully directional illumination. This configuration is equivalent to a point source infinitely far away from the sample and allows for better light efficiency. Likewise, we use a camera on the opposite side of the sample, focused at infinity. This setup corresponds to a 4F system around the sample, which we implement using two macro lenses (Nikon 105mm f/2.8D).
For scattering samples, we use slices of chicken breast of thicknesses ranging between 100 − 400μm. We measure the thickness of samples by placing them between two microscope slides of known dimensions and using a caliper to measure the total thickness.
In both the near-field and far-field setups, we translate the laser point source at different locations behind the sample, capturing different images I i k at each location. We then sum these intensity images to form I = k I i k , simulating the input from multiple mutually incoherent sources. Having access to the individual I i k images is useful for analyzing various algorithmic tradeoffs. We also use a second, single-shot setup, consisting of a binary mask illuminated by a wide-area, spatially incoherent 625nm LED.

Far-field Experiments
Local auto-correlation. In Figure 10, we visualize the different structure of local and global auto-correlations. Computing autocorrelation at small subwindows of the speckle image reveals the local orientation of the arc in the latent image. By contrast, when computing the auto-correlation of the full frame, the correlation is considerably noisier even for small displacements. Correlations between far illuminators are even harder to detect due to the limited ME range.
Range and density. As discussed in Section 6, a fundamental limitation of imaging-through scattering algorithms is the density of illuminators they can recover. To demonstrate this, in Figure 11, we compare recovery results for illuminator patterns of the same range and layout, but at different densities. In each case, we display the densest subset at which the full-frame and our local autocorrelation algorithms successfully recovered the latent pattern, with our local approach often handling order-of-magnitude larger densities. We captured the data by imaging speckle patterns created by individual point sources placed at different locations and summing the speckle images in post-processing, allowing us to form test images at any density of interest. Details on the fullframe phase-retrieval algorithm we used, as well as a comparison to the sparse approach of Chang and Wetzstein [2018] are provided in Section 8.4.
As we have access to the speckle images generated at each illuminator location, we can compute the decay of ME across the frame. Denoting by I i 1 , I i 2 the individual speckle images from illuminator locations i 1 , i 2 , and setting Δ 1,2 = i 2 x,y − i 1 x,y , we evaluate: We plot this correlation at the right of Figure 12, as a function of |Δ|. We note that for the smile pattern, which was captured with a thin tissue layer, the ME range covers the entire frame (empirical correlation does not decrease below 0.8 even for the widest displacement). Even under these favorable conditions for the fullframe auto-correlation algorithm, our local algorithm recovers a denser set of illuminators.
In Figure 12 and additionally in Figure 24 of Appendix A.4, we demonstrate the increased range that our algorithm provides. To achieve this, we select a few local subwindows from the patterns in Figure 11 and display the maximal window for which the fullframe auto-correlation was successful-each pair of subwindows demonstrates a small window with reasonable reconstruction and a slightly bigger one where reconstruction already failed. Overall, in Figure 11 our local algorithm successfully handles patterns that Fig. 11. Comparison of our local and the full-frame auto-correlation algorithms. For each example, we show the densest arrangement of illuminators for which the full-frame auto-correlation algorithm succeeded. In the top example, our algorithm successfully recovered ×32 more illuminators. Even in the lower example where the ME extends over the entire frame (see correlation plots in Figure 12), our local approach outperforms the full-frame one. The tissue thickness of each example, from top to bottom, is 330μm, 340μm, and 200μm. Fig. 12. Full-frame auto-correlation algorithm applied to small crops of the patterns in Figure 11. The yellow and cyan sub-windows demonstrate areas where reconstruction roughly succeeds, and the magenta and green ones a slightly larger window where reconstruction fails. To the right, we plot correlation as a function of displacement length |Δ |, as measured for the corresponding tissue slice. Tissue thickness from top to bottom, are 340μm and 200μm. Fig. 13. Classical setup illustration. In the classical full-frame autocorrelation setup, the latent image is usually significantly smaller than the sensor width, and the number of speckle features it includes is about ×10 4 higher than the number of independent illuminators.
are an order of magnitude wider than the maximal patterns recovered by the full-frame approach in Figure 12. There is no inherent limit preventing us from handling an even larger range, except that in the specific experimental setup we used, increasing the range would exceed the aperture width.
For the small images of Figure 12, our algorithm is equivalent to the full-frame auto-correlation algorithm, as the images are small enough that they do not fit more than one w τ window. While both approaches fail on small images (Figure 12), the local approach is successful when applied to larger images, where the full-frame algorithm still fails ( Figure 11). This is because, when considering a larger image, our algorithm computes correlation with speckle patterns at other parts of the frame, providing additional constraints.
Contrasting with classical setup. The patterns recovered in our implementation are very different from the ones used in previous full-frame auto-correlation implementations [Katz et al. 2014]. The patterns in Figure 11 included about 10 3 illuminators spread nearuniformly across the area of a 2-megapixel sensor. By contrast, Figure 13 shows a typical input for previous full-frame implementations, where the target pattern is concentrated within a small area of about 100 × 100 pixels. Yet, the speckle support is much larger, covering the entire sensor. To achieve this wide speckle spread, previous implementations either imaged the sample with a lensless sensor rather than a focused one or used scattering layers that are thicker or have wider phase functions. In synthetic simulations of such a full-frame setup, the phase retrieval algorithm by Fienup [1982] usually fails if more than 100 sources are included; The number of sources can slightly increase with a better phase retrieval approach.
Finally, it is worth noting that, as illuminator density increases, the local approach eventually fails as well. We include an example of such a failure case in Figure 22 of Appendix A.4.
LED illumination. In Figure 14 and in Appendix A.4, we show reconstructions from the single-shot setup of Figure 9, where the entire area of a target mask is illuminated by spatially incoherent LED light. The main challenge in this case arises from the fact that the illumination is no longer purely monochromatic: Different wavelengths are diffracted in slightly different angles, blurring  Despite the seemingly small degradation, the full-frame approach fails unless provided an input composed of a considerably sparser set of illuminators. To the right, we plot correlation decay as measured using images of individual sources through the corresponding tissue slice. speckle contrast. In the mid-order scattering examples that we are considering, this is mitigated by the fact that the speckle support size is limited, meaning that speckle patterns are less affected by blur. To reduce this effect, we placed a 10nm band-pass filter at the LED output.

Near-field Experiments
In Figure 1 and Figures 15 and 16, we show reconstruction results from our near-field setup. The first set of examples (Figures 1  and 15) come from thinner tissue layers (L = 100 − 150μm), for which scattering is modest and the latent pattern may be Fig. 16. More near-field comparisons. We now use a thicker tissue slice (L = 200μm). The lower rows zoom on two of the digits in the second row, demonstrating reconstruction at a few different densities. Reconstruction degrades as density increases. The full-frame approach was successful at a considerably lower density than our local approach.
recognizable from the degraded input image. Our algorithm still improves the pattern quality significantly and reconstructs fine details obscured by the speckle. However, even this modest degradation is challenging for the full-frame approach, which fails unless applied on a significantly sparser set of illuminators (last rows of Figures 1 and 15).
For the second set of examples in Figure 16, the sample thickness is larger at L ≈ 200μm, and thus the degradation is stronger. This reconstruction task is more challenging for two reasons: First, the ME range is limited, as seen by the correlation curve at the right of Figure 16. Second, the illuminator density is large. Our local algorithm outperforms the full-frame auto-correlation approach, but its reconstruction is not free of artifacts either. In the lower rows, we zoom on two of the digit patterns, showing reconstructions for different illuminator densities. Our local approach performs worse as density increases, but still outperforms the full-frame approach, which is successful only at considerably smaller densities.
To demonstrate the difference between the far-field and nearfield settings, in Figure 17, we compare reconstructions of a line pattern of the exact same length, from far-field and near-field measurements. We can see that the far-field speckles cover a much larger area of the sensor than the near-field ones. These additional speckle-covered pixels help improve the SNR of the far-field correlation estimates. The explanation for this difference is the same as for the difference between the small and medium supports in Figure 7(b) and (c): The near-field image in Figure 17, due to the many pixels that do not receive light, has a larger effective illuminator density than the far-field image. This larger density results in worse correlation SNR, and thus reduced reconstruction quality. As another way to see this, we show in Figure 18 speckle patterns from a single point source captured under near-field and far-field conditions and the corresponding auto-correlationS S . Both auto-correlations resemble an impulse, up to noise. However, we observe that the near-field speckle image includes much fewer speckle features than the far-field one. Consequently, there is more noise in the near-field auto-correlation image than in the far-field one. This increased noise results in the reduced quality of the nearfield reconstructions.
Previous near-field implementations. The only reported attempt to apply speckle auto-correlation techniques in the near-field setting that we are aware of is by Chang and Wetzstein [2018]. A direct comparison between our results and theirs is not possible, as the two sets of experiments use very different scattering samples. In particular, the speckle images captured by Chang and Wetzstein [2018] are similar to the ones used by Katz et al. [2014], comprising a target pattern concentrated with a small sensor area, Fig. 17. Comparison of near-field and far-field settings for a line pattern. Even though the ground truth pattern is similar, the near-field speckle pattern has a much smaller support size. As more speckle pixels are provided, the far-field correlation is less noisy, improving reconstruction quality.
producing non-localized speckle patterns covering the entire sensor. The ME range reported by Chang and Wetzstein [2018] was approximately equal to a 10μm displacement. Accordingly, their experiments recovered illuminator patterns of size 10μm × 10μm. By contrast, the size of the near-field illuminator patterns we recovered scales up to 65μm × 65μm. Below, we additionally compare with their robust phase-retrieval algorithm, using measurements captured with our imaging setups.

Comparison to Alternative Algorithms and Limitations
In all of the previous figures, the reconstructions for the fullframe approach are achieved using an ADAM optimization procedure [Kingma and Ba 2014] with non-negativity constraints. This is analogous to Wirtinger flow optimization for phase retrieval [Candes et al. 2014;Chakravarthula et al. 2019]. We found that this approach works better than the classical optimization algorithm by Fienup [1982]. We show comparisons between these two algorithms in Figure 19, using sparse and dense sets of illuminators for both the far-field and near-field examples from Figures 11 and 16. For these examples the Fienup-based full-frame variant (Figure 19(b)) did not converge, whereas the ADAM-based full frame variant converged on the sparse set and failed to converge on the dense one. We additionally attempted to optimize the full-frame approach with the ADMM-based phase retrieval algorithm of Chang and Wetzstein [2018], which uses an L 1 regularization term. As shown in Figure 19(c), this performed better than the Fienup-based variant, but provided results very similar to the ADAM-based variant for most examples. We believe this is because the non-negativity constraints we enforce during ADAM optimization already leads to sparse solutions. We note in Figure 19 that all variants of the full-frame approach fail as we further increase the illuminator density. By contrast, our local correlation approach is successful in the higher density case, as shown in Figures 11 and 16. The auto-correlation of speckle images due to a single sourceS S resembles an impulse plus noise, in both the near-field and far-field cases. However, in the near-field case, this auto-correlation is noisier, as fewer speckle features are averaged. The single-source speckle imagesS are shown in the insets.
Given the local extent of the ME, another option one may consider is cropping local windows from the full speckle image, running the full-frame auto-correlation approach on each local window, and then seaming the individual local solutions. However, as shown in Figure 19(e), the independent solutions are rarely consistent and the seamed result has strong artifacts. We also note that the solution of the full-frame approach in each local window is only defined up to an arbitrary flip or shift. For the result in Figure 19(e), we favored this algorithm by flipping and shifting each window to best match the groundtruth. Even under this simplification, this algorithm is inferior to our approach that jointly optimizes all local windows.
Finally, in Figure 19(f), we compare against the approach of Wang et al. [2019] for extending the range of imaging-throughscattering algorithms. Their model assumes the latent image O can be decomposed into two parts O 1 , O 2 of a smaller extent, and the ME applies in each window separately. Mathematically, they model the image formation as I = O 1 S 1 + O 2 S 2 , where S 1 and S 2 are the speckle patterns from a single illuminator in each region, and which are assumed to be decorrelated, S 1 S 2 = 0. They then try to simultaneously solve for two smaller support images O 1 , O 2 satisfying O 1 O 1 + O 2 O 2 = I I . As seen in Figure 19, this approach was successful on the far-field sparse examples, but failed on the denser and near-field examples. The reason for this is that, in our examples, the correlation decays gradually, and thus the assumption by  that speckle patterns at different parts of the image are completely decorrelated does not hold. Additionally, the range of the latent illuminators is much larger than twice the ME range.
Runtime. Compared to the full-frame auto-correlation approach, one disadvantage of our algorithm is increased computational cost. In particular, our unoptimized Matlab implementation, when running on an NVIDIA Quadro RTX 8000 GPU, requires a few hours to converge for each of the results in this section. Gradient evaluations make up the bulk of this runtime. Gradient evaluations are essentially a sequence of convolution operations, which we implement using Matlab's GPU-based fast Fourier transform function. These gradient evaluations can potentially be accelerated using more sophisticated GPU-based convolution libraries. We also note that our optimization procedure only needs to be run once, as its  Figures 11 and 16, we evaluate a few alternative strategies. (b)-(d) show different full-frame phase retrieval approaches. Classical Fienup optimization [Fienup 1982] is rather noisesensitive. Chang and Wetzstein [2018] proposed a better algorithm introducing a sparse prior and ADMM optimization. In this work, we used gradient decent update with ADAM step selection size, with similar results. These algorithms can solve the full frame phase retrieval on a sparse subset of the sources, but fail on a denser one. In contrast, our local cost led to much better results, as presented in Figures 11 and 16. (e) Solving the standard phase retrieval problem on independent local windows and seaming the results in post-processing leads to noticeable artifact. (f) Seeing beyond the ME range ] by decomposing the speckle auto-correlation into two independent parts is sometimes successful on the sparse data, yet fails on the dense one.
results are insensitive to initialization. By contrast, iterative phase retrieval algorithms such as the algorithm by Fienup [1982] typically require multiple runs with different initializations.

DISCUSSION
We provided a comprehensive study of algorithms using the speckle ME to image through scattering. Using theory, simulations, and real experiments, we investigated the inherent limits of these approaches. In particular, we explored whether these approaches can be applied to practical biomedical imaging scenarios, where illumination sources are located inside, rather than far behind, a scattering sample such as a tissue layer.
We reported the following important findings of our study: First, we showed that the ME is affected by the angular difference between illumination sources and not by their actual displacement. Second, we showed that the correlation of speckle intensities can be improved through simple design choices in the imaging setup; in particular, we found that the ME range is maximized when using a lens focused at the illuminator plane. Third, we showed that it is important to closely replicate near-field imaging conditions by ensuring that illuminators are accurately placed exactly behind the scattering sample, as doing otherwise can artificially increase the ME range. Fourth, we found that, in the near-field setting, the angular displacement for which significant correlations exist can correspond to actual displacements smaller than the illumination wavelength. As a consequence, ME approaches are only applicable to the near-field setting when considering scattering samples of modest thickness, where mid-order scattering is dominant. This thickness range still corresponds to penetration depths considerably beyond those achievable by a standard microscope. Therefore, analyzing and developing new ME approaches for imaging through scattering in this range can benefit biomedical imaging applications.
Our study additionally highlighted an important property of speckle intensity patterns formed due to samples where mid-order Fig. 20. Reconstructing fluorescent beads. We demonstrate singleshot near-field imaging of fluorescent beads attached at the back of â [100]μm chicken breast slice and the reconstruction produced by the fullframe and our local auto-correlation approach. Speckle contrast is low due to the weak emission of fluorescent beads (see a close-up on imaging noise in the inset) and due to the fact that the emission is not fully monochromatic.
scattering is dominant when imaged with a focused lens. These patterns have a small support size, typically much smaller than the sensor size. We showed theoretically and experimentally that using this local support to create a matched filter when computing speckle correlations can boost the SNR of latent illuminator detection by orders of magnitude. We additionally developed an algorithm that takes advantage of this property by operating as a local version of classical full-frame auto-correlation techniques. Through experiments on real measurements captured using both far-field and near-field imaging configurations, we showed that our algorithm provides an order-of-magnitude improvement in terms of both the range and the density of recoverable illuminators.
Furthermore, our study shed light on two fundamental challenges associated with the near-field case. The first challenge is the fact that the ME holds for very small displacements. The second challenge relates to the fact that, even after exploiting the local support, only sparse latent patterns are recoverable. These challenges still leave ample room for applications in medical imaging settings where sparse targets are considered, for example STORM imaging of blinking fluorescent molecules, sparse nuclei, or other cell components.
Additional challenges can arise due to reduced speckle contrast and signal-to-noise ratio in measurements captured under real near-field fluorescent imaging conditions. To highlight these challenges, in Figure 20, we show captured speckled patterns generated by a sparse set of fluorescent beads placed at the back of a tissue sample, as well as the reconstructions produced by the full-frame and our local auto-correlation techniques. Details on the experimental setup are provided in Appendix A.6. The quality of the input speckle images is severely affected by two factors. First, as in Figure 14, the light emitted by the fluorescent beads is spectrally broadband, reducing speckle contrast. Second, the fluorescent emission is very weak, leading to noisy images. Despite these challenges, we observe that our local auto-correlation algorithm significantly improves reconstruction quality compared to the full-frame baseline. Therefore, our results showcase both the strong potential of local auto-correlation techniques compared to full-frame variants, and the need for further research towards ME techniques that are fully robust to real experimental conditions in applications such as fluorescence microscopy.
Last but not least, by drawing attention to the local support characteristics of near-field speckle images, our results open the door for future research on using different image processing approaches for imaging through scattering, such as local deconvolution [Wu et al. 2020] and sharpening operations. A particularly promising direction is adapting the large array of mature blind deconvolution techniques to the imaging-through-scattering setting by developing appropriate prior models for the scatter-free images and (spatially varying) speckle blur kernels. Our study on the statistics of speckle patterns, and their dependence on the scattering layer geometry and material properties, can help inform the development of such priors.