Light field (LF) imaging technology has shown great significance in recent years for its high-temporal-resolution 3D imaging feature through simultaneously capturing the 2D spatial and 2D angular information of light [four-dimensional (4D) LF information]^[1–5]. Especially, the LF imaging method based on the wave-optics model and LF point spread function (LFPSF) allows 3D deconvolution for the high-quality single-shot volumetric reconstruction^[6,7]. However, in some applications of LF imaging, a scattering medium is present in the scenes, such as biological tissue, fog, and turbid water^[8–16]. In these imaging scenes, signal light could still be captured, and 3D reconstruction could be conducted, but scattered light produced by the scattering medium introduces blur and scattering background artifacts to the 3D reconstruction image, which lead to low resolution and low contrast. For the LF imaging method based on 3D deconvolution, sequence recorded frames are utilized to extract ballistic light and undo the effect of scattering on the LF in real-time localization of neuronal activity^[10,11], but the demand of multiple frames rather than single-shot reduced the flexibility of the method. Signal light and scattered light are separated and reconstructed separately^[14], but this will increase the amount of calculation and the solving difficulty. In a general method, scattering is incorporated into the LF imaging forward model to deal with the case of mismatching between the scattering-free model and the LF with scattering, and then the scattering background artifacts of 3D reconstruction can be removed^[8,9]. However, scattering introduces degradation of the LF image and produces blur patterns of LFPSF. These make the inverse problem of the 3D deconvolution more ill-conditioned and lead to a noise-sensitive result^[17].

In recent years, deep learning has been widely studied due to its remarkable ability in information processing and has shown great value in the field of computational imaging^[18–24], especially in LF imaging^[25,26] and scattering imaging^[27–30].

In this paper, a deep learning-based method of removal of the scattering effect on LF imaging is proposed to deal with the case of scattering imaging in LF imaging applications mentioned above. We term this method deep learning-based scattering LF imaging (DeepSLFI). In DeepSLFI, we build the LF imaging forward models and generate training samples of LF images with and without scattering by simulation. By the means of simulation, the huge experimental workload for samples capturing can be avoided. Then, a neural network trained with the simulation samples is utilized to remove the effect of scattering on the LF image captured experimentally. Finally, the high-resolution and high-contrast 3D reconstruction can be conducted with the processed LF image. In this paper, a validation experiment is conducted, and the 3D reconstructions of the scene of imaging through a single scattering layer are obtained. Compared to the results of conventional methods, the 3D reconstruction using DeepSLFI has a higher quality in terms of resolution and contrast.

The overview of DeepSLFI is shown in Fig. 1. The building of the LF imaging forward model is conducted with the theory of LF imaging^[6–9,31]. The LFPSFs and the forward model, shown in Fig. 2, which refer to the projection of the 3D volume in the object space to the 2D image space on the sensor, can be obtained with the determined scene and imaging system. The discretized forward model can be expressed as a matrix equation, $$\mathbf{f}=\mathbf{H}\mathbf{g},$$(1)where $\mathbf{f}$ represents the LF image captured by the sensor, $\mathbf{g}$ is the object space volume to be reconstructed, and $H$ is the measurement matrix determined by the LFPSFs (reflecting the light propagation). The element $h_{ij}$ in row $j$ of column $i$ in $\mathbf{H}$ is $$h_{ij}={\iint }_{{\alpha }_j,{\beta }_i}I(\mathbf{x},\mathbf{p})\mathrm{d}\mathbf{p}\mathrm{d}\mathbf{x},$$(2)which refers to the proportion of the light arriving at pixel $j$ propagating from emitters in voxel $i$, where $I(\mathbf{x},\mathbf{p})$ is the intensity at $\mathbf{x}=(x_1,x_2)$ on the sensor, responding to the light emitted from the point emitter located at $\mathbf{p}=(p_1,p_2,p_3)$, ${\alpha }_j$ is the area for pixel $j$, and ${\beta }_i$ is the volume for voxel $i$^[6]. Generally, a micro-lens array (MLA) is utilized to encode the 3D intensity distribution of the volumetric object into LF information and determine the LFPSFs.

In the absence of scattering, light propagation can be described by theories of wave optics and Fresnel propagation^[6,7,31]. The wavefront at $\mathbf{x}=(x_1,x_2,x_3)$ of light emitted from the point emitter located at $\mathbf{p}=(p_1,p_2,p_3)$ is $$U(\mathbf{x},\mathbf{p})=\frac{1}{r}\exp (ikr),$$(3)where $k$ is the wave number, and $r={\Vert \mathbf{p}-\mathbf{x}\Vert }_2$. Fresnel propagation can be expressed by the Fourier transform operator $\mathcal{F}\{•\}$ as $$U_{\text{output}}(\mathbf{x})={\mathcal{F}}^{-1}\{\mathcal{F}\{U_{\text{input}}(\mathbf{x})\}\exp [-\frac{i}{4\pi }\lambda d({\omega }_{\mathit{x}}^2+{\omega }_y^2)]\},$$(4)where $U_{\text{output}}$ and $U_{\text{input}}$ are the wavefront of output plane and input plane (after and before propagation), respectively, $\lambda$ is the wavelength in the vacuum, $d$ is the propagation distance, and ${\omega }_x$ and ${\omega }_y$ are spatial frequencies along the $x$ and $y$ directions, respectively. The lens/lenslet and MLA can be modeled as phase masks ${\varphi }_{\text{lens}}(\mathbf{x})$ and ${\mathrm{\Phi }}_{\mathrm{MLA}}(\mathbf{x})$, respectively: $${\varphi }_{\text{lens}}(\mathbf{x})=\exp \left(\frac{-ik}{2f}{\Vert \mathbf{x}\Vert }_2^2\right),$$(5)$${\mathrm{\Phi }}_{\mathrm{MLA}}(\mathbf{x})={\varphi }_{\text{lens}}(\mathbf{x})*\mathrm{comb}(\mathbf{x}/d_{\text{pitch}}),$$(6)where $f$ is the focal length of the lens/lenslet of the MLA, $d_{\text{pitch}}$ is the pitch of the lenslet of the MLA, $*$ means convolution, and $\mathrm{comb}(•)$ is comb function. The intensity on the sensor plane of light is $I_{\text{sensor}}(\mathbf{x})={|U_{\text{sensor}}(\mathbf{x})|}^2$.

In the presence of scattering, light propagation can be described by an analysis method based on phase-space theory^[8,9,32]. This allows the scattering to be described with a homogeneous scattering approximation to obtain an analytic forward model, which could not be realized by the theory of Fresnel propagation. The 4D phase-space Wigner function is $W(\mathbf{r},\mathbf{u})$^[8], where $\mathbf{r}=(x,y)$ is the transverse spatial coordinates, and $\mathbf{u}=(u_x,u_y)$ is the spatial frequency coordinates. Then, a 2D intensity image can be described as $I(\mathbf{r})=\iint W(\mathbf{r},\mathbf{u}){\mathrm{d}}^2\mathbf{u}$. An incoherence point emitter located at $\mathbf{p}=(p_1,p_2,p_3)$ could be described by the Wigner function of the plane located on depth $p_3$ as $W_s(\mathbf{r},\mathbf{u})=\delta (\mathbf{r}-{\mathbf{r}}_{\mathit{s}})$, where ${\mathbf{r}}_{\mathit{s}}=(p_1,p_2)$. The propagation of distance $\mathrm{\Delta }z$ in Wigner space and the effect of the lens can be modeled as Eqs. (7) and (8), respectively^[32]: $$W_{\text{output}}(\mathbf{r},\mathbf{u})=\iint W_{\text{input}}({\mathbf{r}}^{\prime },\mathbf{u})\delta (\mathbf{r}-{\mathbf{r}}^{\prime }-\frac{\lambda \mathrm{\Delta }z\mathbf{u}}{n_r}){\mathrm{d}}^2{\mathbf{r}}^{\prime },$$(7)$$W_{\text{output}}(\mathbf{r},\mathbf{u})=\iint W_{\text{input}}(\mathbf{r},{\mathbf{u}}^{\prime })\delta (\mathbf{u}-{\mathbf{u}}^{\prime }+\frac{\mathbf{r}}{\lambda f}){\mathrm{d}}^2{\mathbf{u}}^{\prime },$$(8)where $W_{\text{output}}$ and $W_{\text{input}}$ are the Wigner functions of the output plane and input plane, respectively, ${\mathbf{r}}^{\prime }=(x_1^{\prime },x_2^{\prime })$ is the spatial coordinates of the input plane, ${\mathbf{u}}^{\prime }=(u_x^{\prime },u_y^{\prime })$ is the spatial frequency coordinates of the input plane, and $n_r$ is the refractive index. Especially, the intensity at point $\mathbf{r}$ on the back focal plane of the lens/lenslet of the MLA can be approximately represented as $I(\mathbf{r})=W({\mathbf{r}}_{\text{lens}},{\mathbf{u}}_{\text{lens}})$, where $W({\mathbf{r}}_{\text{lens}},{\mathbf{u}}_{\text{lens}})$ is the Wigner function of the lens plane, ${\mathbf{r}}_{\text{lens}}$ is the center of the lens/lenslet, and ${\mathbf{u}}_{\text{lens}}$ approximately corresponds to $\lambda {\mathbf{u}}_{\text{lens}}/n_r=(\mathbf{r}-{\mathbf{r}}_{\text{lens}})/f$. The scattering of the diffusing plane can be described as $$W_{\text{output}}(\mathbf{r},\mathbf{u})=\iint W_{\text{input}}(\mathbf{r},{\mathbf{u}}^{\prime })\mathit{D}(\mathbf{u}-{\mathbf{u}}^{\prime }){\mathrm{d}}^2{\mathbf{u}}^{\prime },$$(9)where $\mathit{D}(\mathbf{u}-{\mathbf{u}}^{\prime })$ is the phase function, which refers to the effect of scattering on directions of the outgoing photons^[33]. A homogeneous scattering approximation is utilized, and the phase function could represent the statical intensity distribution of the outgoing angle information relative to the incident angle information. In Ref. [8], scattering is modeled as an in-place Gaussian spreading of angle information. Here, we propose a modified version of this phase function. The modified version is also based on Gaussian distribution but is more flexible, especially because it can well describe the property of the scattering layer in our experiment later in this paper: $$\mathit{D}(\mathbf{u}-{\mathbf{u}}^{\prime })={\sum }_m{A}_m\exp [-\frac{{(\mathbf{u}-{\mathbf{u}}^{\prime })}^2}{{\sigma }_m^2}],m\in {\mathbb{Z}}^{+},$$(10)where ${\sigma }_m$ is the scattering coefficient, and $A_m$ is the weight factor.

After determining the light propagation in the imaging scene, the intensity proportion on the sensor of the light emitted from point emitter in the object space, $I(\mathbf{x},\mathbf{p})$, can be calculated by Eqs. (3)–(10). Then, $\mathbf{H}$ and the forward model can be built by Eqs. (1) and (2). With the $\mathbf{H}$ and the LF image $\mathbf{f}$ captured experimentally, the 3D reconstruction of $\mathbf{g}$ can be obtained by solving the inverse problem of the matrix equation through 3D deconvolution. Figure 3 shows the blur and the scattering noise introduced by the scattering layer used in the experiment later in this paper. The signal light could still be captured, but the scattered light would result in serious background artifacts of the 3D reconstruction. Once introducing scattering into the forward model, the blurred LFPSFs would make the inverse problem more ill-conditioned and lead to serious noise sensitivity. These would seriously reduce the quality of 3D reconstruction, especially in a complex scenario.

In order to remove the effect of scattering on the LF image, a convolutional neural network (CNN) is used in DeepSLFI. Figure 4 is the overall architecture of the CNN, which is adapted from the residual network (ResNet)^[34] and U-shaped network (U-Net)^[35]. The network consists of eight convolution blocks, each composed of a convolution, a leaky rectified linear unit (LeakyReLU), a batch normalization (BN) layer, and a final tanh activation. The CNN has shown the ability of image deblurring in our preliminary simulation work^[36].

Considering the huge workload and the serious difficulty for the experimental capturing of samples of 3D images or LF images in scattering scenes, the training sample pairs will be generated by simulation with the LF imaging forward models. This method also makes DeepSLFI more flexible. The LF imaging forward models have been applied in research^[6–9] and have good performances in describing the practical LF propagation. Therefore, it is believed that the forward models can be utilized to generate samples. Trained with the LF images with scattering (as input) and the scattering-free LF images (as output), the CNN will effectively remove the scattering effect on the input LF image.

The final reconstruction will be realized by 3D deconvolution using the scattering-removed LF image and the scattering-free forward model (which has sharp patterns of LFPSFs). In this paper, CNN is utilized to remove the scattering of the LF image, not to produce 3D image results of the inverse problem end-to-end^[25,26]. The 3D reconstruction will be conducted by an iterative algorithm, such as the Richardson–Lucy algorithm^[6], which is the proven technique, and we are more adept at.

An experiment is conducted to verify the validity of DeepSLFI. Compared with the 3D reconstruction results obtained from different LF imaging methods, the advantages of DeepSLFI in LF imaging through scattering could be demonstrated.

The experimental system is shown in Fig. 5, which consists of a CCD sensor (BFLY-PGE-50H5M-C, Point Grey, 3.45 µm pixel size and $2\times 2$ pixels binning employed, the LF image will be a $10\times 10$ image array with a $730\times 730$ resolution), an MLA (#64-480, Edmund Optics, 0.5 mm lenslet pitch and 15.3 mm focal length), a lens (focal length 100 mm), a square aperture with the width of 3.27 mm (used to prevent aliasing of subimages of the LF image), and a 0.5 mm thick scattering layer (introduced as a single diffusing plane that can be added or removed as demanded). The imaging object will be a white light source. The wavelength of 632.8 nm is selected to build the forward models. In this experiment, the light mainly propagates in free space, thin optical devices, and thin scattering layer. In addition, the numerical aperture of the system is small. Therefore, the dispersion of the experimental system is very small, and the mismatch between the monochromatic forward models and the white light imaging scene can be ignored. In Ref. [6], the monochromatic forward model also has a good performance in imaging with white light sources.

The effect of the scattering layer on the LF image is shown in Fig. 3. Considering that the intensity curves in Fig. 3(c) are smooth in general, and the homogeneous scattering approximation in research of LF imaging through scattering has demonstrated its validity^[8,9], the homogeneous scattering is a very good approximation and will be used in our modeling process. In addition, a homogeneous scattering approximation allows for an analytic forward model for scattering.

By calculating the light intensity distributions on the sensor after the propagation shown in Fig. 6, the LFPSFs and the forward models can be obtained with Eqs. (1)–(10). A set of parameters of the scattering model, which can well describe the scattering effect of the scattering layer in the experiment, is used. Generally, the effect of aberration can be negligible in the forward model^[6–9,31], especially because the numerical aperture and field of view of the imaging system in this experiment are small. The field of view of the system is a $5.1\mathrm{mm}[\mathrm{height}(\mathit{H})]\times 5.1\mathrm{mm}[\mathrm{width}(\mathit{W})]\times 25\mathrm{mm}[\text{depth}(\mathit{D})]$ cube, where the front plane of the depth of field is 145 mm away from the lens, as shown in Fig. 5. The resolution of the voxel image is $211(\mathit{H})\times 211(\mathit{W})\times 25(\mathit{D})$, and the sampling rates are uniform at different depths.

The simulation training samples are generated based on the scene: single-layer planar objects located on 25 depths (1 voxel thickness) with intensity maps of the Modified National Institute of Standards and Technology (MNIST)^[37] images. For each depth, 100 pairs of training sets and 10 pairs of test sets (each pair consists of the LF images with and without scattering of the simulation object) are generated. In total, 2500 groups of training samples and 250 groups of test samples are generated. To avoid overfitting and enhance the generalization of the neural network, random factors are added to the samples to simulate more imaging conditions. The details of the generation process of the samples are shown in Fig. 7.

Before the neural network training process, all of the generated LF images are scaled to [−1,1] (the range of tanh activation function). Each convolution layer of the network has a stride of 1 and kernel size of $3\times 3$, and the number of filters is shown in Fig. 4. The slope parameter of all LeakyReLU activations is set to 0.3. The model is trained end-to-end with a batch size of four. An Adam optimizer with an initial learning rate of 0.01 and decay rate of 0.5 is utilized. Mean squared error (MSE) is used as the loss function. The training is performed with one graphics processing unit (GPU, NVIDIA Tesla T4) using Keras/TensorFlow. The network training is converged after 200 epochs and about 40 h.

LF images of the object are captured experimentally with and without scattering, respectively. The experiment is conducted under the same condition except for the presence or absence of the scattering layer. The 3D deconvolutions (using the Richardson–Lucy iteration scheme) of the same object are performed in four ways: using the LF without scattering and the scattering-free forward model (recorded as ${\mathbf{R}}_{\mathrm{ns}}$), using the LF with scattering and the scattering-free forward model (recorded as ${\mathbf{R}}_{\mathrm{s}\mathrm{0}}$), using LF with scattering and the forward model with scattering (recorded as ${\mathbf{R}}_{\mathrm{s}\mathrm{1}}$), and using LF with scattering and DeepSLFI (recorded as ${\mathbf{R}}_{\mathrm{DeepSLFI}}$). Several 3D reconstructions of the same scattered LF image in the different ways are compared, where the result of ${\mathbf{R}}_{\mathrm{ns}}$ is taken as the reference value. We adapted the MATLAB codes provided by Refs. [7,31] for building the forward models and the 3D deconvolution according to the demands of this experiment. The iteration number to achieve convergence of deconvolution is 40. The processing is completed on MATLAB R2019b.

Figure 8 shows the results of a United States Air Force (USAF) 1951 target (Group 0, Element 6, attached to a white light source to get even illumination) located in the depth of field with a distance of 10 mm from the front plane of the depth of field. To quantitatively analyze the reconstruction results, the peak signal to noise ratio (PSNR) and structural similarity (SSIM)^[38,39] are utilized to evaluate the similarity between the results of scattered LF (scattered LF, deblurred LF, and the 3D reconstructions of them) and the scattering-free ones (serving as the reference results). The network effectively removes the blur of the LF led by scattering. This result demonstrates the validity of the homogeneous scattering approximation and the simulation method of CNN training. In addition, the trained CNN using LF images of MNIST as the samples can also deal with the LF images of the USAF target, which indicates the generalization ability. As can be observed from the 3D reconstruction results, the result of ${\mathbf{R}}_{\mathrm{ns}}$ has the lowest degree of spread in both lateral and axial directions. The result of ${\mathbf{R}}_{\mathrm{s}\mathrm{0}}$ has serious scattering background artifacts, severely disturbing the observation. When incorporating scattering into the forward model, the reduction of scattering background artifacts is obvious in the result of ${\mathbf{R}}_{\mathrm{s}\mathrm{1}}$. However, the spread in both lateral and depth directions is serious, leading to a low resolution. In the result of ${\mathbf{R}}_{\mathrm{DeepSLFI}}$, the scattering background artifacts are removed obviously, and the spread of the 3D image both in the lateral and depth directions is reduced significantly. These indicate that the results of ${\mathbf{R}}_{\mathrm{DeepSLFI}}$ have both the highest resolution and the highest contrast and are closest to the result of ${\mathbf{R}}_{\mathrm{ns}}$.

Figure 9 shows the 3D reconstruction results of two slits localized at different depths in object space. It shows the ability of 3D reconstruction through scattering of DeepSLFI.

In this paper, we propose DeepSLFI to remove the effect of scattering on LF imaging. DeepSLFI utilizes a CNN, trained by simulation samples, which are generated by the LF imaging forward models, to remove the effect of scattering on the LF image captured experimentally and then realizes high-quality 3D reconstruction using the scattering-removed LF image and the scattering-free forward model. The advantage of DeepSLFI in high-quality imaging in terms of resolution and contrast is demonstrated experimentally. The future research of DeepSLFI will focus on the applications with more different and complex imaging scenes, such as scenes with volumetric scattering and stronger scattering. Future research will also focus on specific applications, such as 3D imaging in biological tissue and fog.