Parametrization of speckle intensity correlations over object position for coherent sensing and imaging in heavily scattering random media

A statistical treatment of speckle correlations as a function of the position of a moving object is shown to provide access to object information through thick and heavily scattering random media. Experiments for a patchlike object of varying size and for varying degree of background scatter are explained using a model, and an experimental study allows evaluation of key attributes. Given a sufﬁcient signal-to-noise ratio, adequate coherence, and developed ﬁeld statistics, measured speckle intensity patterns from a set of object positions can allow high-resolution imaging deep into an obscuring medium and the medium’s scattering strength can be gauged quantitatively with calibration. This enables new opportunities in application domains such as optical sensing, material inspection, and deep tissue in vivo imaging.


I. INTRODUCTION
Coherent optical sensing and imaging methods offer high spatial resolution and spectroscopic information.When coherent light interacts with randomly scattering structures, the resulting constructive and destructive interference between the wave fronts forms speckle, often perceived as problematic.For thin scattering domains, the memory effect preserves information about the incident wave front within a range of incident angles [1], and this has enabled a series of demonstrative imaging experiments [2][3][4].Inversion of the measured transmission matrix [5,6] allows focusing inside a scattering medium.Alternatively, optimizing the incident wave front [7], facilitated by using a guide star in the scattering medium, has enabled focusing to that point [8][9][10].Random speckle intensity patterns also reveal useful statistical information through correlations over appropriate variables.For example, correlations over laser frequency can be used to image hidden inhomogeneities [11] and to characterize the scattering medium [12].The temporal correlation of optical measurements of dynamic material has led to various applications such as diffusing-wave spectroscopy [13].Light in random media has also been a fertile ground for the investigation of localization, a regime where circular Bessel field statistics hold [14].Conversely, while random scatter, measured through a polarizer, generally results in a zero mean circular Gaussian complex field distribution and Rayleigh magnitude, phase control through a spatial light modulator has resulted in synthesized speckle statistics [15].Despite such multifaceted Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
efforts related to light in random media and the important suite of applications, coherent imaging in randomly scattering domains has faced substantial obstacles in reaching beyond one transport length.
It was discovered that the Fourier spectrum of the field incident on a random medium could be extracted from speckle intensity patterns measured in transmission and as a function of spatial translation of the field [16].This understanding that deterministic information surviving heavy scatter through analysis using speckle correlations as a function of incident field motion was extended to the imaging of apertures located between two randomly scattering slabs [17].This led to complicated aperture shapes being imaged with a simple phase retrieval step and the possibility that general objects moving within a random medium might be sensed [18].Here we present a general and rigorous theory for an object moving inside a heavily scattering random medium and evaluate an approximate form using the key degrees of freedom (object size and level of background scatter) with a set of experiments.These results provide insight into a new avenue for basic physical studies of light in random media and for applications related to imaging and material characterization.

II. THEORY
We write the detected (temporal frequency domain phasor) field at r d , measured through a polarizer, as d .The second-order averaged phasor field correlation (at circular frequency ω) over the translated position (from r to r + r) is * d (r) d (r + r) , where • is mathematically an average over the background scatterer reconfiguration while measuring at a point detector (r d ).Experimentally, an average is formed by computing average cross-correlation coefficients using captured speckle intensity data from the pixels of a camera with stationary statistics (where a small region on the surface of the randomly scattering domain is imaged onto the camera).An exact superposition of the field as the result of the illumination and the background randomly scattering environment without the moving object ( db ) and the additional scattered field that fully characterizes the effect of the moving object ( ds ) gives d = db + ds .The second-order field correlation can thus be expanded as * d (r) d (r + r) = I db + I db 1/2 I ds (r + r) 1/2 g (1)  bs ( r) + I ds (r) 1/2 I db 1/2 g * ( 1) bs (0) + I ds (r) 1/2 I ds (r + r) 1/2 g (1)  ss ( r), where (with and is associated with the object, and the normalized secondorder field correlations (with compact argument notation) are defined by g (1)  bs (r, r + r) = * db ds (r + r) / [ I db 1/2 I ds (r + r) 1/2 ]=g (1)  bs ( r), g (1)  sb (r, 0)= * ds (r) db / [ I db 1/2 I ds (r) 1/2 ] = g (1)  sb (0) = g (1) * bs (0), and g (1)  ss (r, r + r) = * ds (r) ds (r + r) /[ I ds (r) 1/2 I ds (r + r) 1/2 ] = g (1)  ss ( r), where field correlations do not depend on the reference position r with the normalizations used.
A Green's function describing the background random scatter can be used to form an expression for ds in terms of a wave-equation-based object function (that mathematically describes the dielectric constant of the moving object, but here is normalized, providing shape information, as treated in a simplified manner previously [18]).Therefore, without the moving object and with only the background field, ds = 0. Assuming that | r| is large compared to the wavelength (and short-range effects are neglected) and that there is adequate random background scatter for d to be Gaussian, so that only the joint spatial support of the object and the translated object contributes to g (1)  ss , we can write g (1)  ss ( r) = dr Õ * (r) Õ(r + r), where Õ is the normalized object function yielding dr Õ * (r) Õ(r + r) = 1 for r = 0. Õ represents a scaled parametrization of the object, proportional to the object's spatial complex dielectric constant distribution.
The normalized intensity pattern at the detector is Ĩd = (I d − I d )/ I d , with I d the mean.With zero-mean circular Gaussian statistics for the field detected through a polarizer, use of a moment theorem [19] leads to where scaling by the means removes the dependency of Ĩd on r.Using (1), we write the expansion of the numerator of (3) as (1)  ss ( r) + D 2 (r, r) g (1)  ss ( r) 2 , ( where D 0 , D 1 , and D 2 are given by D 0 (r, r) = I db 2 + 2 I db 3/2 I ds (r) 1/2 Re g (1)  bs (0) + I ds (r + r) 1/2 Re g (1)  bs ( r) + I db I ds (r) 1/2 I ds (r + r) 1/2 2Re g (1)  bs (0)g (1)  bs ( r) + I ds (r) g (1)  bs (0) 2 + I ds (r + r) g (1)  bs ( r) 2 , ( Re{D 1 (r, r)} = I db I ds (r) 1/2 I ds (r + r) 1/2 + I db 1/2 I ds (r) I ds (r + r) 1/2 Re g (1)  bs (0) + I ds (r) 1/2 I ds (r + r) Re g (1)  bs ( r) , ( Im{D 1 (r, r)} = I db 1/2 I ds (r) 1/2 I ds (r + r) Im g (1)  bs ( r) − I ds (r) I ds (r + r) 1/2 Im g (1)  bs (0) , (7) with Re{•} the real part and Im{•} the imaginary part.Note that D 0 , D 1 , and D 2 contain terms that are not directly obtainable from measurements, for example, I ds and g (1)  bs .Writing (1) for r = 0, we can infer that g (1)  bs (0) must in general be nonzero to account for absorbing or occluding objects, and have a negative real part for a reduction in I d through D 0 .Supported by substantial experimental evidence in which heavy background scatter is present so that I ds can be assumed slowly varying for | r| d, we approximate D 0 , D 1 , and D 2 as constants [assuming a fixed r, so that Equation ( 9) is our key finding, presented in a suitable level of abstraction to treat experimental data with simplicity and a useful parametrization.Each of the mean intensities in the denominator of ( 9) can be measured.Therefore, for known objects, g (1)  ss ( r) can be obtained through (2), and the validity of the approximation that D 0 , D 1 , and D 2 are constants over the length scale corresponding to the object can be evaluated and established.Use of (9) therefore facilitates a tractable inversion of measured intensity correlations to form a geometric image of a hidden, moving object.
By performing a nonlinear least-square fit based on (9), we find excellent agreement with the experimental results we have, as we now describe.From (2), g (1)  ss decreases from 1 at r = 0 to 0 at r = d.Thus, (9) indicates that D 0 is revealed by the minimum of the speckle intensity decorrelation at d, the size of the object.The two other constants in (9) (D 1 FIG. 1.The experimental setup for measuring the speckle intensity correlation of a moving circular patch translated along the y axis between two scattering layers.Speckle patterns are collected at the camera after passing through a spatial filter and a polarizer.and D 2 ) can be determined by fitting the experimental data using a known object function.Supported by both theory and experimental results, we found that D 1 is small relative to both D 0 and D 2 .By using ( 6), (7), and ( 8 + I db 1/2 I ds (r + r) −1/2 Re g (1)  bs (0) + I ds (r) −1/2 Re g (1)  bs ( r) , (10) bs ( r) bs (0) .( 11) With substantial background scatter between the object and the detector and a relatively strongly scattering object, I db and I ds are expected to be a similar order of magnitude.For a nonaperture type object, Re{g (1)  bs } < 0. Jointly, this leads to the conclusion that |D 1 | D 2 .Also, for situations where I ds I db , D 1 is negligible.However, with fixed background scatter and small object scatter, such that I ds is small relative to I db , D 1 /D 2 will increase.These conclusions are supported by (10) and (11), as well as the experimental data.With negligible D 1 , (9) becomes Ĩd (0) Ĩd ( r) Equation ( 12) produces accurate predictions for a strongly scattering moving object for which a relatively large I ds is expected.Consequently, a simple re-normalization of ( 12) after subtraction of D 0 (indicated by the minimum) provides direct access to |g (1)  ss ( r)| 2 , and with phase retrieval, imaging becomes possible.For a weakly scattering object in a strongly scattering environment, (9) should be used.We consider the variables related to application of ( 9) and ( 12) based upon experimental data, under the assumption that Õ is real, and how D 0 , D 1 , and D 2 change in various situations.

III. EXPERIMENT
Figure 1 shows our experiment.A 59-mW 850-nm laser diode with a linewidth less than 10 MHz was used for illumination.Two layers of scattering material are separated by a small distance (about 5 cm) to allow a pair of stages to move the objects of interest in the transverse two-dimensional plane between these scattering slabs.A 4F system is used to filter the resultant speckle patterns, so that the camera pixels have adequate resolution and there are a sufficient number of random samples.Detection is through a polarizer, to provide circular Gaussian scalar field statistics.An area of approximately 1.8 mm × 1.8 mm on the back of the second scattering layer is imaged by the camera using magnifying optics.The speckle images have stationary statistics that are used to form averages.Also, to improve the estimates, multiple intensity correlation coefficients that correspond to the same r but different r are averaged to form an estimation of the ensemble average over many scatterer reconfigurations.
In our experiments, ground glass diffusers and acrylic slabs [shown in Fig. 2(a)] were used as the scattering material.The acrylic slabs contain TiO 2 scatterers that have a mean diameter of 50 nm, with a reduced scattering coefficient estimated to be μ s = 4 cm −1 , and negligible absorption.Stacks of ground glass diffusers also provide heavy scatter and the measured speckle patterns become highly decorrelated over object displacement, evident in Fig. 2(b).Five different scattering layers were used: 3-and 6-mm-thick acrylic slabs (14 cm × 14 cm), a stack of four ground glass diffusers, a stack of six ground glass diffusers (each individual ground glass piece was 10 cm × 10 cm and 0.2 cm thick, and the stacks were taped together at the edges), and a single ground glass slide (1500 grit).The four-piece stack consisted of two 120-grit and two 1500-grit ground glass slides and the six-piece stack consisted of two 120-grit and four 1500-grit ground pieces.The moving objects were circular absorptive patches of different diameters (3.7, 5, and 6 mm).These patches were formed with adhesive black tape attached to transparent plastic windows (10 cm × 13 cm × 0.15 cm), making for a binary object (either completely absorptive or transparent).

IV. RESULTS
Figure 3(a) shows the measured normalized speckle correlations for varying object (patch) size and hence scattering strength, with fixed scatter on the laser and detector sides.For all data we present, the (wavelength-scale) short-range intensity correlation has been neglected (see, e.g., previous results for this super-resolution regime [17]) and the macro-FIG.3. (a) Measured intensity correlations over object position with various absorbing patch sizes with a fixed slab configuration (a four-piece ground glass diffuser at the laser side and a 6-mm-thick μ s = 4 cm −1 acrylic slab at the detector side) and circular patches having diameters of 3.7, 5, and 6 mm.The larger the size of the object, the deeper the decorrelation dip.(b) The scaled decorrelation of each circular patch [using the numerator minus D 0 in ( 9)] agrees well with our prediction, with the use of ( 2) and ( 9), for r smaller than the object's size.From nonlinear least-squares fitting, Re{D * 1 } is negligible for the largest circular patch (6 mm in diameter), and more than one order of magnitude smaller than D 2 for the smaller patches.scopic measurement data has been extrapolated to r = 0 for normalization.A four-piece ground-glass-diffuser stack was used on the laser side and a 6-mm-thick acrylic slab on the detector side.The motivation for using the glass diffuser stack on the laser side was to minimize possible decorrelations due to heating, as might occur with the acrylic scattering slabs.As Fig. 3(a) indicates, the normalized intensity correlations over scanned object position decrease from unity, reach a minimum, and then increase (ultimately to close to 1 again, because the speckle patterns become more similar when the object is displaced distances large relative to its size).The larger the object, the greater the dip in the intensity correlation.This can be understood as the increasing patch size producing larger I ds and D 2 becoming more dominant.The minima in the correlations occur in each case at the diameters of each circular patch and reveal the values of D 0 from (12). Figure 3(a) indicates that I ds / I db increases with patch size and hence object scattering strength.This is reflected in the lower minimum (D 0 ) with increasing patch size.
Computing the numerator of the right-hand side of (9) using measured data, subtracting D 0 (from the respective minima), and rescaling so the result is unity at zero displacement, we obtain the dotted points in Fig. 3(b) for the three patch sizes (shown by the symbols).With prior information about g (1)  ss , fitting the processed data according to (2) using 2Re{D * 1 g (1)  ss ( r)} + D 2 |g (1)  ss ( r)| 2 and rescaling in the same way, we obtain the predictions shown as the solid lines in Fig. 3(b).Over the regime where our assumptions hold and the speckle is decorrelating, the agreement between theory and experiment is excellent.For the 6-mm-diameter patch, Re{D * 1 } is negligibly small.For the two smaller patches, Re{D * 1 } values, obtained from nonlinear least-square fitting, are more than one order of magnitude smaller than D 2 (assuming Õ a real function).
Figure 4 shows experimental intensity correlation results for the 3.7-mm-diameter patch and varying levels of scatter on the detector side [(Figs.4(a) and 4(c)] and the laser side  9)] for different amounts of scatter on the detector side ("gg" represents the number of ground glass slides used, and "mm" indicates the thickness of scattering acrylic slab).(b) and (d) show the measured and rescaled correlations for different amounts of scatter on the laser side.The scaled speckle correlations agree well with our prediction using ( 2) and ( 9) for r smaller than the object's size.With nonlinear least-square fitting, Re{D * 1 } is negligible except for when the four ground glass slides and 6-mm-thick acrylic slab were used (presenting a heavily scattering environment).
[Figs.4(b) and 4(d)].For the results shown in Fig. 4(a), the scattering layer on the laser side was fixed (four-piece glass diffuser stack) while the scattering layer on the detector side was varied (six-piece diffuser stack, 3-mm acrylic slab, and 6-mm acrylic slab, in order of increasing scatter).Increasing the amount of background scatter on the detector side reduces the dip in the intensity correlation and hence makes D 0 more prominent relative to D 2 .Using the same data processing approach described for Fig. 3(b), we show in Fig. 4(c) how well the theory described in (9) matches the experimental results, with the use of (2).In another set of measurements, fixing the scatter on the detector side (6-mm-thick acrylic slab) and varying the scatter on the laser side (the single ground glass diffuser, the four-piece stack of diffusers and the 3-mm-thick acrylic slab, in order of increasing scatter) produced the results in Fig. 4(b).Increasing the background scatter on the laser side also increases D 0 /D 2 .The exception is the 3-mm-thick acrylic slab, for which heating becomes the dominant source of change, overshadowing the displacement of the small circular patch.We show in Fig. 4(d) that the prediction from (2), with the use of ( 9), again matches the experimental results nicely (except for the 3-mm-thick acrylic slab).With nonlinear least-square fitting, the curves in Figs.4(c) and 4(d) have negligible Re{D * 1 } except for the case with more scattering four ground glass pieces and 6-mm-thick acrylic slab, where Re{D * 1 } is still more than one order of magnitude smaller than D 2 , indicating that the representation in (12) is suitable.This can again be understood using the comparisons in (10) and (11), where a heavily scattering environment reduces the detected intensity without the moving object, denoted as I db .
(13) Equation ( 13) describes our conclusion from the two sets of experiments shown in Fig. 4, where with an increase in background scatter on either the laser or detector size, D 0 /D 2 increases, and from (13) this implies an increase in I db / I ds .An increase in the background scatter increases the spatial spread of the speckle pattern exiting the scattering medium.With the transmission arrangement of Fig. 1, where the spot being imaged is fixed in size, this increase in background scatter results in a smaller detected mean intensity.However, the intensity correlation data presented is normalized, thereby accounting for this reduction.Increasing the amount of scatter on either the laser or the detector side has a similar impact in the results of Fig. 4, but we should consider the impact on I db and I ds .The distance between the moving object and the detection spot, |r d − r| = r do , is less than that between the laser excitation spot on the scattering medium and the detector, |r d − r s | = r ds .Thus, the relative influence of I db ∼ r −2 ds and I ds ∼ σ r −2 so r −2 do , with σ describing the scattering cross section of the moving object and r so the distance from the source to the object, scale differently.With a more absorbing patch, σ reduces and the measured intensity reduces, and in the experiments this occurs with an increase in absorbing patch size.For the axial situation with r do = αr so and r ds = r so + r do , we find In the experiment, the background scatterer density is considered fixed.With larger thickness of the scattering medium on the detector side (with r so fixed), the increasing α results in a smaller I ds / I db .When the thickness of the scattering medium on the laser side is increased (with r do fixed), the resulting decreasing α also leads to a reduction in I ds / I db .
We have thus established consistency between the model and the experimental results shown in Fig. 4, both in terms of the coefficients (D 0 and D 2 ) and how their ratios relate to the mean intensities used in the model.The model we present based on (9) describes decorrelation of the speckle patterns with motion of the object.In the special case of an aperture in a screen, this model contracts to that presented earlier [17,18].The impact of changing the scattering strength of the moving object and the surrounding environment, predicted by ( 9) with ( 5)-( 8), is experimentally verified through the results shown in Figs.3(a), 4(a), and 4(b), demonstrating quantitative and predictive character.Moreover, the agreement between our experimental results and the predictions using ( 2) with (9), shown in Figs.3(b), 4(c), and 4(d), indicates that in such situations D 0 , D 1 , and D 2 can be treated as constants.This leads to a simple means to reconstruct an object function by fitting the values of the constants (D 0 , D 1 , and D 2 ) and performing phase retrieval using the Fourier magnitude of the object obtained from its autocorrelation.The dimensions available for this reconstruction are commensurate with those from object movement.Also, the ratio between the constants D 0 , D 1 , and D 2 is indicative of the relative strengths of the scattering environment and the moving object.This could lead to applications in which the character of the background scattering environment is sensed with a calibrated measurement, and direct comparison of the scattering strengths of different hidden objects is also feasible.
There is an assumption that the background random scatterers are fixed or that their motion (whether displaced or natural) is either negligible or can be calibrated.We expect that (9) will also explain the increase in the speckle intensity correlations of Figs. 3 and 4 after the minima.Considering (9) with ( 5)- (8), this is explained by the dominance of D 0 and is a regime where g (1)  ss [hence the second term in the numerator of (9)] and D 2 [and therefore the third term in (9)] are both small.In the limit of large r, D 0 (r, r) → I db 2 and thus Ĩd (r) Ĩd (r + r) → 1.For this regime of larger translation distance relative to the object size, the r dependency of D 0 must be considered, and this domain may thus yield more information about the nature of (5).

VI. CONCLUSION
The use of speckle correlations over object position provides for imaging inside and through an unprecedented level of background random scatter, subject to detector noise and extraneous motion.It also allows probing of the amount of environmental scatter and quantification of the scattering strength of the hidden moving object.With speckle measurements over small translation distance, far-subwavelength spatial information about the moving object becomes available [17], indicating the promise of motion in structured illumination for super-resolution imaging [20].The measurement can be adapted into a reflection geometry for a wider range of practical applications.The velocity or the relative position of the moving object could be estimated by the speckle temporal contrast, the Doppler shift [21], or use of localization in a diffusion model [22].While fluorescent imaging has proven useful in biological samples, our model allows for coherent imaging at high resolution (larger than one wavelength, with neglect of short-range correlations) and without the need to introduce fluorescent reporters.By combining the accuracy of localization by emission [22] and coherent imaging based on motion, complementary and supporting information becomes available.For example, one could track a moving cellular cluster labeled with quantum dots inside deep tissue using emitted diffusive fluorescent light and image at high spatial resolution using intensity correlations over object with the coherent excitation light, as described here, to achieve in vivo imaging of biomolecules.

FIG. 2 .
FIG. 2. Example heavily scattering material used in our experiments.(a) Two 6-mm, μ s = 4 cm −1 slabs of acrylic are placed on top of a page with printed stripes.Through one 6-mm-thick slab, one can no longer distinguish individual stripes, and through a total of 12 mm, it is not possible to distinguish the striped area.(b) The centrally cropped speckle patterns for a moving 4-mm circular patch placed between a 4-ground-glass stack and a 6-ground-glass stack are highly decorrelated over an object displacement of about 2.5 mm.

FIG. 4 .
FIG.4.Measured intensity correlations over object position with various amounts of scatter on either side of a 3.7-mm-diameter circular absorbing patch.(a) and (c) show the measured correlations and the rescaled correlations [with the numerator minus D 0 in(9)] for different amounts of scatter on the detector side ("gg" represents the number of ground glass slides used, and "mm" indicates the thickness of scattering acrylic slab).(b) and (d) show the measured and rescaled correlations for different amounts of scatter on the laser side.The scaled speckle correlations agree well with our prediction using (2) and (9) for r smaller than the object's size.With nonlinear least-square fitting, Re{D * 1 } is negligible except for when the four ground glass slides and 6-mm-thick acrylic slab were used (presenting a heavily scattering environment).