Optical imaging through turbid media has various applications, from long-range observations through turbulence to imaging inside living tissues [1–5]. In turbid media, light waves are distorted due to inhomogeneities, resulting in a loss of spatial resolution and reduced imaging depth. Low-order aberrations due to random refractive index fluctuations can be overcome by adaptive optics [6] and turbulence-free ghost imaging [7], while the problem becomes intractable for optically opaque media in which strong light scattering scrambles the spatial information conveyed by light fields. Early experiments using holographic imaging [8] demonstrated that scattering by stationary media does not erase spatial information carried by light fields [9]. Speckle patterns appear random but are essentially deterministic, and information about the optical input can be retrieved. Recent advances in wavefront shaping exploiting transmission matrices [10–14] and optical phase conjugation [15–19] have enabled focusing and imaging through scattering media. However, they are either invasive and require holographic or interferometric measurements or need prior knowledge of the scattering properties of the media. Recent developments taking advantage of the angular correlation of speckle patterns [20–22], i.e., “memory effect,” enable non-invasive imaging through scattering layers using the auto-correlation of speckle patterns and a phase retrieval algorithm [23–26]. Even though these techniques do not rely on the scattering properties of the media, they share some shortcomings. For example, the memory-effect range restricts this approach to thin scattering layers, the small single speckle grain requires a high-resolution camera to resolve, and the iterative phase-retrieval algorithm suffers from falling into local optimal solutions.

In addition, the dynamic nature of some media, such as fog and biological tissues, introduces another aspect to the challenge: the mapping between input and output fields becomes time dependent, resulting in a rapid decorrelation of optical output information. The vision through the turbid medium between the object and camera has been realized through speckle correlation and shower curtain effects [27]. However, this method fails when the medium between the source and the object is also dynamic. Alternatively, ultrasound-modulated light correlation can be used to image an object hidden inside dynamic media [28]. Nevertheless, the measurements are carried out using raster scans, thus causing a very slow imaging speed. Therefore, fast imaging modalities of objects completely immersed in dynamic media are still in high demand.

In this study, we describe a technique for computational imaging, called spatial–temporal encoded pattern (STEP) illumination, which allows noninvasive imaging through scattering media. We demonstrate experimentally that an image of an object can be reconstructed from a 1D time series of light intensity measured by a photodetector (PD), using ground glass diffusers and slices of chicken breast (1.2 mm thick each) as scattering media.

Enlightened by the concepts of intensity modulation and Fourier transform (FT)-based discrimination [29], we design a sequence of patterns that consists of a bundle of sinusoidal time series with different frequencies. Every spatial location in each pattern is encoded by a unique frequency, a unique feature of the periodic oscillation of the pixel values along the “time” axis (i.e., looking at a single spatial location through different patterns). We illuminate the diffuser–object system with this sequence of patterns and collect the transmitted light by a single-pixel detector. The images of the objects are retrieved based on the fast FT (FFT) algorithm without knowledge of the object or scattering medium.

The design of STEP and the principle of image reconstruction are sketched in Fig. 1. We generate a sequence of grayscale (8-bit, 256 pixel values) patterns of height $H$ and width $W$ ($H\times W$ matrix). For a spatial location $(i,j)$ ($i$th row, $j$th column), the time series is given by $$a_{ijt}=127.5\sin \left(2\pi f_{ij}\frac{t}{r_s}\right)+\mathrm{127.5,}t=\mathrm{0,1,}\dots {,M}_{\mathrm{patt}}-1,$$(1)in which $r_s$ is the sampling rate, $t$ is the discrete time variable (index of the patterns), and $M_{\mathrm{patt}}$ is the total number of patterns. The frequency $f_{ij}$ is defined by $f_{ij}=f_0+(j+iW)\mathrm{\Delta }f$, where $i=\mathrm{0,1,}\dots ,H-1,j=\mathrm{0,1,}\dots ,W-1$ are the row and column indices of the pixels, respectively. $f_0$ is the starting frequency, and $\mathrm{\Delta }f$ is the increment. To avoid signal aliasing, we set $r_s=8f_{\max }$, where $f_{\max }=f_0+(HW-1)\mathrm{\Delta }f$. The patterns are successively projected to the diffuser–object system and synchronically collected by the single-pixel PD. The 1D time series intensity can be expressed by $$I_t=\sum _{i,j}{a}_{ijt}{I}_{ij}+N_t,$$(2)where $I_{ij}$ is the illumination intensity at location $(i,j)$ of the pattern displayed on a digital micromirror device (DMD), and $N_t$ is a white noise term describing the noise of the detector and environment. To reconstruct an image, $I_t$ is transformed to the spectral domain: $$S(\omega )={\mathcal{F}}_t\{I_t\}\propto M_{\mathrm{patt}}\sum _{i,j}{I}_{ij}\delta (f-f_{ij})+\overline{N},$$(3)where $\mathcal{F}$ denotes the FT, $M_{\mathrm{patt}}$ happens to be the number of data points in the discrete 1D time series (integration length) since the measurement is synchronized with the projection, and $\overline{N}$ is the average noise level. For each $f_{ij}$ in the patterns, we find the nearest frequency ${\widehat{f}}_{ij}$ in the spectrum and save its magnitude $S({\widehat{f}}_{ij})$. Finally, an $H\times W$ matrix is filled with all the $S({\widehat{f}}_{ij})$ in their locations $(i,j)$, and a heatmap of this matrix will yield an image of the object. It is the one-to-one correspondence between the frequency $f_{ij}$ and the spatial location $(i,j)$ that allows us to retrieve the spatial information computationally, and therefore we need to measure the transmitted light with only a single-pixel detector. We note here that the noise only contributes a constant term in the spectrum given by Eq. (3).

A schematic of the experiment setup is shown in Fig. 2(a). The grayscale patterns are first decomposed to $8M_{\mathrm{patt}}$ monochrome (1-bit, two pixel values) patterns to be compatible with the input format of the DMD. A solid-state laser (633 nm) with a power of $\sim 30\mathrm{mW}$ is used to illuminate the DMD, which spatially modulates the incident laser beam and generates a set of illumination projections with a spatial structure similar to the input monochrome patterns. Each projection has $20\times 60\mathrm{pixels}$, and each pixel is maintained by a $10\times 10$ array of DMD micromirrors (each micromirror has a size of $10.8\mathrm{\mu m}\times 10.8\mathrm{\mu m}$). We define such an array of mirrors as a single “DMD pixel.” The patterns are then imaged by a lens (L1) with a magnification of two. High-order images due to diffraction are filtered out by an iris (I) such that only the zeroth-order image with the strongest intensity is formed in the image plane. A high-contrast object (O) with letters “IQSE” ($3\mathrm{mm}\times 10\mathrm{mm}$) being transparent is placed at the image plane of L1, where the zeroth-order images of the patterns are directly projected on the region of “IQSE” without scattering media. Then, two ground glass diffusers (D1 and D2, 220 grit) are inserted to block the view of the object. The 1 mm thick ground glass diffusers with a 220 grit size are considered as standard highly scattering media [24,30]. The distances between the object and D1 and D2 are about 1 mm and 5 mm, respectively. The diffusers can be kept stationary or moved back and forth at a random speed by a motorized stage. A lens (L2) followed by a PD is placed behind D2 to collect the transmitted light. The system is synchronized. The PD will record the intensity each time a pattern is applied to the DMD.

The image of the object is first captured by a CMOS camera (replacing the PD) under three conditions: without scattering media (ND), with stationary diffusers (SD), and with dynamic diffusers (DD). As shown in Fig. 2(b), the object is invisible in the camera images when the diffusers are present (SD and DD), while the STEP imaging scheme can retrieve the images through both SD and DD [Fig. 2(c)]. We define $\beta =M_{\mathrm{patt}}/M_{\text{pixel}}$ as a scaled number of patterns, in which $M_{\mathrm{patt}}$ is the number of grayscale patterns, and $M_{\text{pixel}}=120$ is the number of pixels in one pattern. The measurements in Fig. 2(c) are performed with $\beta =8$, i.e., 9600 grayscale patterns ($8\times 9600$ monochrome patterns), and $\beta =64$. To reconstruct an image, the frequency resolution of the FT $\delta f=r_s/M_{\mathrm{patt}}$ must satisfy the condition $\delta f\ge \mathrm{\Delta }f$, which determines the minimum value of $\beta$. Our experimental parameters ($f_0=\mathrm{\Delta }f=0.1$, $H=20$, $W=60$) give $\beta \ge 8$. Images in Fig. 2(b) are obtained with $\beta =8$. Raw pixels in the reconstructed image are always the same as those in the patterns. Such a small number of pixels leads to pixelated images. Nevertheless, the pixelation effect can be eliminated computationally by applying bilinear interpolation to the raw pixels, and interpolated images of size $400\times 1200$ are given in the second column of Fig. 2(c) for each $\beta$.

The size of an individual pixel determines the imaging resolution in the patterns projected onto the object: the smaller the pixel is, the more details of the object can be resolved. The pixel size at the image plane of L1 is determined by the size of the DMD pixel and the magnification of the imaging system defined by L1 when there are no scattering media. Decreasing the DMD pixel size will increase the resolution for a given magnification. However, the optical power reflected by each DMD pixel will be reduced due to the shrinkage of the reflective area (fewer micromirrors), resulting in a decreased signal-to-noise ratio (SNR) of the measured light. Conversely, increasing the DMD pixel will provide better SNR at the expense of a lower resolution. If the scattering layers are inserted, the degradation of pattern quality and the attenuation of optical power due to scattering must also be considered. Therefore, the trade-off between resolution and SNR should be decided according to the configuration of a specific setup and the optical properties of the media. In our experiment, the resolution is $\sim 0.2\mathrm{mm}$, which is adequate for resolving the object of size $3\mathrm{mm}\times 10\mathrm{mm}$, and the transmitted light intensity is far above the shot noise level of the detector under the scattering of 220 grit ground glass diffusers and 1.2 mm thick chicken breast slices.

To further demonstrate that STEP is also insensitive to the motion of the scattering centers in the media, we replaced the ground glass diffusers in Fig. 2(a) with two slices of chicken breast and performed similar measurements with $\beta =96$. In general, chicken breast has the following optical properties for light at a wavelength of 633 nm [18]: scattering coefficient 23 mm⁻¹, reduced scattering coefficient 0.8 mm⁻¹, and absorption coefficient 0.01 mm⁻¹. These parameters imply that scattering is the dominant process that prevents the tissue from being transparent. As shown in Fig. 3(b), the object cannot be resolved in the camera image, whereas the image can be reconstructed by STEP [Fig. 3(c)]. The larger the $\beta$, the better the spectrum’s ability to distinguish target frequencies from background noise. A much larger $\beta$ to perform the computation implies stronger tissue scattering than in the ground glass diffusers. In the reconstructed image, some bright spots appear due to the uneven texture of the tissue, which results in less evenly distributed transmitted light intensity than for the ground glass diffusers. Nevertheless, the results suggest that STEP can image through dynamic scattering media.

We define the visibility of the reconstructed image to be $v=({\overline{p}}_s-{\overline{p}}_b)/({\overline{p}}_s+{\overline{p}}_b)$, in which ${\overline{p}}_s$ and ${\overline{p}}_b$ are the average pixel values of the signal (“IQSE” regions) and background (other regions), respectively. Different from conventional structured illumination techniques in which the information required for image reconstruction is encoded only in the spatial structure of the patterns [31–35], the STEP illumination also encodes the information (frequencies) in time. One advantage of introducing the time-domain encoding is that it allows further improvement of the visibility via halving a successive segment of the data and calculating the cross-spectrum (CS) of the two halves. The cross-correlation technique has been demonstrated very effective in weak signal detection [36–38]. For images reconstructed with a specific $\beta$, higher visibility is obtained if the searching of target frequencies is performed in the CS of the two halves of the data [Fig. 4(b)] rather than in the Fourier spectrum of all the data [Fig. 4(a)]. As shown in Fig. 4(b), we divide the whole data set into two segments of the same length, and their cross-correlation is found to be $$I_1\star I_2=\frac{M}{2}{I}_{ij}{I}_{ij}^{\prime }\sin \left(2\pi f_{ij}\frac{\tau }{r_s}\right)+\frac{M}{2}{\overline{N}}_1{\overline{N}}_2,$$(4)where $M$ is the total length of the two signal segments, and ${\overline{N}}_1$ and ${\overline{N}}_2$ are the average noise levels for the two segments. It follows that the CS is $$S_{12}(f)={\mathcal{F}}_{\tau }\{I_1\star I_2\}\propto \frac{M}{2}\sum _{ij}{I}_{ij}{I}_{ij}^{\prime }\delta (f-f_{ij})+\frac{M}{2}\delta (0).$$(5)As we have seen in Eq. (3), the noise contributes a constant background equally at every frequency in the spectrum. However, in the CS [Eq. (5)], the noise is concentrated to zero frequency [indicated by $\delta (0)$] and never contributes to the cross-spectral magnitude at any other frequencies. Therefore, the SNR is enhanced, leading to sharper peaks of the target frequencies and less noisy reconstructed images.

The image reconstruction can also be implemented by calculating the correlations between the measured intensity data $I_t$ and the time series $a_{ijt}$ in the original patterns [31,32,35]: $$C_{ij}=\sum _{t=0}^{M-1}{I}_t{a}_{ijt}.$$(6)Filling an $H\times W$ matrix with all the $C_{ij}$ at their locations $(i,j)$ will yield an image of the object. However, the visibility of the image reconstructed by correlation is lower than that of the image produced by FFT-based methods, because in this case the correlation is done without a shift in time, leading to a non-trivial contribution of white noise. Figure 4(d) compares the visibility of the images generated by the three algorithms for different $\beta$ in the cases of SD and DD. It is worth mentioning that the performances of STEP-FT and STEP-CS are almost independent of $\beta$: they have similar visibility over the investigated range of $\beta$ and show saturation behaviors. This means high-quality images may be obtained with small data sets, thus consuming less time on measurement and computation. On the other hand, the correlation has the worst overall performance and is sensitive to the value of $\beta$. Therefore, high-quality images may be acquired only with a large number of data points.

The correlation method also suffers from low computational efficiency. For $M$ patterns of size $H\times W$, both the time and space complexity of the image reconstruction via correlation are $O(HWM)$, in terms of the big-O notation, whereas the time complexity of FFT-based reconstruction (FT and CS) is $O(M{\log }_{2}M)$ since the FFT algorithm is used to compute the spectrum, and the space complexity is $O(M)$, as there is no need to store the original patterns. A benchmark of the computing time for FFT-based and correlation algorithms is given in Fig. 4(e), where the rapid increase in computing time for the correlation algorithm provides a sharp contrast with those of FFT-based methods. The high computational efficiency of the FFT-based algorithm may enable fast image processing with devices having limited computing resources.

In conclusion, we have developed a computational imaging method named “STEP” to realize non-invasive imaging through scattering media with a single-pixel PD. We have demonstrated the ability of imaging an object sandwiched between two opaque scattering glass diffusers, as well as for the same object embedded in chicken breast layers, where the light experienced multiple scattering. This method is insensitive to the motion of the scattering centers in the media. The design of STEP removes the requirement of a high-resolution camera. It allows an elegant FFT-based image reconstruction algorithm that is more computationally efficient than correlation-based methods, which may be more favorable in many application fields. Our technique provides new perspectives for peeking through multiple scattering turbid media and enables potential fast images to diagnose biological tissues.