Conventionally, an optical imaging system can be regarded simply as a one-to-one mapping system, as a spherical wavelet emitted from a local point on the object plane is converged by the imaging optic and forms a unique point (subject to diffraction) on the imaging plane. In this way, a clear image of the object can be formed. However, the presence of scattering media between the object and the imaging optic prohibits such a clear image from being formed because a significant part of the light propagating all the way from the object to the optic is scattered many times, producing a noise-like scattered pattern on the image plane. As the optical thickness of the scattering medium is larger than a certain value, the scattered light arriving at the camera becomes more intensive than the unscattered or ballistic one; further, the image is submerged by the noise-like scattered pattern. Thus, an intuitive method for image acquisition under such hazed conditions is to select the ballistic light by using gating techniques in the time [1], space [2], polarization [3], or even coherence [4] domains.

Actually, the propagation and transport of coherent waves in disordered media has been one of the central problems in many different disciplines of science and engineering, ranging from microwaves to electron waves. It has been shown that coherent scattering in a static scattering medium is linear and deterministic and can be described by a transmission matrix [5]. One of the lines of recent studies then focuses on the ways of reversing the deterministic coherent scattering process. Physically, this can be done by generating a phase conjugated version of the scattered wave [6,7], so that it can trace the way back through the scattering medium in the first place and form a clear image. An alternative way is to measure the transmission matrix [8] or precompensate the phase distortion introduced by the scattering process by the technique of wavefront shaping [5]. When applying these methods in dynamic scattering media, the wavefront measurement and playback must be made within the decorrelation time [9]. Therefore, faster wavefront modulators, parallel processing, and more effective optimization algorithms are highly demanded in these cases [10,11].

Alternatively, one can utilize some intrinsic statistical properties of the scattered light instead of reversing the scattering process. For example, the short-range second-order correlation of a speckle pattern formed by the scattered light, or the memory effect [12,13], has been found to be particularly interesting in this case, as it allows computational reconstruction of an object image from the autocorrelation of the scattered light intensity by using a regular phase-retrieval algorithm [14,15]. However, the memory effect offers a field of view (FoV) that is inversely proportional to the effective thickness of the medium. Thus, this method is usually applicable for imaging through optically thin scattering media [14,15] or around the corner [16]. We note that efforts have been made to enlarge the FoV by using multiplexing or scanning [17,18]. But these techniques have limited effect in the case when the scattering medium is not static, as the dynamic characteristics of the media and the spectral width of the light in many practical environments will further shrink the correlation length of the scattered light [15,19] and thereby the FoV of the imaging system.

An interesting and practically important question then naturally arises, i.e., if it is possible to reconstruct an image when the coherence of the scattered light is completely lost. At least two factors may lead to the loss of coherence: multiple scattering of the light during its propagation inside the scattering medium [19] and the coherence nature of the light itself.

To answer this question, we attempt to demonstrate that it is still possible to reconstruct the image of an object hidden behind a nonstatic and optically thick diffusive medium and under the illumination of incoherent light. To mimic a realistic environment, we also switch on an ambient light so as to introduce a strong interference noise that is incoherently superimposed onto the scattered light from the object. The method we use to address this problem is deep neural networks.

Deep neural networks have shown great potential in solving many computational optical imaging problems [20], ranging from digital holography reconstruction [21,22], phase imaging [23,24], computational ghost or single-pixel imaging [25,26], to coherent imaging through optically thin [27,28] and even thick scattering media [29,30]. As almost all the deep neural networks for computational optical imaging were trained supervisedly (except, for example, the work reported in Ref. [24]), it is therefore required that each pair of labeled images in the training set should have exclusively unique features [31]. In the application of imaging through scattering media, when the illumination is coherent [27–30], what is recorded by the camera is a high-contrast speckle pattern, the morphology of which is explicitly object-dependent. However, when the objects are illuminated with an incoherent light, as in our case, the coherent effect will be smeared out, yielding milky patterns that are almost visually indistinguishable with respect to one another. Therefore, additional care has to be taken to identify unique features to distinguish each milky pattern.

This article is organized as follows. The problem of incoherent imaging through highly dynamic turbid media is formally described in Section 2. Measurements of the parameters that characterize the turbid medium are presented in Section 3. The experimental demonstration of the proposed method is presented in Section 4.

The optical thickness of the tank of fat emulsion suspensions was about 16 (refer to Section 3 for the details of calibration). As the absorption of the intralipid is low, we can conclude that it is the scattering that mostly accounts for the attenuation of the light. Then, the object’s scattered light intensity S{Io} can be expressed as S{Io}=SW{Io}+SL{Io},(1)where SW{Io} and SL{Io} are the intensity of weak light and large-angle scattered light, respectively, without considering the ambient light.

In our experimental demonstration, we also take the ambient scattered light Ia into account. Thus, the final intensity pattern captured by the camera can be written as Is=SW{Io}+SL{Io}+Ia,(2)where the interference noise In=SL{Io}+Ia, and SW can be approximately regarded as the signal component of the recorded image Is. Thus, the signal-to-interference ratio SIR=10log10SW{Io}/In=−17dB.

Thus, the purpose is to reconstruct the object image Io from Is.

Fat emulsion is similar to milk in many ways and mainly consists of soybean oil, water, glycerin, and eggphospholipid. Due to emulsification in the presence of phospholipids, the oil is suspended in water as small droplets with a momolayer lipid membrane. As a large amount of fat droplets are randomly distributed in the emulsion suspensions, a beam of light that propagates through it is distorted owing to multiple scattering, making it difficult to predict the propagation trajectory for each light ray.

Owing to the advantages such as low absorption, easily adjustable scattering coefficient by dilution, standardized particle size distribution, nontoxicity, and low cost [32], intralipid has been frequently used as an optical phantom that mimics a turbid medium, in particular, for biomedical applications. For example, it has been widely used in the calibration of clinical application systems [33,34] and as a diffusive reference standard for optical tissue modeling [35,36]. Many studies have been carried out to measure the scattering properties of intralipid suspensions [32,37]. However, when it is used as a scattering medium in imaging, the decorrelation time and optical thickness are quite important parameters related to the scattering, which are rarely measured or quantified. Thus, we need to calibrate these two parameters in this section.

The optical thickness (OT) is a parameter deduced from the Lambert–Beer law: I=I0e−OT=I0·e−μ·c·L,(3)where I0 and I are the ballistic light intensities measured in front of and behind the tank of intralipid suspensions, μ is the attenuation coefficient, c the concentration of intralipid dilution, and L the path length through the suspensions. The attenuation coefficient μ is the sum of scattering coefficient μs and the absorption coefficient μa, i.e., μ=μs+μa.

Since we used the 530 nm LED (M530L3, Thorlabs) in our demonstrations, the optical thickness of the fat emulsion should be measured at or at least closed to this wavelength. We thus used a semiconductor pumped green laser with a wavelength of 532 nm as the light source in the calibration process. An optical power meter (PM200, Thorlabs) with the probe diameter of 9.5 mm was used to measure I0 and I. Then, the optical thickness can be calculated according to Eq. (3). Note that the contributions to the absorption of fat emulsion are mainly from water and soybean oil, while the influence of glycerin and egg lipid can be ignored due to the low concentration. For the light at around 530 nm, the absorption coefficient μa is about 10−4mm−1 according to a previous study [32], so that its influence can be ignored; therefore, μ≈μs.

The averaged diameter of the droplet is around 0.3 μm [38]; at a scale comparable with the wavelength of the light, Mie scattering occurs when light passes through the fat emulsion suspension. During our experiment, we observed the Tyndall effect from the side of the tank, as shown later in this paper, which was taken by a mobile phone camera at the experiment site. Every droplet in the fat emulsion can be seen as a secondary source that re-emits a wavelet, the superposition of which forms a speckle pattern. However, we observed an interesting phenomenon, i.e., the appearance of the speckle pattern changes with respect to the density of the intralipid suspensions. When the density is low, the speckle grain has a shape like a thin stripe [Fig. 2(b)–2(e)]. The reasons for this can be addressed as follows. First, the density of fat emulsion is low so that only a few scattering events occur to the light that propagates through it. That is why the contrast of the speckle in this case is very high compared with Figs. 2(i) and 2(j). Second, the droplets are moving downward owing to the gravity of the Earth during the course of data acquisition. Finally, the camera plane may be tilted with respect to the incoming light [39]. As the intralipid density increases, many more scattering events occur as the light propagates through the tank. As a result, the speckle pattern is formed by the interference of light waves that experienced various times of scattering. The speckle grain becomes more isotropic, but the contrast is apparently reduced (to 0.08).

As noted, a large number of fat droplets are randomly distributed in emulsion suspension. Owing to gravity and Brownian motion in the surrounding liquid [40], these droplet scatterers are not static but moving randomly, resulting in fluctuations of light propagation trajectory in the turbid suspensions. For analyzing the influence of particle motion on the light scattering in this case, one usually attempts to measure the temporal correlation by using diffusing wave spectroscopy (DWS) [40–42]. In this paper, we propose an alternative method to characterize the temporal correlation property of the intralipid diluents with different concentrations in the 33.6 cm thick glass tank.

In order to determine the decorrelation time, we need to fit the data of C(Δτ) to a theoretical model. Note that there are various contributions to the correlation function that are determined by the different sets of multiple scattering trajectories, which are classified in terms of the crossing number of two diffusive paths [44]. Each contribution has discrepant time-dependence owing to different pairings of two trajectories obtained at time t0 and t0+Δτ, respectively. One dominant contribution, C(1), the one without crossings, decreases exponentially with Δτ. Further, it specifies the memory effect [12,13]. Another considerable contribution, C(2), has one crossing, which decreases algebraically [45]. Generally speaking, for one more path crossing, the corresponding contribution is reduced by a factor g, which is a large dimensionless number with a value at the order of 102 for a visible light in liquid suspensions [44,45]. As the contributions of other correlation terms with multiple crossings can be neglected, we only consider the effect of C(1) and C(2) here. Specifically, these two terms can be expressed as [19] C(1)(Δτ)=[L/Lssinh(L/Ls)]2,(6)C(2)(Δτ)=1g1sinh2(L/Ls)[sinh(2L/Ls)2L/Ls−1],(7)where L is the geometric thickness of sample, and Ls=Dτe·f(Δτ),(8)where D is the diffusion coefficient that describes wave diffusion and is related to the elastic mean free path le, τe is the elastic collision time, and f(Δτ)=[e−Δτ/(2τb)1−e−Δτ/(2τb)]1/2,(9)where τb describes the Brownian motion of the scatterers. One should note that τb is the characteristic time for a scatterer to move a distance in the order of wavelength λ.

We fit the experimental data with the weighted sum of C(1) and C(2) defined in Eqs. (6) and (7), respectively, with the fitting parameters a, b, m, and n as $$\begin{gathered} C(Δτ)=a{m/f(n·Δτ)sinh[m/f(n·Δτ)]}2 \\ +b1{sinh[m/f(n·Δτ)]}2{sinh[2m/f(n·Δτ)]2m/f(n·Δτ)−1}. \end{gathered}$$(10)Apparently, the parameters a and b represent the related contributions of C(1) and C(2) to C. And m is a parameter related to L/Dτe, while n is associated with τb. The fitting results are plotted as solid lines in Fig. 5. From these fitted curves, one can easily determine the decorrelation time τd, which is defined as the time where the correlation coefficient drops to 1/e of its maximum. We found that τd for the suspensions with 0.8, 1.2, 1.6, and 2.4 mL intralipid 20% are 13, 4.8, 2.6, and 1.2 ms, respectively. Unfortunately, owing to the light being scattered too many times when we applied 3.6 mL intralipid or more in the tank, the decorrelation time drops dramatically, and the cameras are not fast enough to capture this process. We can conclude that the decorrelation time in this cases must be less than 100 μs, which is the shutter speed limit of our cameras. Therefore, we assume that the decorrelation time of the turbid suspensions is around several to several tens of microseconds, corresponding to an optical depth of around 16.

For the demonstration of the proposed method, we developed an experimental system, schematically shown in Fig. 1, to acquire all the data. In the first experiment, the optical thickness of the fat emulsion suspensions was 8 (2 mL intralipid 20%). In the second set of experiments, the optical thickness of the suspensions was about 16 (4 mL intralipid 20%), and decorrelation time was at the order of microseconds.

Here, we propose an end-to-end deep neural net to reconstruct object images from the acquired data. For convenience, we used 10,050 images altogether from the MNIST data set [46] as the training and test data. We first resized the original MNIST images to 512×512 pixels and displayed them sequentially onto the central pixels of the 1024×1080 pixel SLM. An incoherent LED light was used to illuminate the system. Reflected from the SLM, the incoherent light was then projected onto the front surface of the tank by using a delayed optic. As the light transmitted from the front surface all the way through the turbid suspensions to the back surface of the tank, it encountered many scattering events, as discussed in Section 3. The scattered patterns at the back surface of the tank were then captured by the sCMOS camera with an imaging optic. Note that there was also strong interference noise produced by an ambient light, with SIR=−17dB, as described by Eq. (2). Only the central 512×512 out of the 1024×1024 pixels of the raw image taken by the camera were used in the reconstruction process. Our previous study [29] shows that this does not affect the quality of the reconstructed image. Owing to the limited frame rate of the SLM and the sCMOS camera, it took about 3.5 h to capture all 10,050 scattered patterns. Then, we partitioned them exclusively into two groups, one with 9900 patterns as the training set after pairing up with their corresponding ground true images (those displayed on the SLM), and the other as the test set. Typical visible features of the scattered patterns are shown in Fig. 1.

The objective, then, is to reconstruct the object Io from the corresponding recorded scattered pattern Is. For the supervised training of a neural network, we need to construct a training set by pairing up every MNIST image in the training set together with the corresponding scattering pattern Is, i.e., Straining={(Is(1),Io(1)),(Is(2),Io(2)),…,(Is(N),Io(N))}. Then, we can converse the problem to learning the probability distribution P(Io|Is) [47] or the mapping relation Rlearn: Is→Io from the training set. Once it is learned, an object Io in the test set should be predicted from the corresponding scattered pattern Is [29]. The objective function then can be defined as Rlearn=argminRθ,θ∈Θ∑n=1NL(Io(n),Rθ{Is(n)})+g(θ),(11)where the set Θ contains handcraft parameters that specify the network structure and weighting parameters that are automatically adjusted during the training process [48], L(·) is the loss function that evaluates the error between label Ion and predicted output Rθ{Isn}, and g(θ) is a regularizer defined on the parameters with the aim of avoiding overfitting [20].

We adopted the convolutional neural network (CNN) [46] to solve the problem described by Eq. (11). Refer to Appendix A for more details about the structure of the neural net we designed.

To train the neural net, we must specify the loss function in Eq. (11). In this work, we adopted the mean square error (MSE) as the loss function to evaluate the predicted network outputs and ground truth images: MSE=min1WHN1∑n=1N1∑(u,v)W,H[Ip(n)(u,v)−Io(n)(u,v)]2,(12)where Ip(n)(u,v) is defined as the reconstructed image from the nth scattered pattern Is(n)(u,v), Io(n)(u,v) is the corresponding ground-truth resized MNIST image, W and H are the width and height of the reconstructed image, and N1 is the batch size, which was set to 5 in this work. We used the stochastic gradient descent method [49] to evaluate the MSE and Adam [50] to adjust the weights. The training step number was set as 30,000; in practice, however, we were able to use less training steps.

As mentioned above, we also performed a second set of experiments, in which the optical thickness of the turbid suspension is 16. In this case, we muse block out a significant part of the large-angle scattered light by placing two apertures between the back surface of the tank and the camera. The captured scattered patterns then look dimmed, as shown in panel (2), Fig. 1. The corresponding reconstructed images from the proposed neural net are plotted in Fig. 6(c). With respect to the ground truths, the RMSE was about 4.66 and the SSIM was around 0.98, at the same level as those shown in Fig. 6(b).

In this subsection, we concisely analyze the robustness of the proposed technique. The optical thickness of the turbid suspension is 16 in all the analysis here.

The robustness to the scaling and rotation of the object/camera is shown in Fig. 8. In both cases, the image can be reconstructed well. Again, it is more robust to the change in the object (both scaling and rotation) owing to the reason mentioned above.

In conclusion, we have demonstrated an end-to-end deep-learning-based incoherent imaging method through optically thick and dynamic scattering media under a condition of strong interference noise (SIR=−17dB). The proposed method allows direct reconstruction of an object image from the captured scattered pattern with high quality (with an averaged RMSE around 5 and averaged SSIM about 0.86). We have quantitatively measured the optical thickness and decorrelation time of the turbid suspension [which was made by mixing fat emulsion (intralipid 20%) with purified water], with respect to the intralipid concentrations, and calibrated the optical thickness of the tank of intralipid suspensions in our experiments, which is about 16 in the severest case. Although the results were obtained in the lab environment, it is expected that the proposed method can be used in a wider scope of other scattering environments.

Acknowledgment. We thank the anonymous reviewers for their constructive comments.