## Abstract

Random scattering of light by a turbid layer prevents conventional imaging of objects hidden behind it. Angular correlations in the scattered light, created by the so-called optical memory effect, have been shown to enable computational image retrieval of hidden sources. However, basic memory-effect imaging contains no spatial (*x*) information, as only angular (*k*-space) measurements are made. Here, we use windowed Fourier transforms to record scattered-light images in the full {*x*,*k*} phase space. The result is the ability to discriminate size and depth of individual sources that are hidden behind a thin scattering layer.

© 2014 Optical Society of America

## Introduction

The scattering of light by optically inhomogeneous random media is encountered in many imaging contexts, such as clouds and fog in photography, atmospheric turbulence in astronomy, and biological tissue in microscopy. Typically, scattering ruins image quality, degrading the correspondence of the spatial distribution of detected light with that of its source. Progress in adaptive optics, wavefront shaping, and computational imaging has led to methods of reestablishing correspondence, both optically and digitally [1–7]. Methods to overcome scattering rely on understanding of the scattering process to varying degrees. At one extreme, knowledge of the medium can be leveraged to undo the effect of scattering, for example through inversion of the transmission matrix [8], phase conjugation [9], or nonlinear instability [10,11]. However, in cases where microstructural knowledge is unobtainable or measurement is impractical, correlative or statistical properties of the scattered light may still be sufficient to enable image recovery [12].

Light scattered by random media retains information about the illuminating source in the form of correlations. One remarkable example is the optical memory effect [13,14] in which the angular spectrum of the scattered light co-varies with the angle of incident light, a correlation which persists even into the multiply scattered regime. The memory effect has enabled imaging through scattering by optical phase correction with wavefront shaping [15,16], as well as by computational image retrieval with iterative algorithms adapted from classic techniques in astronomy [17,18].

A limitation of image retrieval studies based on memory-effect correlations alone is that the recovered images have been functions of illumination angle only. As such, pixels are intensities associated with angles of incidence relative to the scattering surface, and this data alone does not specify where in space the light originates. The problem has a simple analogy to that of perspective in photography, in which each pixel collects parallel rays of a given direction. In such photographs, the human viewer can interpret relative size and distance based on the content of the image, such as location of the horizon, but the image acquisition itself cannot separate these factors. Here, we supplement the intensity image retrieved from memory-effect correlations with phase-space measurement [19] of spatio-angular correlations to discriminate axial depth, analogous to taking photographs from multiple viewing angles, to obtain the distance and size of the source.

## Phase-space measurement of scattered light

As shown in Fig. 1, we construct an {*x*,*k*} phase-space microscope using a spatial light modulator (SLM). As a simple lens at the f-f Fourier condition would transform all position coordinates {x} to angular coordinates {k}, we use amplitude modulation on the SLM to produce windowed (local) Fourier transforms of the field in the microscope focal plane (Fig. 1a). More specifically an amplitude mask on the SLM blocks light outside of a transmission window (details can be found in Ref [19].). The field within the modulated region is optically Fourier transformed by a 2f lens, and the intensity in the Fourier plane (*k*-space) is imaged by a CCD sensor (7.4 μm pixel size, 4 MP, Photometrics CoolSNAP K4). A second CCD sensor placed in a 4f image plane of the SLM permits imaging of the spatial (*x*-space) distribution of intensity within the windowed region. As a demonstration of a square-windowed Fourier transform, a collimated laser beam produces the image of the window on the *x*-space sensor and the corresponding far-field intensity of a plane wave on the *k*-space sensor (Fig. 1b). After inserting a ground glass diffuser (120-grit, Thorlabs DG20-120) in the focal plane, we observed scattering of the beam as coherent speckle in *k*-space (Fig. 1c). Phase variation due to the diffuser randomizes the correspondence between the *x*-space and *k*-space intensities. Without microstructural knowledge of the diffuser, the *k*-space intensity cannot be predicted from the *x*-space intensity alone, and without compensation or processing, the random distortion of the wavefront prevents direct imaging.

## Correlations in scattered light

To demonstrate memory-effect correlations between neighboring source points, we focused a laser through a long working distance, 0.7 NA objective to a focus 1mm behind the diffuser, creating a point source which we then translated laterally on a manual stage. Lateral translation produced a corresponding shift in the *k*-space speckle (Fig. 2(a)). We compared the cross-correlation of the shifted and reference speckle with the results of numerical simulation of wave scattering by a simulated diffuser (Fig. 2(b)). Theoretically, the normalized cross-correlation of the speckle produced by incident illumination differing by an angle, $\Delta \theta \approx \lambda \Delta k$, depends on the magnitude of the difference, $\left|\Delta k\right|$, and the thickness of the scattering medium, *L*, according to $C\left(\left|\Delta k\right|\right)={\left(\left|\Delta k\right|L/\mathrm{sin}h\left(\left|\Delta k\right|L\right)\right)}^{2}$ [13]. In our simulation, we used an idealized scattering layer with *L* = 0; correspondingly, the simulated speckle remained fully correlated (*C* = 1) along the line $k={x}_{o}/\left(\lambda z\right)$. Experimentally, we found that *C* declined to 80% over a source translation distance of 300μm, corresponding to an angular shift of 0.5μm^{−1}. Thus, the random speckle remains strongly correlated between source points separated by hundreds of microns, though subject to noise from speckle degradation. Shift-invariance of the speckle with respect to source position implies that the *k*-space intensity, *I _{k}*, of scattered light from a planar, incoherent source can be represented as a convolution of the source shape,

*O*, with the intensity pattern of coherent speckle from an individual point source,

*S*, so that ${I}_{k}=O\otimes S$. As was noted, however, this convolution occurs in

*k*-space, which mixes the independent factors of lateral position and axial depth.

Correlations that depend on propagation distance, and therefore axial depth alone, can be detected in phase space. Whereas memory-effect correlations originate in the scattering process itself, correlations arising from the spatial coherence of the incident light prior to scattering may be retained in the scattered light. To demonstrate preservation and detection of spatial coherence in the scattered light, we scanned the phase-space window over the surface of the diffuser. For a fixed position of the source, lateral translation of the window produced a shift in the average *k*-space distribution, revealing partial spatial coherence retained through scattering (Fig. 2(c)).

Spatial coherence can be represented formally by the equal-time, two-point correlation of the complex field, $J\left({r}_{1},{r}_{2}\right)=\u3008E\left({r}_{1}\right){E}^{\ast}\left({r}_{2}\right)\u3009$, where *J* is the mutual intensity and $\u3008\u3009$ denotes an ensemble average over the statistical field, *E*(r), defined over the spatial coordinate, *r* [20]. Our phase-space measurement does not sample the mutual intensity directly but rather measures its equivalent representation in phase space, the optical Wigner distribution function (WDF), $W\left(x,k\right)={\displaystyle \int J\left(x-\xi /2,x+\xi /2\right)}{e}^{ik\xi}d\xi $ [21], which has the form of a Fourier transform of *J* over the transverse spatial variable *ξ* centered on *x*. A primary advantage of working in phase space is that parameterization in spatial and angular coordinates, {*x*,*k*}, makes phase space a natural setting for describing and analyzing spatio-angular correlations. In fact, the scanning window used here samples a filtered version of the WDF, due to the finite extent of the aperture, known as the spatial spectrogram. The spectrogram, *P _{S}*, can be derived from the WDF by filtering with an appropriate phase-space kernel [22] through a double convolution over spatial and angular coordinates, so that ${P}_{S}\left(x,k\right)={\displaystyle \iint H\left(x-x\text{'},k-k\text{'}\right)W\left(x\text{'},k\text{'}\right)dx\text{'}dk\text{'}}$, where

*H*is the WDF of the window function. Implicit in this convolution is a tradeoff between spatial and angular resolution which depends on the shape of the window. Window size determines both the spatial and angular width of the filter kernel,

*H*. In particular, the finite spatial width of the window limits angular resolution in

*k*-space measurement, and limited

*k*-space resolution can lead to blurring of small, distant objects which require fine angular discrimination. With the square window, the limit on object lateral resolution is$\Delta {x}_{o}>\lambda z/D$, which for an axial depth

*z*~1mm and window size

*D*= 64μm limits image recovery to planar sources of sizes $\underset{\u02dc}{>}$ 10μm.

We simulated the phase space of the scattered light by calculating *P _{S}* for the field at the diffuser surface. In the absence of scattering,

*P*of the freely propagating field is a sheared line with a slope determined by the axial distance of the source to the focal plane. In both simulation and the experimental data, we found scattering leads to vertical (

_{S}*k*-axis) line broadening, but the slope of the average

*P*is unchanged, reflecting the depth of the point source behind the scattering layer (Fig. 2(d)).

_{S}## Image retrieval from angular correlations

To demonstrate image retrieval of incoherent samples at different depths, we created an incoherent source in the shape of a letter by imaging a lamp-illuminated “P” or “U”, separately measured in independent experiments, to planes at ~1.4mm and ~1mm, respectively, behind the surface of the diffuser (Fig. 1(a)). Individual images taken by the CCD sensor displayed complex *k*-space intensity, *I _{k}*, resulting from convolution of the source shape with random speckle from the diffuser (Fig. 3(a)). By moving the phase-space window to a non-overlapping region of the diffuser, we obtained statistically independent speckle. We combined images from up to 15 window positions to estimate the ensemble-averaged autocorrelation, $\u3008{I}_{k}\ast {I}_{k}\u3009$ (Fig. 3(b)).

The convolution relationship $\u3008{I}_{k}\ast {I}_{k}\u3009=\left(O\ast O\right)\otimes \u3008S\ast S\u3009$ shows that the ensemble-averaged autocorrelation of speckle, $\u3008S\ast S\u3009$, acts as a linear filter on the object autocorrelation, $O\ast O$, to produce the average autocorrelation of *k*-space intensity, $\u3008{I}_{k}\ast {I}_{k}\u3009$. Though *S* is random for each individual image, $\u3008S\ast S\u3009$ is a sharply peaked function that allows for direct estimation of $O\ast O$ through measuring $\u3008{I}_{k}\ast {I}_{k}\u3009$. When described in the Fourier domain, this method is conceptually identical to the speckle interferometry developed by Labeyrie for stellar imaging through atmospheric turbulence [23].

Through the Weiner-Khinchin theorem, the power spectral density of the source distribution, ${\left|\tilde{O}\right|}^{2}$, can be obtained from the Fourier transform of $O\ast O$. In order to specify the source distribution, the Fourier modulus, $\left|\tilde{O}\right|$, requires missing phase information. The phase can be recovered computationally through Fienup-type retrieval algorithms which enable iterative reconstruction of phase from the Fourier modulus and satisfaction of constraints, typically that intensity is real-valued and non-negative [24]. We computationally recovered the source distribution from the average autocorrelation by an iterative phase retrieval procedure slightly modified from that of Bertolotti, *et al*. [17] (Fig. 3(c)). Comparison of the retrieved images with images directly obtained of the source letters showed good agreement of shape and relative contrast (Fig. 3(d)). However, the quantitative size of the retrieved image cannot be determined without knowledge of source distance.

## Phase-space measurement of spatio-angular correlations

By itself, the spatial (*x*-space) distribution of intensity of the scattered light from the hidden source does not reveal source shape or source depth (Fig. 4(a)). However, as described above, the spatial dependence of the angular spectrum, captured in the spatial spectrogram, *P _{S}*, is determined by the depth of the hidden source. To measure

*P*of the scattered light, we recorded the

_{S}*k*-space intensity as a function of window position (Fig. 4(b)). This measurement includes contributions from both the heterogeneity of random scattering as well as spatio-angular correlations in the incident light. The variation from random scattering adds uncertainty to the measurement of the correlations, which we addressed by acquiring and fitting to data from many window positions, as well as averaging over an additional independent parameter (the

*y*coordinate in this case) to statistically reduce the effect of random tilts of individual scattered wavefronts. We analyzed

*P*by plotting the line profile of intensity along the

_{S}*k*-axis against the

_{x}*x*position of the window (Fig. 4(c)-4(d)). To obtain the slope, we averaged

*k*-axis intensities for the same

_{x}*x*coordinate over 5 windows with different

*y*coordinates and fit the average

*k*-axis intensity with a double-Gaussian function to determine the location of the peak. The relation between peak location and window position was well fit by a straight line, consistent with the correlation expected from propagation-induced shear of

_{x}*P*from a point source (Fig. 4(e)). Although the source distribution is not a single point in our case, the presence of additional emitters did not affect the individual correlations from each source point; the observed correlation is thus the superposition or average of the individual correlations. Though the linear correlation broadens over the angular extent of the source, thus limiting the resolution in depth, the

_{S}*k*-coordinate of the peak average intensity allowed determination of the correlation slope. This method is analogous to spatial point-source localization, as in astronomy and single-molecule imaging, generalized to the {

*x*,

*k*}-coordinates of phase space.

## Depth and size measurement of hidden source

From the fitted slope, *a*, we calculated the source distance from the diffuser as $z=1/\left(a\lambda \right)$. The calculated values (0.9 mm and 1.3 mm) match well the known source displacements (1 mm and 1.4 mm). In turn, the source distance, *z*, allows calibration of the angular image data. Converting angular pixel size, *k _{p}*, to spatial pixel size,

*x*, according to ${x}_{p}={k}_{p}\lambda z$, and resizing the recovered image allows a direct comparison of the size of the source with that of the image retrieved from its scattered light (Fig. 5).

_{p}## Conclusion

Phase-space measurement of scattered light provides simultaneous spatial and angular resolution, which is the essential feature for resolving spatio-angular correlations for depth and size calibration. In the experiments above, the spatial dependence of the angular spectrum resulted from the spatial coherence of the incident light. Each non-overlapping region of the diffuser was an identically distributed and statistically independent scattering configuration, which permitted us to form averages and recover an image from the data acquired in a single phase-space measurement. The statistics of scattering media need not be spatially homogeneous, and one potentially interesting application of this technique would be in performing image recovery in media of spatially varying scattering statistics.

Although we considered image recovery of only individual planar sources, the ability to resolve depth raises the possibility of separating source planes in a composite setting. Solving the 3-dimensional problem will likely require new computational tools, in particular computational methods operating in phase space.

## References and links

**1. **E. N. Leith and J. Upatnieks, “Holographic imagery through diffusing media,” J. Opt. Soc. Am. **56**(4), 523 (1966). [CrossRef]

**2. **W. B. Bridges, P. T. Brunner, S. P. Lazzara, T. A. Nussmeier, T. R. O’Meara, J. A. Sanguinet, and W. P. Brown Jr., “Coherent optical adaptive techniques,” Appl. Opt. **13**(2), 291–300 (1974). [CrossRef] [PubMed]

**3. **I. Freund, “Looking through walls and around corners,” Physica A: Statistical Mechanics and its Applications **168**(1), 49–65 (1990). [CrossRef]

**4. **I. M. Vellekoop and A. P. Mosk, “Focusing coherent light through opaque strongly scattering media,” Opt. Lett. **32**(16), 2309–2311 (2007). [CrossRef] [PubMed]

**5. **Z. Yaqoob, D. Psaltis, M. S. Feld, and C. Yang, “Optical phase conjugation for turbidity suppression in biological samples,” Nat. Photonics **2**(2), 110–115 (2008). [CrossRef] [PubMed]

**6. **N. Ji, D. E. Milkie, and E. Betzig, “Adaptive optics via pupil segmentation for high-resolution imaging in biological tissues,” Nat. Methods **7**(2), 141–147 (2010). [CrossRef] [PubMed]

**7. **S. Popoff, G. Lerosey, M. Fink, A. C. Boccara, and S. Gigan, “Image transmission through an opaque material,” Nat Commun **1**(6), 81 (2010). [CrossRef] [PubMed]

**8. **S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, “Measuring the transmission matrix in optics: an approach to the study and control of light propagation in disordered media,” Phys. Rev. Lett. **104**(10), 100601 (2010). [CrossRef] [PubMed]

**9. **A. Yariv, “Phase conjugate optics and real-time holography,” IEEE J. Quantum Electron. **14**(9), 650–660 (1978). [CrossRef]

**10. **D. V. Dylov and J. W. Fleischer, “Nonlinear self-filtering of noisy images via dynamical stochastic resonance,” Nat. Photonics **4**(5), 323–328 (2010). [CrossRef]

**11. **D. V. Dylov, L. Waller, and J. W. Fleischer, “Nonlinear restoration of diffused images via seeded instability,” IEEE J. Sel. Top. Quantum Electron. **18**(2), 916–925 (2012). [CrossRef]

**12. **J. W. Goodman, *Speckle Phenomena in Optics: Theory and Applications* (Roberts & Company, 2010).

**13. **S. C. Feng, C. Kane, P. A. Lee, and A. D. Stone, “Correlations and fluctuations of coherent wave transmission through disordered media,” Phys. Rev. Lett. **61**(7), 834–837 (1988). [CrossRef] [PubMed]

**14. **I. Freund, M. Rosenbluh, and S. Feng, “Memory effects in propagation of optical waves through disordered media,” Phys. Rev. Lett. **61**(20), 2328–2331 (1988). [CrossRef] [PubMed]

**15. **I. M. Vellekoop and C. M. Aegerter, “Scattered light fluorescence microscopy: imaging through turbid layers,” Opt. Lett. **35**(8), 1245–1247 (2010). [CrossRef] [PubMed]

**16. **O. Katz, E. Small, and Y. Silberberg, “Looking around corners and through thin turbid layers in real time with scattered incoherent light,” Nat. Photonics **6**(8), 549–553 (2012). [CrossRef]

**17. **J. Bertolotti, E. G. van Putten, C. Blum, A. Lagendijk, W. L. Vos, and A. P. Mosk, “Non-invasive imaging through opaque scattering layers,” Nature **491**(7423), 232–234 (2012). [CrossRef] [PubMed]

**18. **O. Katz, P. Heidmann, M. Fink, and S. Gigan, “Non-invasive single-shot imaging through scattering layers and around corners via speckle correlations,” Nat. Photonics **8**(10), 784–790 (2014). [CrossRef]

**19. **L. Waller, G. Situ, and J. W. Fleischer, “Phase-space measurement and coherence synthesis of optical beams,” Nat. Photonics **6**(7), 474–479 (2012). [CrossRef]

**20. **M. Born and E. Wolf, *Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light* (New York: Cambridge University Press, 2002).

**21. **M. J. Bastiaans, “The Wigner distribution function applied to optical signals and systems,” Opt. Commun. **25**(1), 26–30 (1978). [CrossRef]

**22. **L. Cohen, “Time-frequency distributions – a review,” Proc. IEEE **77**(7), 941–981 (1989). [CrossRef]

**23. **A. Labeyrie, “Attainment of diffraction limited resolution in large telescopes by Fourier analyzing speckle patterns in star images,” Astron. Astrophys. **6**, 85–87 (1970).

**24. **J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. **21**(15), 2758–2769 (1982). [CrossRef] [PubMed]