Scattered light imaging beyond the memory effect using the dynamic properties of thick turbid media

Abstract. Scattered light imaging through complex turbid media has significant applications in biomedical and optical research. For the past decade, various approaches have been proposed for rapidly reconstructing full-color, depth-extended images by introducing point spread functions (PSFs). However, because most of these methods consider memory effects (MEs), the PSFs have angular shift invariance over certain ranges of angles. This assumption is valid for only thin turbid media and hinders broader applications of these technologies in thick media. Furthermore, the time-variant characteristics of scattering media determine that the PSF acquisition and image reconstruction times must be less than the speckle decorrelation time, which is usually difficult to achieve. We demonstrate that image reconstruction methods can be applied to time-variant thick turbid media. Using the time-variant characteristics, the PSFs in dynamic turbid media within certain time intervals are recorded, and ergodic scattering regimes are achieved and combined as ensemble point spread functions (ePSFs). The ePSF traverses shift-invariant regions in the turbid media and retrieves objects beyond the ME. Furthermore, our theory and experimental results verify that our approach is applicable to thick turbid media with thickness of 1 cm at visible incident wavelengths.


Introduction
Improving the quality of optical imaging is a primary goal in optical research. However, most objects and their surroundings, such as biological tissues and smog, are optically complex and turbid, inducing strong scattering light and drastically degrading image quality. Thus imaging through opaque diffusive media has become an important technical challenge, and many methods have been proposed to address this fundamental problem. Straightforward strategies that use only unscattered, "ballistic" photons selected by optical coherence tomography, 1 two-photon microscopy, 2 or holographic gating 3 have proven to be useful but not efficient because the number of ballistic photons that can be used is intrinsically limited. Wavefront shaping 4-8 through phase compensation with a spatial light modulator is an alternative solution that allows focusing and imaging through stronger scattering samples. Furthermore, a variety of improved approaches have been proposed with the initial assistance from guide stars 9,10 and known objects. 11,12 Matrix factorization combined with linear fluorescence has shown its capability to reconstruct the transmission matrices in biological microscopy. 13,14 The constructive two-photon interference formed by the thin dynamic scattering layer also enables direct imaging of the object. 15 With the development of machine learning, neural networkbased imaging methods have demonstrated their ability to handle dynamic scattering and additional noise. 16,17 A breakthrough approach that exploits inherent angular correlations in scattered speckle patterns [18][19][20] was used to successfully eliminate prior calibration with reference sources; however, the sequential scanning illumination is applicable only to stationary objects. 21 In addition, the intensity autocorrelation between the object image and the object can be reconstructed using a phase retrieval algorithm. [22][23][24] Optical imaging using the Fourier-domain shower-curtain effect has also yielded results in dynamic turbid media. 25 The core concept of these approaches is to consider the scattering medium as a "lens" so that the medium can be described by a point spread function (PSF). Based on this technique, high-speed, full-color images in scattering media can be reconstructed over large ranges. 26 Moreover, the PSF of a scattering medium can be extended to a depth that can be used for three-dimensional (3D) imaging. 27,28 The DiffuserCam, which uses a well-designed compress sensing algorithm, was proposed for 3D light-field imaging, and applies a single diffuser rather than a lens array. 29 Important results were also achieved using the noninvasive optical imaging by speckle ensemble technique, 30,31 demonstrating that the ensemble average of the PSF can be separated into a normal imaging part and a variance part for weak scatters, allowing targets to be reconstructed by microlens arrays or rotating scatters. However, the PSF is angular shift invariant over a very limited range of angles. 19,20 This phenomenon is known as the memory effect (ME). An increase in the thickness of the scattering medium leads to a limited expansion of the ME range, which hinders the use of these techniques in wider applications. Imaging through diffusers with larger ranges has attracted increasing attention, and many breakthroughs have been achieved through various algorithms and optical designs. [32][33][34][35] Therefore, the development of a general approach that is applicable to thick turbid media is essential for imaging scenarios that are currently inaccessible.
Based on the above description, we found that the PSF concept can be modified and extended for imaging through a thick turbid medium as long as the medium is time variant. Throughout the slowly changing scattering medium, a series of PSFs can be obtained and combined as an ensemble point spread function (ePSF) to achieve ergodicity. 36 The ePSF traverses the shift-invariant regions of the turbid medium and can retrieve objects beyond the ME range. In our experiment, we successfully applied this approach and reconstructed objects in a 1-cm-thick container filled with turbid liquid and dispersed calcium carbonate (CaCO 3 ) particles at visible incident wavelengths. The reconstruction range approaches that of a scattering system and is beyond the ME range. Compared with conventional single-exposure image recovery methods in static scattering media, the dynamic media allow a larger reconstruction range and a longer lifetime for the measured scattering properties. The dual advantages in time and space introduce the possibility of developing more applications in fields ranging from biology to medicine.

Principle
A PSF is an impulse response of a focused optical system. The intensity distribution of the light field under incoherent illumination is described by a convolution, I ¼ O Ã PSF, where I and O are the intensity distributions on the image and object planes, respectively, and the symbol * denotes the convolution operation. Several studies have proposed that thin turbid media can be considered lenses, 37,38 with the scattering properties determined by the PSF: where ðx; yÞ are the spatial coordinates and the subscripts i and o denote the image and object planes, respectively. In this case, I and the PSF are both complicated speckle patterns, and Eq. (1) is valid only within a shift-invariant region known as the ME range. By regarding the imaging system as a "black box," object O can be restored through a deconvolution process with the PSF measured by a reference point. If the turbid property changes or the object O moves beyond the ME range, the PSF and image I must change simultaneously to ensure that Eq. (1) remains valid. For a time-variant scattering medium, Eq. (1) is valid only when the exposure time is shorter than the decorrelation time of the scattering patterns. If the exposure time is longer than the decorrelation time, the captured PSF and image I are the summation of different scattering states. Each scattering state describes a certain scattering channel for a group of incident photons. The concept of summing different scattering states is also applicable to multiple capturing applications. In these cases, Eq. (1) is rewritten as where n is the serial number of each pair of I and PSF. As the time varying process is physically continuous, the summation notation in Eq. (2) can be moved inside the integral and rewritten as where ePSF ¼ P n PSF n and eI ¼ P n I n are called the ePSF and ensemble image, respectively. In contrast to the normal PSF, the ePSF is superposed by each PSF n that traverses different shiftinvariant regions in the turbid medium.
Thus the ePSF describes almost all the scattering properties and optical information of the medium for the closed scattering system, which contains finite ergodic states, resulting such that the objects can be retrieved outside the ME range of any single scattering state. Higher summation times yield more complete ensemble scattering states that can be recorded by ePSF and eI; therefore, the retrieved object information is closer to the original information. Figure 1 illustrates this physical concept.

Proof of Concept
The ME is the first order of the correlation function when the estimated transport mean free path within the scattering medium is substantially larger than the wavelength and is defined as 18,20 CðqLÞ ¼ ½qL∕ sinhðqLÞ 2 ; (4) where L is the thickness of the scattering medium, q ¼ 2πθ∕λ, θ is the scattering angle of the incident light, and λ is the incident wavelength. The ME can then be determined according to the full width at half-maximum of the correlation curve. To prove the concept of scattering state traversal, we devised an experiment to investigate a thin scattering medium. As shown in Fig. 2(a), a physical pinhole (∅ ¼ 200 μm) illuminated by an incoherent light source (528 AE 10 nm, Daheng Optics GCI-060403) was placed on the optical axis of the imaging system. The diffractive light was projected onto a standard scattering medium (Newport 5-deg circular light shaping diffuser) and diffused. The object plane was set at a distance of d o ¼ 10 mm from the diffuser. In this experiment, the diffuser was suspended and swung as a pendulum. The string was l ¼ 300 mm long, and the mechanical frequency was 0.91 Hz according to a simple calculation. It is worth noting that the mechanical frequency was intentionally set at a low value to ensure that the captured images were blurred. A higher frequency may weaken the diffusion property, leading to a scenario in which the "scattering photon" no longer dominates. 39 A monochromatic CMOS (Basler acA2040-90um) equipped with a main lens (Zeiss Milvus 2/100M ZF.2) was used to capture a series of similar PSF n s at a shutter speed of 24 fps by placing a pinhole in the object plane. Then the pinhole was replaced by a 1.5-mm-long letter "G" as an object, and a series of images (I n s) were captured at the same speed. The letter was hollow, while the outside was completely opaque.
Figures 2(b) and 2(c) show the ePSF and eI at different superposition times N. Because the PSF n s and I n s are captured at different times, the object cannot be reconstructed for N ¼ 1, as shown in the first image in Fig. 2(d). However, as the number of superpositions increases, the reconstructed object becomes clearer, and the object is well distinguished when the number of superpositions reaches N ¼ 500, as shown in the last image in Fig. 2(d). This finding demonstrates that the ePSF is effective for reconstructing unknown objects through time-variant scattering media. It is worth noting that the ballistic light is not significant when the number of superpositions exceeds 500, demonstrating that the target light is sufficiently scattered and recovered by the scattered light in a dominant manner.

Imaging Through Thick Turbid Media
According to Eq. (4), the ME range becomes narrower as the scattering medium becomes thicker. For example, when a typical static 1-cm-thick scattering medium is illuminated by incident light with λ ¼ 532 nm, the ME angle is in the range of 0.1 mrad, 38,40 resulting in an extremely narrow ME range for practical applications. To prove that our proposed method can retrieve targets outside the ME range in thick scattering media, the thin diffuser is replaced by a 1-cm-thick cuvette (inner dimension) containing a CaCO 3 turbid suspension.
The optical attenuation length can effectively characterize the scattering property of the media, and it can be obtained by measuring the transmittance of ballistic light. 17 As the demonstration experiments use a 528-nm LED, a laser light source with similar wavelength (532 nm) is used for the measurement. An optical power detector (Newport 818-ST2/DB) is used to measure the output light intensity with and without the CaCO 3 turbid suspension. A ∅ ¼ 1 mm iris was used to screen the ballistic light during the measurement. According to the measurement, the optical attenuation length of the suspension for the following experiment reaches 13.3, corresponding to a transport mean free path l Ã ¼ 0.75 mm.
The experimental setup is shown in Fig. 3(a). The object plane was set at a distance of d o ¼ 12 cm from the right surface of the cuvette. The illumination was provided by the same LED and shaped by a ∅ ¼ 100 mm, f ¼ 200 mm convex lens to generate a uniform and sufficiently large spot. Before the experiment, the cuvette was treated with an ultrasonic cleaner to remove any CaCO 3 particles attached to the surface to prevent the formation of a static scattering layer. To prevent precipitates of the particles and ensure that the scattering property changed over time, the liquid was slowly agitated with a dropper purge. To obtain individual single states (or ensembles of partial states), the state-to-state differences must be within the numerical precision limit of the system. Thus to obtain speckle data with high numerical precision, a pco.edge 4.2 sCMOS with a dynamic range of 37,500:1 was used to capture the images. A Thorlabs MVL7000 lens was used to collect scattered light.
The reconstruction range of the scattering system is measured by calculating the normalized intensity of the recovered image when translating a round hole in the horizontal direction. 19,26,38 The same procedure is followed to evaluate the reconstruction range that can be achieved by our recovery method. The irradiance of an optical system decreases as the off-axis angle increases; 41 accordingly, the yellow curve in Fig. 3(b) shows this vignetting characteristic in the experimental setup, and the whole reconstruction range is measured to be AE205 mrad. Then a thick turbid medium is introduced, and the round hole image is reconstructed at various angles by our proposed method (the number of superpositions is 1000). The blue dots in Fig. 3(b) illustrate the normalized intensity of each reconstructed image. Figure 3(e) demonstrates that our method can consistently reconstruct images at different angles. The reconstruction range is close to the range of the original scattering system, indicating that the ePSF is effective for addressing ME limitations. The scattered results in Figs. 2 and 3 indicate that the experimental system experienced strong scattering, as the characteristics of our results are significantly different from those of the ballistic light-based imaging process.
To directly visualize the reconstruction range beyond the ME range and evaluate the effect of our method with different objects, a series of letters are simultaneously placed in the object plane. The speckles in the image plane are captured with the above procedure. The letters and their reconstructed images are illustrated in Fig. 4, demonstrating the powerful ability of the ePSF to image beyond the ME range. According to Fig. 4(b), as the number of superpositions increases, the traversal of the scattering state is gradually completed. eI and ePSF both achieve ergodic scattering processes in a specific medium, consistently allowing the scattered target to be revealed at different positions in the reconstruction range. The quality of the reconstructed image is degraded compared to the raw image due to the aberrations from the imaging system and the signal-to-noise ratio of the speckle.

Discussion
It is worth noting that an ensemble of speckles is necessary for achieving ergodicity in reconstructed images in thick turbid media. Theoretically, each scattering state can be collected separately to determine the ensemble average. When the motion of a thin scatterer is regarded as a quasistatic process, the correlations among the speckles decrease with each small movement d, satisfying d > r ME (r ME corresponds to the radius of the ME range, which is the same for every possible state, and d is described by the shifting velocity v and characteristic time τ as d ¼ vτ). The extraction of the speckles in each state requires fine control over the movement speed and exposure parameters. With specific coordination, it is possible to achieve image reconstruction for each state. 42,43 Moreover, shifting the diffusive media is an effective method for destroying the speckle structure and suppressing speckles. [44][45][46] For time-variant thick turbid media, the state changes caused by the overall macroscopic movement are transformed into various state changes caused by countless microscopic movements. Since the decorrelation time, which measures the state transition speed, is typically much smaller than the exposure time (from 10 to 200 ms in our experiments), 47 a single image contains an ensemble of multiple speckles in various states. A single image of speckles with a sufficiently long exposure time is theoretically equivalent to an ensemble of speckles after a series of images acquired over the same number of states. Both methods allow us to capture information about a large number of scattering states. However, imaging sensors cannot usually achieve a enough long exposure time due to insufficient dynamic ranges, leading to overexposure and a failure to reconstruct the target through thick turbid media, while an ensemble of speckles is effective at preserving as much scattering information about the object as possible. In addition, the proposed ergodic method is only applicable to the systems that are closed within the acquisition time. In the closed systems, the total ergodic states are finite, and the exact same sequence of PSFs can be traversed for several times to be matched with the measurement. However, for the open systems where the total ergodic states are infinite, the same sequence of scattering states cannot be traversed for another time, resulting in the failure of using the proposed method.
In practical application scenarios, objects originally existing in the background are usually known, such as bones, blood vessels, or fluorescent structures. Imaging targets are often emerging objects around known objects, such as abnormal biological tissue structures. The known objects O R and the emerging targets O T could form their corresponding ensemble speckle patterns eI R and eI T through the scattering media after reaching ergodicity: As the convolution in the spatial domain is a multiplication in the spatial frequency domain, imaging targets can be recovered by the processing of the Fourier transform: where "F" represents the Fourier transform, and × indicates multiplication. Therefore, the a priori information can be used to image unknown targets over a large area. 48 Using the same experimental setup in Fig. 3(a), as shown in Fig. 5, we know the letter "H" as a prior information. After the collection of its ensemble speckle, two unknown letters are  inserted and form their corresponding ensemble speckles. With the process proposed in Eqs. (5) and (6), the unknown letters can be reconstructed to be "G" and "F."

Conclusions
We demonstrate that the PSF concept can be modified and extended to image through thick turbid media as long as the media are time variant. In this paper, the change in scattering media over time is used to reach ergodicity and obtain the ePSF, which allows the reconstruction range to exceed ME limits. The concept of the ePSF was developed to explain scattering state traversal and is capable of retrieving objects beyond the ME range in time-variant thick turbid media. We successfully applied this approach in experiments to reconstruct objects through a 1-cm-thick container filled with a stirring turbid liquid dispersion of CaCO 3 particles at visible incident wavelengths. The reconstruction range reaches the spatial limit of the scattering system and exceeds the ME range. Compared with static scattering media, time-variant scattering media exhibit temporal and spatial advantages in image recovery. We believe that this approach has great potential in various applications, such as biomedical imaging under strong diffusive circumstances.