Single scattering modeling of speckle correlation

—Coherent images of scattering materials, such as biological tissue, typically exhibit high-frequency intensity ﬂuctuations known as speckle. These seemingly noise-like speckle patterns have strong statistical correlation properties that have been successfully utilized by computational imaging systems in different application areas. Unfortunately, these properties are not well-understood, in part due to the difﬁculty of simulating physically-accurate speckle patterns. In this work, we propose a new model for speckle statistics based on a single scattering approximation, that is, the assumption that all light contributing to speckle correlation has scattered only once. Even though single-scattering models have been used in computer vision and graphics to approximate intensity images due to scattering, such models usually hold only for very optically thin materials, where light indeed does not scatter more than once. In contrast, we show that the single-scattering model for speckle correlation remains accurate for much thicker materials. We evaluate the accuracy of the single-scattering correlation model through exhaustive comparisons against an exact speckle correlation simulator. We additionally demonstrate the model’s accuracy through comparisons with real lab measurements. We show, that for many practical application settings, predictions from the single-scattering model are more accurate than those from other approximate models popular in optics, such as the diffusion and Fokker-Planck models. We show how to use the single-scattering model to derive closed-form expressions for speckle correlation, and how these expressions can facilitate the study of statistical speckle properties. In particular, we demonstrate that these expressions provide simple explanations for previously reported speckle properties, and lead to the discovery of new ones. Finally, we discuss potential applications for future computational imaging systems.


INTRODUCTION
When a camera images a scattering volume illuminated by coherent light, captured images are characterized by pseudo-random patterns called speckle. Despite their noiselike appearance, speckle patterns have strong statistical properties that provide rich information about the scattering material that is being imaged. For example, speckle images are approximately shift-invariant with respect to small perturbations of imaging parameters, a property known as the memory effect (ME) (Fig. 1) [1], [2], [3], [4]. These statistical properties have been studied extensively in optics [5], [6], [7], [8], and form the basis for many imaging techniques in application areas where scattering is important, such as medical imaging and remote sensing. Examples include seeing through a scattering layer and behind corners [9], [10], [11], as well as adaptive-optics focusing of light through highly-aberrated materials [12], [13], [14].
The large number of applications have motivated the creation of simulation tools and mathematical models for analyzing statistical properties of speckle statistics, without the need for painstaking lab measurements [15]. In particular, simulation tools include numerical wave-equation solvers [16], [17], [18], as well as Monte Carlo rendering algorithms [19], [20]. These simulation tools produce physically-accurate results, but they do not lend themselves to mathematical analysis, due to the complexity of the underlying equations that are being simulated. A different approach, that aims to facilitate analysis, is to replace the wave equations and multiple scattering process with a simpler mathematical model: such a model is only approximately correct, holding true only asymptotically under certain conditions, but is easier to analyze as it can be used to derive closed-form expressions of speckle statistics. This approach is well-established not only for speckle, but also when modeling incoherent intensity due to scattering materials. In the intensity case, the radiative transfer model of subsurface scattering [21], [22] is commonly replaced by approximations based on one of two extreme assumptions: The first assumption, termed single scattering, is that the scattering material is very optically thin, and thus most light paths do not scatter more than once [23], [24]. On the other extreme, the second assumption, termed diffusion, is that the scattering material is very optically thick, so that light paths scatter a large number of times, producing angularly-invariant radiance distributions [25], [26]. Models based on single scattering and diffusion have been used extensively for both rendering and material acquisition applications [23], [24], [25], [26], [27], [28].
Motivated by the challenges associated with simulating speckle statistics, there has been extensive attempts to go beyond models for smooth intensities, and develop approximate analytical models for speckle statistics [2], [3], [13], [14], [29], [30], [31], [32], [33]. Despite their approximate nature, these models have made it possible to better understand qualitative statistical properties of speckle, and have become valuable tools for developing speckle-based imaging algorithms [13]. However, most of these models are based on diffusion assumptions.
Our goal in this paper is to showcase the utility of an analytical model based on single scattering for explaining and analyzing experimentally-observed statistical properties of speckle patterns due to scattering. Compared to the incoherent intensity case, where single scattering models are only valid for very optically-thin materials, we show that the single scattering model for speckle statistics remain applicable for much thicker materials, where light paths scatter many times. As we explain, this is due to the fact that light paths that scatter multiple times are decorrelated and do not contribute to speckle correlation statistics. Additionally, we use simulations and real measurements to demonstrate that, for many practical situations where speckle-based imaging algorithms are applicable, predictions made by the single scattering model are more accurate than corresponding predictions by alternative models based on diffusion [3], [29] and Fokker-Planck assumptions [13], [14].
We use the single scattering model to derive closedform expressions for speckle statistics, and show that these expressions can be used to explain previously-reported statistical speckle properties, and explore previouslyundiscovered ones. We additionally discuss some potential applications of our model. We hope that, by facilitating a better understanding of speckle statistics, the single scattering model we introduce will lead to the development of better speckle-based imaging techniques.

Speckle correlation applications
Speckle correlations are important for many applications. For example, Katz et al. [10] showed that it is possible to recover a scattering-free image of a sparse set of mutuallyincoherent illuminators, observed through a thick scattering layer. They achieve this by computing the auto-correlation of the observed speckle image, and performing phase retrieval [34]. Analyzing expected speckle correlations can help understand under what imaging and material conditions this and similar imaging-through-scattering algorithms apply, and expand their applicability [35].
Another application of speckle correlations (Sec. 5.3) is using adaptive optics [36] to focus incident coherent wavefronts at points deep inside tissue samples. Determining the incident wavefronts is challenging, and typically involves using a guide star [12]. Computing expected speckle correlations makes it possible to adapt incident wavefronts for one point to focus at a neighboring point, without the need for a new guide star. It can also help understand over what ranges this refocusing is possible [13], [14], [20].
A third important application (Sec. 6) of speckle correlations is in scattering material acquisition [37]. Measuring speckle correlations, and inverting the relationship between material parameters and these correlations, can make it possible to infer these parameters for unknown materials.

BACKGROUND ON SPECKLE STATISTICS
We begin by reviewing speckle statistics and their path integral model [19]. We will use this model in Sec. 3 to present the single scattering approximation. Notation and setting: We use bold letters for vectors (e.g., points o, i, v), with a circumflex for unit vectors (e.g., directionsω,î,v). We consider scattering volumes V ∈ R 3 that satisfy the assumptions of radiative transfer for isotropic [38] media: each volume comprises a set of scatterers at statistically independent locations. Scatterers are assumed to be small enough relative to the wavelength of light to be considered infinitesimal points. They are also assumed to scatter incident light waves in a rotationallyinvariant way. We model speckle fields arising in such volumes due to monochromatic, fully-coherent and unpolarized incident illumination. These fields are a function of the volumes' bulk properties, which we describe next. Bulk material properties: We use a statistical description of the optical properties of scattering volumes. The scattering and absorption coefficients σ s and σ a model, respectively, the portion of energy that is scattered and absorbed upon interaction with scatterers. Their sum is the extinction coefficient σ t ≡ σ a + σ s , and its inverse is the mean free path, MFP ≡ 1/σ t , which is the average distance light travels between two scattering events. We often express the geometric dimensions of a volume V relative to its MFP . For example, a volume has optical depth OD = 2 if its thickness equals 2 · MFP , meaning that light in the volume undergoes on average two scattering events. The phase function ρ(î ·v) describes the amount of light from directionî scattered towards directionv. It is commonly characterized by an anisotropy parameter −1 ≤ g ≤ 1, equal to the average cos θ: g = 0 and g = 1 correspond to directionallyuniform and fully-forward scattering, respectively. Tissue is characterized by very forward scattering (g > 0.9) [39]. We assume that scattering volumes are homogeneous, meaning that the bulk parameters are spatially constant. Geometry: We focus on the geometry shown in Fig. 2. Without loss of generality, we assume that the optical axis is parallel to the z axis. We assume that the scattering sample is infinitely-wide along the x, y axes, has a thickness L along the z axis, and is located in the range z ∈ [0, L].
We consider two different configurations of this geometry: In a far-field configuration ( Fig. 2(a)), the sample is illuminated and imaged using directional light sources and sensorsî,v. Directional sensing can be realized by placing a 2D sensor at the focal plane of a lens as illustrated in Fig. 7. In a near-field configuration (Fig. 2(c)), the illumination is a point source i on plane z i close to or inside the sample (e.g., fluorescent particles inside tissue), or an illumination beam focused at the sample (e.g., confocal microscopy scanning). The sample is imaged by a camera focused at a plane z v inside or close to the sample. We follow Alterman et al. [35] and mostly consider a camera focused at the illuminator (a) In a far-field transmissive geometry, the sample is illuminated and viewed by directional beams from two different directions. (b) In a far-field reflective geometry, the illumination and viewing directions are at the same side of the sample. (c) In a near-field transmissive geometry, point illuminators are at the edge of the sample, and a camera viewing the sample from the opposite edge is focused at the illuminator plane. (d) The displacements ∆, τ between the two illuminators and between the illuminator and viewing point. We denote by∆,τ the angles that illuminators at these displacements form with a point at plane zo. (e) A visualization of the speckle pattern on a directional sensor (equivalently, a standard sensor at the Fourier plane of an imaging lens), where each point corresponds to the field observed from a different viewing directionv. The sample is illuminated simultaneously by two directional sourcesî 1 ,î 2 . plane z v = z i , so that in the absence of the sample, it would see a sharp image of the illuminators. We denote by v the 3D point where the camera is focused.
We denote by i x,y , v x,y , andî x,y ,v x,y the x, y coordinates of the illumination and viewing points and directions. The sample is illuminated and viewed from opposite sides in a transmissive geometry ( Fig. 2(a)), and from the same side in a reflective geometry ( Fig. 2(b)).
For this section and Sec. 3-4 we first consider the farfield case. We then adapt our results to the near-field case in Sec. 4.1.

Modeling far-field speckle covariance
We denote by uî(v) the complex field generated when light from sourceî propagates through the scattering sample, and is observed from directionv. Such a field produces a highfluctuation noise-like speckle intensity pattern (Fig. 1). We then consider a sample illuminated by two monochromatic and mutually-coherent sources atî 1 ,î 2 , and observed from two directionsv 1 ,v 2 . We define the speckle covariance as where * is complex conjugation, and expectation is taken with respect to multiple realizations of random media with the same statistical properties (e.g., tissue layers of the same type and thickness). Whenî 1 =î 2 =î,v 1 =v 2 =v,, the covariance Cî typically considered in computer graphics.

Covariance path integral
Bar et al. [19] showed that the speckle covariance of Eq. (1) can be expressed as a path integral involving only the bulk material properties. In particular, consider a path where B = 0, . . . , ∞. Light traveling along this path contributes energy equal to the path's complex throughput, describing amplitude and phase variations along the path. The amplitude term accounts for attenuation along the path, and is derived in App. A.1. The phase accumulated is proportional to the traveled path length, where accumulated phase between successive points is and the phase from the incoming direction to the first point or from the last point to the outgoing direction is where k = 2π/λ and λ is the wavelength of light. Bar et al. [19] derive their path integral expression by first considering the correlation of fields that travel along all possible pairs of paths ( x 1 , x 2 ), x 1 fromî 1 tov 1 , and x 2 from i 2 tov 2 . The contribution from a pair of paths has phase where ∝ in Eq. (6) denotes equivalence up to a scale factor. The exact path contribution is provided in App. A.1 and includes an additional amplitude term. A first path integral expression for covariance can be obtained by integrating over all such pairwise path contributions ( Fig. 3(a)). However, Bar et al. [19] made the following observation: Given the large variation between the lengths of different paths, summing path contributions over all path pairs involves summing complex numbers with phases that are essentially uniformly random. As the expected average of complex numbers with uniformly random phases is zero, paths with random phase do not contribute, in expectation, to the covariance. Based on this observation, Bar et al. [19] showed that the path integral can be simplified to use only pairs of paths that coincide everywhere, except for their connections toî 1 ,v 1 ,î 2 ,v 2 , visualized in Fig. 3(b). In particular, consider the space P of sub-paths x s = o 1 → · · · → o B , B ≥ 1, where each vertex o b ∈ V. These vertices correspond to the shared part of two full paths x 1 =î 1 →o 1 →. . .→o B →v 1 ,   [19] show that most such pairs do not contribute to the correlation, and simplify the path integral by restricting it to pairs of paths sharing all their nodes except for the start and end ones. (c) We show that even this restricted path space includes many decorrelated path pairs, and speckle covariance can be mostly attributed to paths of length 1.
sub-path toî 1 ,v 1 andî 2 ,v 2 . Then, the speckle covariance of Eq. (1) can be expressed as the covariance path integral: where the far-field path contribution function cî 1 ,î 2 v 1 ,v 2 equals the correlation of the fields that travel along x 1 , x 2 . The phase of cî 1 ,î 2 v 1 ,v 2 reduces to the phases of only the four segments connecting the shared path to the start and end nodes: (8) This is due to the fact that all other segments of the path are shared between x 1 , x 2 and have equivalent lengths.
Bar et. al. [19] have comprehensively evaluated this model to establish its accuracy. We replicate in Fig. 4 one of their comparisons between covariance matrices obtained with a wave equation solver [16] (accounting for all pairwise path contributions as in Fig. 3(a)), versus a Monte Carlo estimate of (7) (restricting integration to only joint paths as in Fig. 3(b)). The two matrices are identical, confirming that disjoint path pairs can be removed from the path integral. We performed this simulation in 2D, as the wave equation solver [16] is restricted to this setting. We fixed two illuminationsî 1 ,î 2 displaced by an angle of 4 • . In 2D, viewing directions form a 1D space spanned, e.g., by the angle they form with the main axis. Thus, Fig. 4

THE MEMORY EFFECT PROPERTY
With the simplified covariance path integral of Eq. (7) at hand, we can now ask the question: what pairs of illumination and viewing conditionsî 1 ,î 2 ,v 1 ,v 2 result in large covariance values Cî The memory effect (ME) property of speckle patterns refers to the existence of pairs of illumination and viewing conditions producing highly-correlated speckle fields and intensity images, and a large part of the Waves solver Full M.C. Single Scat.  [19]. (c) Covariance from our single scattering approximation. All approaches are in good agreement, even though the simulated sample has OD = 2. Non-zero correlation is present only along a diagonal of viewing directions satisfying the ME conditions of Eq. (9). literature on speckle-based imaging techniques is dedicated to discovering and characterizing such pairs.
The covariance path integral of Eq. (7) can guide us in searching for conditions under which covariance attains a large value. In particular, we note that the pairwise path contribution function cî (8)) can be complex: This is because the phase of the four source and sensor connection segments ξ(î 1 →o 1 ), ξ(î 2 →o 1 ), ξ(o B →v 1 ), ξ(o B →v 2 ) can differ significantly, as the path vertices o 1 and o B vary in the volume. As a result, for most selections of illumination and viewing pairsî 1 ,î 2 ,v 1 ,v 2 , the phases of the path contributions cancel out and the covariance integral of Eq. (7) becomes zero in expectation. Therefore, configurations where the memory effect applies will correspond to cases where many cî 1 ,î 2 v 1 ,v 2 terms have similar phases. The simplest such configuration is when we selectî 1 =î 2 andv 1 =v 2 , which corresponds to computing intensity.
Outside of the trivial intensity case, the best-known pairs of illumination and viewing conditions producing high covariance correspond to the translational memory effect: these are pairs of illumination and viewing directions separated by the same sufficiently small displacement |∆|, This type of correlation has been studied extensively in the literature [1], [2], [3], [4]. To demonstrate it, we fixed in Fig. 4 two illuminationsî 1 ,î 2 in 2D, and evaluated the covariance matrix Cî 1 ,î 2 v 1 ,v 2 for all pairwise selections of the viewing directionsv 1 ,v 2 . We see that non-zero correlations are present only along a diagonal of viewing directions, which are exactly the direction pairs satisfying Eq. (9), given the fixed illuminators. Based on the above discussion, we expect that these illumination and viewing conditions result in a large set of paths for which the phases in Eq. (8) are similar. In Sec. 4, we use a single scattering approximation to the covariance path integral to investigate these conditions and the memory effect. Before we do so, we review previous diffusion theory models for explaining the memory effect [3], [29]. Even though this is somewhat different than how these models were previously derived, we argue that diffusion theory also attempts to study paths with similar phase contributions.

Diffusion theory models
We denote byî 1 z ,î 2 z ,v 1 z ,v 2 z the z coordinate of the illumination and viewing directions. In diffusion theory, we assume that these directions are close enough to the optical axis for their z coordinate to approximately equal 1. Then, We now consider a sub-path starting at o 1 and ending at o B . The phase contribution of Eq. (8) becomes where in Eq. (11) we substitute the definition of phase between a directional source or sensor and a path vertex (Eq. (5)), and in Eq. (12) we use the assumption that the z coordinate is approximately zero (Eq. (10)), so we only need to consider the xy coordinates of the start and end nodes o 1 , o B , and the xy component of the difference vector nDl (Eq. (9)). As we argued previously, if the complex numbers in Eq. (12) can have varying phases, the sum of complex contributions from multiple random paths becomes zero in expectation. The summation is non zero only if the phase difference between different paths is negligible; from Eq. (12), this corresponds to k|∆||o 1xy − o Bxy | < π, or equivalently This bound explains why ME correlation exists only for a small range of illumination displacements∆. Diffusion theory models [3], [29] can be used to approximately evaluate the expectation E[|o 1xy −o Bxy |] of Eq. (11) as a function of the scattering material depth and bulk optical properties, and thus derive analytic bounds on the memory effect. Intuitively, under diffusive conditions, the average distance between an entrance point and an exit point on the scattering sample scales as the optical depth and the sample thickness L increase. This difference is smaller when the phase function of the material is primarily forward scattering.
In Fig. 5 we plot correlation as a function of∆ as predicted by diffusion theory models, and compare it with predictions produced using the Monte Carlo rendering algorithm of Bar et al. [19]. For these comparisons, we use material parameters that have been reported in the literature for tissue (phase function anisotropy g = 0.97, mean free path 50µm). We consider material thicknesses L = 50, 150, 5000µm corresponding to optical depths OD = 1, 3, 100 respectively. From the comparisons, we observe that diffusion theory predictions are accurate only for materials of very large optical depth. As previously reported in the literature [15], for materials of smaller optical depths, the decay predicted by diffusion theory is significantly more conservative than what is measured experimentally.  We compare the diffusion theory [3], [29] model for the ME with the exact prediction of the Monte-Carlo simulator [19], and our single scattering approximation. The different subfigs visualize volumes of three different thicknesses corresponding to different optical densities. We use scattering parameters representing tissue, in particular a narrow forward scattering phase function with g = 0.97. The single scattering approximation is invalid at low angles∆ where intensity dominates, but for larger ones it explains all correlation even for material thickness of OD = 3 where light paths scatter more than once. This suggests that multiple scattering paths decorrelate quickly and do not contribute to the correlation. In contrast, due to the narrow phase function, the diffusion model becomes accurate only for very high optical depths as visualized in (c). However, for high optical depth, the memory effect is only valid for tiny ranges (note the different ∆ ranges at which the three subfigures are plotted), and most useful ME applications are in the low and medium OD regimes.
Additionally, we note that for very thick materials, where the diffusion approximation is accurate, ME correlation exists only for very small, almost negligible, displacements |∆| (see the different∆ range of the subfigures in Fig. 5). By contrast, speckle-based imaging algorithms are more practical in situations where ME correlation exists for significantly larger displacements |∆|; thus these algorithms are most applicable to materials of modest optical depth [35]. Therefore, we note that the diffusion approximation becomes inaccurate under the conditions most suitable for speckle-based imaging algorithms. This motivates us to develop in the next section an alternative approximate model, based on a single scattering assumption, that can be used to study the memory effect under these conditions.

FAR-FIELD SINGLE-SCATTERING COVARIANCE
We now introduce our single scattering model of speckle covariance. We first present this model for the far-field geometry, and discuss the near-field geometry in Sec. 4.1. Our starting point is the simplified covariance path integral of Eq. (7). As we argued in Sec. 3, strong correlations will exist only under illumination and imaging conditions for which there is a large subset of path pairs where the contributions (Eq. (8)) of the four start and end nodes have similar phase.
Our argument is that, for many practical imaging configurations resulting in strong correlations, the subset of paths with similar phases mostly consists of single-scattering paths.
As we discuss below, approximating the solution using single scattering paths is equivalent to what is known in the optics literature as the first Born approximation [40]. Effectively, this allows us to further simplify the covariance path integral of Eq. (7), by only considering subpaths of length B = 1, resulting in pairs of paths of the form Fig. 3(c).
To provide some qualitative justification for this single scattering model, we first compare it with the exact covariance in Fig. 4 showing good agreement. While Fig. 4 considered a single illuminators displacement∆, for better understanding in Fig. 5 we plot correlation as a function of the displacement∆ between illuminators, this time averaging over all viewing directions which satisfy the ME conditions in Eq. (9) (i.e. viewing directions located on the diagonal of the matrices in Fig. 4). As |∆| increases, the single scattering model is remarkably accurate (in Fig. 5(a,b) the blue and red plots match when |∆| is above 5 • ), despite the fact that it considers only single-scattered light. The single scattering model remains accurate even when the optical density is greater than one, when most paths scatter multiple times (e.g. OD = 3 in Fig. 5(b)). Contrasting the predictions by the single scattering and diffusion models, we note that the diffusion model is accurate only at very large optical depths where the ME is negligible, whereas the single scattering model is accurate at intermediate optical depths where significant correlation exists and the ME is more applicable [35]. The phase-matching of single-scattering paths: We can now use our single-scattering model to study the memory effect. To this end, we note that the single scattering model makes it easier to characterize illumination and viewing configurations that will result in paths whose contributions (Eq. (8)) have a fixed phase, and thus that will produce strong correlations as argued in Sec. 3.
We first consider pairs satisfying the memory effect condition in Eq. (9). As we aim to explain correlation behavior in angles wider than the ones addressed by diffusion theory, we avoid the assumption that the z coordinate of the difference vector is zero, as was applied in Eq. (10). Rather, we use the notation: For single scattering paths where o = o 1 = o B , and for configurations satisfying the ME conditionî 2 x,y −î 1 x,y = v 2 x,y −v 1 x,y , the x, y coordinates of the exponent in Eq. (11) reduce to zero. As a result, the phase of the contribution from a path in Eq. (11) reduces to where o z is the scalar of the z coordinate of o. This phase is independent of the xy position of the scattering point o.
As a result, when the scattering layer is very thin and all scattering points o lie on a plane with the same z coordinate, the phase of all path contributions is constant, producing a large covariance. When the scattering layer is thick, o z can vary significantly, and thus the term ω z o z is not constant; consequently, the path contribution cî can have varying phase and the covariance path integral decays to zero. We observe, then, that for thicker material layers the correlation is low, unless we consider an illumination and viewing configuration for which ω z = 0. Apart from the simple case considered by diffusion theory (Eq. (10)) where all directions are close to the optical axis, below we identify illumination and viewpoint selections leading to ω z = 0, and show that such configurations indeed lead to higher correlations.
Closed-form derivation: An advantage of the single scattering approximation is that it makes it possible to express the covariance path integral of Eq. (7) in closed form. We derive this closed-form expression below, and we will be using it in subsequent sections to analyze the ME.
To this end, we first introduce a more convenient parameterization for illumination and viewing configurations satisfying the ME: We use the notationτ to describe the 2D displacement between the illumination and viewing vectorŝ τ ≡î 1 x,y −v 1 x,y , as visualized in Fig. 2(e). For configurations satisfying the ME conditions of Eq. (9) we have that With this parameterization we show in App. A.2 that thus the phase of the path contribution in Eq. (15) can be expressed as cî In App. A.2 we derive a closed-form expression for the covariance resulting from single scattering paths, this time taking into account the attenuation component of the path contribution, instead of only considering its phase in Eq. (18). The following claim shows the derived expression for the transmissive far-field geometry, and we show the corresponding expression for the reflective case in App. A.3. In this claim, we denote by ρ(τ ) the phase function evaluated at the angle between the illumination and viewing direction, which we approximate asτ (see definition in Eq. (16)).

Claim 1. Consider illumination and viewing pairs satisfying
the memory effect conditions, parameterized by∆,τ as in Eqs. (9) and (16). The single-scattering covariance between fields resulting from scattering points at depth z is C(τ ,∆, z) = ρ(τ )σ s e −σtL e ikz(∆·τ ) (19) and the total covariance where the target is a scattering layer of thickness L located at the range [0, L], and z o = L/2.
We will study the implications of this result in detail in Sec. 5. At a high level, Eq. (19) tells us that the correlation has an amplitude which is proportional to the phase function at angleτ , and decays exponentially with σ t L, which is essentially the optical depth of the scattering layer. The correlation of complex speckle fields is a complex number whose phase varies as exp(ikz(∆ ·τ )). The exact phase of this sinusoid depends on the z plane on which the scatterers lie. In Eq. (20) we integrate the covariance over different z planes. Due to the varying phase the integrated correlation reduces. The amplitude reduction is shown to be proportional to the sinc of the angle (∆ ·τ ), and will thus be higher when the displacements∆,τ are orthogonal. Accuracy of the single scattering approximation: We compare the single scattering approximation against the full Monte-Carlo simulator of [19]. This simulator was exhaus-tively compared against exact speckle covariance evaluated using exact wave solvers [16], [17], [18], establishing its physical accuracy. In Fig. 5, we have shown a first comparison between the diffusion, single-scattering, and full Monte Carlo rendering approaches, by plotting covariance as a function of∆ alone. For this comparison, we fixed i 1 x,y = −0.5∆,î 2 x,y = 0.5∆ and computed In Fig. 6 we visualize all pairs satisfying the ME Eq. (9) for a given selection of illumination directions. We note that while in Fig. 4, which demonstrates a flatland simulation, viewing pairs satisfying the ME conditions lie on a 1D diagonal, for a 3D scene there is a 2D space of such pairs, which we parameterize usingτ (Eq. (16)). We show such 2D correlation maps for a few selections of illuminators displacements∆.
We used optical parameters typical of tissue (phase function anisotropy g = 0.97, mean free path 50µm) and various optical depths. We compare the covariance computed using [19], with the one computed from the single-scattering approximation in Eq. (20). For∆ = 0, we have thatî 1 =î 2 , v 1 =v 2 , and the correlation reduces to the intensity of the scattered light. We observe that the contribution to the intensity from paths that scattered multiple times is significant. By contrast, for large displacements∆ (in Fig. 6, displacements larger than 5 • ), for which correlation and intensity are not the same, the correlation agrees well with the single scattering approximation. This remains true even for OD = 5, that is, for a material thick enough that most light paths scatter more than once. We observe the same phenomenon in both transmissive and reflective geometries. This suggests that the multiple-scattering paths are decorrelated and do not contribute to correlation. Therefore, the single scattering model remains an accurate approximation for describing speckle covariance even for materials of large optical thickness; this is in contrast to its use for describing intensity, where it is an accurate approximation only for very optically thin materials. In addition to the narrow forward scattering phase function evaluated in Fig. 6, in Appendix Fig. 10 we repeat the experiments with a wide isotropic phase function (g = 0). Similar behavior is present, but due to the wider scattering angles the single scattering approximation decays to zero at smaller optical depths.
The availability of a closed-form expression for single scattering covariance can be important for two tasks. First, as we discuss in Sec. 5, we can use this expression to provide simple explanations for various properties of the memory effect, as well as predict new ones. Second, we can use this expression to invert correlation measurements and acquire bulk parameters of unknown materials, in a way that is significantly simpler than previous inverse scattering approaches [37]. We do not pursue this second direction in this paper, but we briefly discuss it in Sec. 6. Single scattering and the first Born approximation: Approximating the solution to the wave equation using single scattering paths is equivalent to the popular first Born approximation [40] from the optics literature. There is, however, an important difference between the two models: The Born approximation models the single scattering component of a speckle field arising from a specific instantiation of the scattering medium (a single, fixed set of scatterer positions). By contrast, our model computes the expected singlescattering correlation arising from a scattering medium with some bulk scattering properties, describing the statistical distribution of scatterers. To obtain our result using the Born approximation, one would need to average many different single-scattering speckle fields (as in Eq. (1)), arising from scatterer instantiations randomly drawn from the bulk scattering properties. Effectively, our model directly accounts for this ensemble averaging process in a mathematically tractable and computationally efficient way.

Near-field single scattering correlation
Up until now, we have applied the single-scattering approximation for far-field configurations, consisting of directional collimated illuminators and sensors. In this section we provide a similar single scattering approximation for speckle covariance in near-field configurations, that is located near or inside the scattering slab or focused at it. To simplify the resulting expression, we assume that the source and sensor are located or focused at the back plane of the sample, so that z v = z i = 0. We denote by x,y the restriction of the 3D illumination and viewing points to this plane. As before we consider ME pairs satisfying x,y . Unlike Eqs. (9) and (16), ∆, τ denote displacements between near-field points rather than directions, so we denote them without circumflexes. Then, in App. A.4 we show that the covariance resulting from scatterers on a single z plane is The total covariance from integrating over all z positions is Unlike the far-field case, in Eq. (22) the phase of the correlation at each z plane depends on the scalar product ( ∆ /z· τ /z), which varies with z. Thus the integration in Eq. (23) cannot be evaluated in closed-form, but we can still evaluate it numerically in an efficient way. We note that the formula does not depend on absolute shifts τ , ∆ of the illumination and viewing points, but rather, on the normalized displacements τ /z, ∆/z. These normalized displacements are equal to a first-order approximation to the angle that illuminators or sensors at distances τ , ∆ form with a scattering point at depth z, as visualized in Fig. 2(d). As a result, we will be denoting the angular displacements corresponding to near-field displacements using circumflexes:τ with z o = L/2 denoting the mean depth of the scattering volume. Fig. 6 simulates near-field covariances in transmissive geometry. As in the far-field case, we observe that for large displacements∆, the single scattering component explains most of the correlation. To better match with the far-field case the visualization in Fig. 6 presents the near field correlation as a function of the angular displacementsτ ,∆ rather than τ , ∆.

MEMORY EFFECT PROPERTIES
In this section we use the single scattering approximation to explain some properties of the memory effect that have been previously reported in the literature. Often this explanation is simpler and more accurate than previous explanations based on diffusion or other approximate models. We also point out some new observations that were not previously reported.

Local support
Revisiting the covariances in Fig. 6, we observe that the correlation can vary significantly as a function of the displacementτ . For modest OD, the correlation has a local support, being high only for lowτ displacements. This is expected given that the simulation in Fig. 6 considered scattering parameters typical of tissue, with a very narrow forward scattering phase function. The result of using such a phase function is that incoming illumination at direction i will, with very high probability, scatter towards viewing directions closer to the incoming directionv ≈î. In [35] it was shown that this local support property of speckle patterns can lead to speckle-based imaging through scattering algorithms [10] with significantly improved performance in terms of the types of hidden illuminator patterns that can be recovered. Therefore, being able to reproduce this property is a significant feature of the single scattering model.
We note that, even though this local support property is predicted by the single scattering model, this is not the case for previous analytic models for the ME based on diffusion and the Fokker-Planck assumptions [3], [13], [14], [29]. Those models describe correlation as a function of only the illuminator displacement∆, and do not take into account the angleτ between illumination and viewing directions. Effectively, these models assume that correlation remains constant for all viewing directions, regardless of their relationship to the illumination directions. The absence of theτ parameter in previous models provides further evidence that, for many settings of practical interest, single scattering is a better approximate model for ME correlations than diffusion or Fokker-Planck models.
As we mentioned in the previous section, predictions made using the diffusion approximation are consistent with the full Monte Carlo simulations only for very thick samples (OD = 30). At such high optical depths, the exact correlation indeed becomes invariant toτ (eachτ square in the OD = 30 column in Fig. 6 is uniform). At this high optical depth, we observe strong values for∆ = 0, which simulates intensity. As the displacement∆ increases even by a small amount, correlation decays sharply (in Fig. 6, allτ squares for angles above 0.05 • are zero). We also note that the single scattering covariance is inaccurate for large optical depths as it always predicts zero correlation, as seen in the OD = 30 column of Fig. 6.

Anisotropic support of memory effect
In Fig. 6 illuminators are displaced on the vertical axis. When the material is of medium thickness (e.g., OD ≤ 5) we observe another interesting property of speckle covariance: in both the near-field and far-field cases, and in both transmissive and reflective geometries, the area with high correlation has an anisotropic support, with the long axis of this area being orthogonal to the displacement between the illuminators. In Fig. 6, as illuminators are displaced on the vertical axis, this leads to longer horizontal lines in the C(τ ,∆) correlation images. The anisotropy increases as the distance∆ between the illuminators increases.
We can use the single scattering model to explain this property. In particular, for every z plane, the covariance in Eq. (19) is a sinusoid with complex phase e ikz(∆·τ ) ; (25) integrating these complex values over different z planes reduces the covariance. However, the phase variation is different for differentτ selections, explaining the anisotropic shape. If the viewpoint displacementτ is orthogonal to the illumination displacement∆, i.e., (∆ ·τ ) = 0, there is zero phase for all depth planes z, resulting in a large integrated covariance value. Ifτ is parallel to the direction of∆, then (∆ ·τ ) is large and the significant phase variations across z planes result in a reduced integrated covariance value. For the far-field case, this integral can be expressed in closedform as in Eq. (20), leading to a sinusoidal phase change as well as a sinc amplitude along the direction whereτ is parallel to∆. In Fig. 6, which displays the real part of the correlation, the sinusoidal phase is visible as multiple cycles along the vertical direction (in this example, the illuminator displacement is vertical, thereforeτ is parallel to∆ along the vertical axis).
As far as we are aware, the anisotropic support property of the memory effect has not previously been reported in the literature. In particular, as anisotropy is a function of the displacement between the illumination and viewing direction τ , previous memory effect models [3], [13], [14], [29] that do not consider τ fail to reproduce it. This new property can be used for more accurate modeling of correlation support in algorithms that take advantage of the local support of speckle [35], or to further improve adaptive optics focusing [13], [14].
In Fig. 7 we have confirmed this property through real lab measurements, using transmissive far-field images of a sample consisting of a chicken breast slice of thickness 500µm. We used a far-field speckle acquisition setup shown in the figure, to capture multiple speckle intensity images Iî(v) under different illuminations. Some of the captured speckle images are visualized in Fig. 1. We compute empirical correlations averaging speckle intensity product values over a small local window. We did this for illumination shifts in both the horizontal and vertical directions. Fig. 7 shows the results, where we observe that correlation decays as the illuminator displacement increases, with the decay being faster in the direction of the illuminator displacement than in the orthogonal direction. Note that Fig. 7 displays correlations of intensity images rather than fields. Classical statistics states that, for zero-mean fields, this is equivalent to C(τ ,∆) 2 . Differences between the measurements in- Fig. 7 and the Monte Carlo and single scattering simulations in Fig. 6 are due to mismatch between the bulk material parameters of the real and virtual samples.

Tilt-shift memory effect
Osnabrugge et al. [13] have proposed a model for the tiltshift memory effect of near-field speckle patterns, which is derived from another popular optics approximation-the Fokker-Planck model. In this section, we first review the tilt-shift memory effect and explain its importance for an application of the memory effect, adaptive-optics focusing. We then show how our single scattering model can provide more accurate predictions about the tilt-shift memory effect than the original Fokker-Planck model of [13].
The tilt-shift memory effect can be used to make it easier to focus light through tissue into a spot [13], [14]. A beam of light passing through tissue is scattered and thus cannot focus into a single point. Adaptive-optics focusing uses a spatial light modulator, to create a coherent wavefront that is conjugate to the aberration due to the scattering inside the tissue. This wavefront is focused into a spot after passing through the tissue. The wavefront correction pattern applied by the spatial light modulator is specific to the tissue sample being imaged and the intended focus point. Finding the exact shape of this wavefront is challenging [12], [41], [42]. The tilt-shift memory effect simplifies this process, by making it possible to modify the wavefront required to focus at a specific point, so as to achieve refocus at a nearby location [14]. This adaptation usually involves small tilting and shifting of the wavefront correction, hence the name tilt-shift memory effect. We refer the reader to [13], [14], [20] for details.
The refocusing tasks described above in the context of adaptive optics scanning can be formulated mathematically in the context of speckle covariance as follows: Consider two complex scattered fields u i 1 (v), u i 2 (v) generated by the near-field illumination points i 1 x,y and i 2 x,y = i 1 x,y + ∆. Can we find a phase correction ψ(v) that will maximize the expected correlation As the correlation is complex, the ideal phase correction is just the conjugate of the expected correlation C i 1 ,i 2 v,v+∆ . The phase of this correlation can be evaluated using Monte Carlo rendering [19], [20]. Instead, Osnabrugge et al. [13] consider a sinusoidal phase correction of the form ψ(τ ) = (θ · τ ), where θ is a 2D vector and τ = i 1 − v 1 . These sinusoidal corrections are attractive in an optical setup, as they can be implemented as a simple tilt of the field by an amount equal to θ. To this end, they define the tilt-shift correlation as where Eq. (27) is the equivalent of Eq. (26) for the specific case of a sinusoidal phase correction, and where for the simplicity of subsequent analysis we express summation over sensor pixels using the parameterization v x,y = i 1 x,y + τ . Assuming the illumination and sensor points are located at the back edge of the medium, at z i = 0, we can adapt an analytic formula derived Osnabrugge et al. [13] to express the tilt-shift correlation: 1 where tr is the transport mean free path tr = MFP /1−g.
Their derivation is based on the Fokker-Planck model, which makes three simplifying assumptions: it uses a multislice layered representation of wave propagation [15]; it assumes forward-only propagation; and it uses a differential equation to integrate over multiple scattering planes. Using Eq. (28), Osnabrugge et al. can also predict that for a displacement ∆, the phase correction maximizing the correlation corresponds to a tilt at angle This optimal tilt angle is essentially the angle∆ defined in Eq. (24), that is, the angle a pair of sources at displacement ∆ form with a scattering point at depth z o = (2/3)L, as illustrated in Fig. 2(d). Revisiting the single scattering correlation formula in Eqs. (22) and (23), this is not surprising. For a thin scattering layer at a single z plane, the correlation is indeed a sinusoid with frequency ∆/z, thus it is reasonable to expect that, when integrating correlations from multiple 1. The formula of [13] assumed the illumination sources are located at the back edge of the medium while the sensor is focused at the front edge. In [35], it is shown that changing the sensor focal plane is equivalent to a simple shearing of the coordinate system. In this paper we pre-applied the shear and the result in Eq. (28) is an adaptation of the result in [13] using a sensor at the left edge. This coordinate choice simplifies the resulting formula.  [13]. (b) Empirical lab measurement of tilt-shift correlation by [13]. (c) Tilt-shift correlation by our single-scattering model, which better captures the double wedge shape of the measured correlation.
(d-f) To explain the double wedge shape of the correlation we note that if we compute the correlation from particles inside a thin z slice we get a narrow line, whose slope varies as a function of the plane depth z. The total tilt shift correlation in (b,c) is the average of all these slopes.
depth planes in the range z ∈ [0, L], we end up with a phase corresponding to one of the intermediate planes.
We note that in Eq. (24) we choose the plane z o with which we compute the angle as z o = L/2, rather than z o = (2/3)L as in the result of Osnabrugge et al. [13]. In App. A.5 we compare these two choices, showing that in practice the difference between them is minor, and explain the motivation for choosing z o = L/2.
In Fig. 8 we compare the tilt-shift covariance predicted by the model of Osnabrugge et al. [13], as expressed in Eq. (28), to the covariance predicted by the single scattering model. We also compare with an empirical lab measurement of this correlation provided by Osnabrugge et al. [13]. We used material parameters equivalent to those reported by Osnabrugge et al. [13], corresponding to g = 0.98, MFP = 296µm and L = 258µm. Both approximate models produce a correlation function with a dominant lobe that has the same orientation as the one in the measured data. However, our single scattering model matches the overall shape of the measured correlation more closely than the Fokker-Planck model of Osnabrugge et al. [13].
Note that, following the definition in Eq. (27), every entry of the correlation image in Fig. 8 is equal to the 2D summation over the τ coordinates of a correlation image as in Fig. 6, multiplied by a sinusoid at a different frequency We also note that, even though the analytic formula of Osnabrugge et al. [13] predicts that the shape of the correlation is a slanted Gaussian lobe, the measured tilt-shift correlation has the shape of a double wedge. The single scattering model not only predicts this double wedge, but can also help us understand its source. In Fig. 8(d-e) we evaluate a configuration where the scatterers are placed only at a thin layer out of the range [0, L]. As suggested by Eq. (22), the correlation from particles at depth z has a phase ∆/z, thus the correlation from a thin layer at depth z is maximized by a tilt angle θ = ∆/z. As a result, the contribution of particles at depth z to the tilt-shift correlation map visualized in Fig. 8(d-e) is a line at slope θ = ∆/z. The wedge in Fig. 8(c) is essentially an integral of all slopes in the range.
To summarize, the single scattering correlation model provides a simple explanation for the near-field tilt-shift correlation derived by Osnabrugge et al. [13]. Single scattering explains the lab measurements better than the Fokker-Planck model they derive. As mentioned earlier, the tiltshift correlation is very important for finding a local tiltshift adaptation of a wavefront shaping pattern. [20] showed that using accurate models for speckle covariance helps improve the refocusing capability. Our single scattering model provides a middle ground between improved accuracy and computational efficiency, and thus we hope it will help proliferate adaptive optics focusing techniques in future imaging systems.

BULK PARAMETER ESTIMATION
The single scattering approximation allows us to derive closed-form expressions relating bulk material parameters with speckle covariance values. These expressions can be used to invert measurements of speckle covariance, to infer from them the bulk parameters of unknown scattering materials. This can potentially lead to significant simplification of the computation required for accurate inverse scattering [37]. We briefly outline this idea below.
In material acquisition systems, one first collects image measurements of unknown scattering materials under different illumination and viewing conditions, then searches for material parameters that, when used for simulation, can explain the captured data. This approach is termed inverse rendering. In intensity-based inverse rendering, the forward simulation process requires using Monte Carlo volumetric rendering; consequently inverse rendering becomes a complex optimization problem [37]. To relax the computational burden, previous works have adopted approximate scattering models that can predict images using closedform expressions, and thus can be inverted efficiently. In particular, if the optical depth of the material is small, one can the single scattering model. As the single scattering model, when applied to intensity, holds only for very low optical depths, acquisition systems based on this model had to dilute samples to reduce their optical thickness [23].
As suggested in this paper, single scattering models applied to speckle covariance can hold for larger optical depths, even if most light paths scatter multiple times. Reviewing the closed-form expression for the single scattering correlation in Eq. (20), we note that we can essentially directly extract values of the phase function of the material at different angles from the empirical speckle correlation. We need to select illumination directionsî 1 ,î 2 where the displacement∆ =î 2 x,y −î 1 x,y is sufficiently large for the single scattering approximation to be accurate. Then, to measure the phase function at angleτ , we need to select viewing directions satisfyingv 1 x,y −î 1 x,y =v 2 x,y −î 2 x,y =τ . The phase function can be retrieved from an empirical averaging of the product Iî 1 (v 1 ) · Iî 2 (v 2 ).

DISCUSSION
We introduced a new model for speckle correlation, utilizing only the single-scattering component of scattered light. Remarkably, even in thick scattering samples where most light paths scatter multiple times, the correlation can be explained by paths that scatter only once, as longer paths decorrelate quickly. We evaluated our model using both synthetic and real data, and showed that for many practical scenarios where the memory effect is useful (e.g., low and medium optical densities), predictions from our model are more accurate than those made by models based on diffusion theory or the Fokker-Planck approximation. An advantage of our single scattering model is that it can be evaluated in closed-form, sidestepping the computationally intensive Monte Carlo process required for an exact evaluation. Moreover, its closed-form expressions can help analyze the properties of speckle correlations. We showed that our model can explain previously reported speckle correlation properties, and also reveals new, previously unexplored, ones. Another important application of our single scattering model is that it can allow estimating material properties from empirical speckle statistics, without the need for complex inverse rendering algorithms [37]. We hope that our model will facilitate the design of future speckle-based computational imaging systems.

APPENDIX A
Our goal here is to derive close-form expressions for the single-scattering components of the speckle covariance. Before addressing the single-scattering case we provide the complete formula for path throughput.

A.1 The full Monte-Carlo path throughput
In Sec. 2.2 we have derived the path throughput and the contribution from a pair of paths. However, we only expressed the phase of the paths. In practice the exact contribution of each path has an amplitude as well, which encodes exponential attenuation along the path. For the simplicity of the exposition we have neglected this attenuation in the main paper and we provide the complete formulas below.
Below we will denote by s(cos θ) the scattering amplitude function, which describes how a field interacts with a scatterer: if a scatterer is illuminated from directionî, the complex scattered field u at directionv is uîv = s(î ·v). The standard phase function used in computer graphics to describe scattering is defined as ρ(cos θ) ≡ |s(cos θ)| 2 .
Rather than considering the space of all path pairs, Bar et al. [19] showed that the path space can be simplified to use only pairs of paths that coincide everywhere, except for their connections toî 1 ,v 1 ,î 2 ,v 2 , visualized in Fig. 3 the direction of the b-th edge of the sub-path. These vertices correspond to the shared part of two full paths formed by connecting the sub-path toî 1 ,v 1 andî 2 ,v 2 . Then, the speckle covariance of Eq. (1) can be expressed as: where the far-field path contribution function cî 1 ,î 2 v 1 ,v 2 equals the correlation of the fields that travel along x 1 , x 2 . For B ≥ 2, this equals: and for B = 1: In the above, f ( x s ) is the real and positive standard radiometric throughput of x s , augmented by scattering coefficients at the first and last vertex, Finally, υ(·) is the complex volumetric throughput, defined as: where d(ω→o), d(o→ω) denote the distance a ray entering or leaving, respectively, o at directionω, travels inside the scattering volume V.
The covariance rendering algorithm of Bar et al. [19] uses a Monte Carlo path sampling approach to evaluate the speckle covariance integral of Eq. (31). This algorithm takes advantage of the presence of the radiometric throughput term in Eq. (32), and samples sub-paths x s using standard volumetric path tracing. Then, for each sampled sub-path, the endpoints o 1 , o B are connected to the far-field illuminationsî 1 ,î 2 and sensorsv 1 ,v 2 , to compute the complex volumetric throughput terms in Eqs. (32) and (33).

A.2 Far field transmissive geometry
We turn our attention to the derivation of closed-form expressions for the single scattering covariance. We start with the far field transmissive geometry case. We want to prove claim 1 stating that the correlation is non zero only for illumination and viewing directions satisfyinĝ For such directions, the contribution to the single scattering correlation from points on a plane at depth z is with∆ To see this, consider a scattering point o and refer to the definition of its single scattering contribution to the correlation as defined in Eq. (33).
To simplify this, note that when the anglesî 1 ,î 2 ,v 1 ,v 2 form with theẑ axis are small, and we image in transmissive geometry We also express with Thus we can simplify Eq. (33) to To compute the contribution of all points on the same z plane we need to compute where o = [o xy , z] is a point on a plane of depth z. Referring to Eq. (44), this is basically the integral of a sinusoid of frequency ω. As the extent of the integration is infinite, the integral is non zero if and only if the first two coordinates of ω are zero, namely the directions satisfy Eq. (37). For such directions we also get following Eq. (40) that the angles betweenî 1 ,v 1 andî 2 ,v 2 are equivalent. To a first order approximation, this angle is justτ . Thus we can write s(î 1 ·v 1 )s(î 2 ·v 2 ) * = ρ(τ ).
where ρ is the phase function. We arrive at A short expansion of Eqs. (53)-(56) leads to Eq. (48). Together with Eq. (47) we conclude at the desired Eq. (38). So far we have computed the single scattering contribution from one z plane. To compute the single scattering correlation of a volume of thickness L we need to integrate the single scattering contribution from all z planes inside the volume, leading to where z o = L/2. Effectively if ω z , the 3rd component of the difference vector is non zero and the material thickness is not zero, integration over multiple z slices will decrease the correlation to zero, in the same way that integration over the x, y dimension leads to zero correlation unless the first 2 components of ω are zero.

A.3 Far field reflective geometry
The reflective geometry case is essentially equivalent to the transmissive case, except for the term d(î→o) + d(o→v).
In reflective geometry, the light arrives from one edge of the material, hits the particle and back-scatters. Thus, the length of the path inside the material is different at different z planes. To this end we denote Assuming the scattering slab is placed at the depth range z ∈ [0, L], it is not hard to show that thus With this notation we can express the single scattering covariance from a single z plane as Integration over all z planes in the slice leads to

A.4 Near field transmissive geometry
We turn to the near field case. We follow the derivation in [20]. This starts from the observation that given a far field, we can compute the field obtained by an illumination source or a camera focused at points i, v inside the sample, using a weighted combination of the far field values at different directions where S 2 is the unit sphere and, assuming an ideal lens, with m(ω) denoting an aperture mask, controlling how much light is passed via aperture directionsω. Using this, they express the single scattering contribution from a scatterer at o as an integral of the far-field correlation terms as in Eq. (33), over all directions in the aperture where Υ denotes integrating the scattering amplitude function over the aperture.
We approximate the viewing aperture function similarly. The scattering function is also approximated as a mixture of von Mises-Fisher functions. For simplicity we will consider here a single mixture component: In this paper we further simplify the von Mises-Fisher model and consider a 2D Gaussian approximation. For that we approximate viewing directions in the aperture as  where z can be a depth inside the slab z ∈ [0, L]. With these notations we can express the integral over illumination and viewing directions, leading to c i 1 ,i 2 v 1 ,v 2 (o) = s · e ao|oxy| 2 +a ∆ |∆ 2 |+aτ |τ | 2 +bo,τ (o·τ )+b ∆,τ (∆·τ ) (77) with where γ G denotes Scattering Tissue The phase of this correlation is and therefore the best tilt angle is When considering the covariance of a thicker volume, spanning the range z ∈ [0, L], we expect that the phase would correspond to some average phase of one of the planes inside the volume. In our derivation in Sec. 4.1 we choose a reference plane of z o = L/2. However, [13] attempts to derive an analytic expression for the tilt shift correlation based on an approximated differential equation. They arrive at a slightly different result, stating that the optimal tilt plane is 2/3 of the way to the end, rather than exactly the mean. Namely, they suggest that the tilt angle should be computed according to the plane z o = (2/3)L rather than according to the plane z o = L/2. The numerical evaluation below suggests that the difference between the 1/2 and the 2/3 rules is really minor. In Fig. 9 we evaluate the effect of the tilt angle numerically using the full version of the MC simulator [19], [20], without relying on the single scattering approximation. We tested the correlation produced by tilt angles of the form θ = ∆ /(zo) for different selections of the plane z o . The highest correlation is obtained when θ is selected at z o = L/2. However, the selection z o = (2/3)L provides results of equivalent quality. Our decision to select the reference plane at z o = L/2 is mostly motivated by the fact that in the far field case, integration can be evaluated in closed-form and the phase indeed corresponds to z o = L/2. Using an isotropic phase function (g = 0), we visualize covariance images C(τ ,∆) as a function of the 2D displacementsτ . Each panel displays 4 different values of vertical illuminator displacements∆ (∆ = 0 corresponds to intensity images). In each case, we compare on the left the correlation produced with the full simulator of [19], and on the right our single scattering approximation of Eq. (20). The comparison includes materials of 3 different thicknesses, corresponding to optical depths 1, 3 and 5. Note that as the displacement between illuminators increases (two lower rows of each panel), the single scattering model is a very good approximation to the correlation, even when the optical depth of the material is much larger than 1. This holds in all 3 panels, visualizing both far and near fields, as well as both reflective and transmissive geometries. At larger displacements∆ we also see the anisotropic support of the correlation. Note that the correlation values are complex and the figure visualizes the absolute value of its real component. As in this example illuminators are displaced vertically, in agreement with Eq. (20) the phase changes along the verticalτ y axis. The real component displays this sinusoidal structure, cycling between low and high values along theτ y direction.