Introduction

Light scattering is an obstacle for optical imaging when our sights are blocked by smoke, mist, anisotropy biological tissues, or air turbulence. As science and technique progress, we have witnessed steady development for optical imaging recovery through scattering media. Among various novel technologies, guide star assisted wavefront shaping1,2,3,4,5, holography6, scattering matrix measurement7 and speckle correlation2,8,9 are the most effective. However, these technologies are valid mainly for detecting the planer shape of the hidden object, where its depth and location information is rarely retrieved. The depth information, or 3-dimensional (3D) imaging, is vital for a wide range of real-world applications. For example, the exact location of the tumor and its axial structure should be determined before any therapy is operated. Furthermore, added axial information will greatly assist the determination of the nature of an object. To achieve this objective, a great deal of effort has been made and new techniques have been introduced in recent years. The following are the most representative ones, to just name a few: ultrafast time-of-flight scattering imaging method10, digital holography technique11, phase-space analysis methods12,13, holographic image sensor with speckle-correlation scattering matrix14, DiffuserCam with a well-designed compress sensing algorithm15, depth-resolved with speckle holography and two-point intensity correlation16, wavefront shaping with axial “Memory Effect” (ME)17,18,19,20, speckle correlation method with an assisted reference point source21, and so on. Among these improved techniques, the speckle correlation method appears very promising, as it does not need any coherent sources, raster-scanning mechanism, nor any time-consuming wavefront shaping devices. As reported in our last paper4, a well-designed optical system generates a Point-Spread-Function (PSF) which is helpful for high speed color imaging through scattering media with a large Field-Of-View (FOV). However, similar to most of the imaging techniques based on ME, they commonly suffer from a limited DOF in axial detection. One straightforward solution is capturing image slices in different axial positions supplemented with a series of corresponding PSFs. However, determination of PSFs for different objective planes complicates the system hence recovery of 3D image appears low and less practical. It is thus fundamentally important to study the physical relationship between the depth of the object and its speckle properties.

Speckle correlation has been well discussed on the imaging plane, from either lateral or axial direction. But the contribution of the incident light filed is seldom considered in textbooks and not detailed depicted in related literatures. Here, we study the speckle properties after a thin scattering medium in different incident conditions and found that a single PSF is sufficient to quantitatively retrieve imaging through the scattering media even beyond the original DOF. Prove-of-concept is demonstrated both with simulation and with experiment that the retrieved imaging is shown with good quality and high fidelity. Experimental results show a 5 times extension of z-axis range with a single PSF compared with original DOF. A rigorous theoretical analysis is presented to support the numerical simulation results.

Principle

For an imaging system, a point at the object plane generates a nearly identical pattern on the image plane (PSF), which represents the properties of the imaging system22,23. The intensity distribution of the light field under incoherent illumination is described with a convolution, I = O*PSF, where I and O are intensity distribution on the image plane and the object plane respectively, the symbol * denotes a convolution operation.

When a thin scattering medium is introduced into the imaging system, the PSF becomes a speckle pattern determined by the property of the scattering medium. This scattering pattern is shown to be shift-invariant within the ME range24,25,26,27. The imaging process can be denoted as a correlation function of the PSF and the object distribution function4,24:

$$I({x}_{i},{y}_{i})=F(k{x}_{i}/{d}_{i},k{y}_{i}/{d}_{i})\iint PSF({x}_{i},{y}_{i};{x}_{o},{y}_{o})O({x}_{o},{y}_{o})d{x}_{o}d{y}_{o},$$
(1)

where k is the wavevector, the subscripts i and o denote the image and the object plane respectively. F(kx i /d i , ky i /d i ) is defined as a form factor function which acts as a field stop placed in the focal plane to limit the observed FOV. The convolution equation can be rewritten in the form of I = F•(O*PSF). Although the PSF in this case is very complicated, the deconvolution process4,27 can restore O by considering the imaging system as a “black box”. Obviously, deconvoluted imaging is valid only if the reference point and the object are on the same plane (focus). The image will become blurred or even invisible as long as the object deviates from the reference point4. The deviated distance between which the image can be restored is determined by the DOF.

For better understanding the relationship between different PSFs, as shown in Fig. 1, the complicated PSF is calculated by treating the thin scattering medium as a random phase mask28 TM(x s , y s ). The light field distribution on the imaging plane is24,

$$h({x}_{i},{y}_{i};{x}_{o},{y}_{o})=Fr{\{TM\cdot Fr{\{\delta ({x}_{o},{y}_{o})\}}_{do}\}}_{di},$$
(2)

where Fr{} di and Fr{} do represent Fresnel diffraction with distance of d i and d o respectively. The intensity of PSF under incoherent illumination is recorded, which can be written as the square of the module of the light field: \(PSF({x}_{i},{y}_{i};{x}_{o},{y}_{o})={|h({x}_{i},{y}_{i};{x}_{o},{y}_{o})|}^{2}\).

Figure 1
figure 1

Schematic of deconvolution 3D imaging beyond DOF limit through a scattering medium. Virtual PSFs from virtual point (green) can be calculated with PSF from a real pinhole (red).

Full size image

By applying the autocorrelation theorem and rewriting equation with a form of Fourier transform. The PSF becomes,

$$PSF({x}_{i},{y}_{i};{x}_{o},{y}_{o})= {\mathcal F} {\{{ {\mathcal F} }_{a}^{-1}{\{A(\alpha ,\beta )\}}_{\alpha \to \frac{{x}_{s}}{2\lambda f},\beta \to \frac{{y}_{s}}{2\lambda f}}\}}_{{x}_{s}\to \frac{{x}_{i}}{\lambda {d}_{i}},{y}_{s}\to \frac{{y}_{i}}{\lambda {d}_{i}}},$$
(3)

where \(A(\alpha ,\beta )=\sqrt{1/4i\lambda f}{e}^{-\frac{ik}{8f}({\alpha }^{2}+{\beta }^{2})}\) * \([TM(\alpha ,\beta )T{M}^{\ast }(\alpha -{x}_{s},\beta -{y}_{s})]\), \( {\mathcal F} \{\,\}\) is a symbol of Fourier transformation, f is the focal length of the “scattering lens” which meets the formula of object-image distance (1/f = 1/d i  + 1/d o ).

When the object distance do is translated to d′ o , the new PSF becomes,

$$PSF^{\prime} ({x}_{i},{y}_{i};{x}_{o},{y}_{o})= {\mathcal F} {\{{ {\mathcal F} }_{a}^{-1}{\{A(\alpha ,\beta )\}}_{\alpha \to \frac{{x}_{s}}{2\lambda f^{\prime} },\beta \to \frac{{y}_{s}}{2\lambda f^{\prime} }}\}}_{{x}_{s}\to \frac{{x}_{i}}{\lambda {d}_{i}},{y}_{s}\to \frac{{y}_{i}}{\lambda {d}_{i}}},$$
(4)

where the transform variables in \({ {\mathcal F} }_{a}^{-1}\) is changing to \(\alpha \to \frac{{x}_{1}}{2\lambda f^{\prime} }\) and \(\beta \to \frac{{y}_{1}}{2\lambda f^{\prime} }\), 1/f′ = 1/d i  + 1/d′ o .

The new PSF is equal to the original one except the scaling in transform variable,

$$PSF^{\prime} ({x}_{i},{y}_{i})={m}^{2}PSF(m{x}_{i},m{y}_{i}).$$
(5)

m is the scaling factor, and it is written by

$$m=\frac{f}{f^{\prime} }=\frac{({d}_{i}+{d}_{o}){d^{\prime} }_{o}}{({d}_{i}+{d^{\prime} }_{o}){d}_{o}}.$$
(6)

A more detailed derivation of equations (2)–(6) can be found in the Supplementary Information.

Based on equation (5), the difference between PSFs in various object planes is only on the size scale. It turns out to be an interesting and important conclusion: once a PSF is obtained, virtual PSFs of other object planes can be calculated by resizing the PSF with a scaling factor m.

Prove-of-concept simulation

In an attempt to verify the above theory, a numerical simulation is implemented without using the above theoretical result. As shown in Fig. 1, the propagation of the light field from a pinhole through a scattering medium to the imaging plane is commonly described as two steps of Fresnel diffraction. The first step is Fresnel diffraction from the pinhole to the plane in front of the scattering medium. The second step is Fresnel diffraction from the scattering medium to the imaging plane. The light field after the scattering medium dU′ s is equal to the one before multiplied by the transmission matrix of the scattering medium as, dU′ s  = dU′ s •TM, where the phase mask TM(x s , y s ) applied in the simulation is created from a random phase matrix convoluted with a 2D Gaussian function and thus, determining the frequencies of its spatial fluctuations27.

By setting the pinhole as a 50 μm diameter round disk and applying to standard Fresnel diffraction algorithms, where both S-FFT for Fresnel diffraction propagation and D-FFT for angular spectrum propagation are tested, the intensity of the light field on the imaging plane is worked out as shown in Fig. 2. Obvious shrinking of speckle pattern can be observed when d o of pinhole is increasing. A scaling factor of m can prevent the speckle pattern from shrinking.

Figure 2
figure 2

The speckles on the imaging plane derived from a pinhole located on different object planes. The first row is the original calculated results with green rectangles selecting the region of interest and blue squares indicating the area similar to the left-most image. The second row shows the images that are m × m times the size of the area in the blue squares.

Full size image

On the other hand, the size of the speckles (also known as resolution) has been considered as relating to the aperture angle that radiation from the scattering medium as29,

$$\sigma \approx 1.2\lambda {d}_{i}/D=0.61\lambda /NA,$$
(7)

where D is the diameter of the scattering medium and NA is numerical aperture of the scattering lens. Here we have to emphasize that σ is also related to the object distances (the input light field). When d o is translating from 10 cm to 40 cm (d i  = 10 cm), the whole speckle pattern is shrinking. One of the brightest spot in green rectangle is shifting towards the image center. The related area (in blue squares) is zoomed in by a factor of m = f/f′ and cut out (on the 2nd row of Fig. 2). The zoomed speckles show the same appearance with the reference speckle pattern (d o  = 10 cm), which is consistent with the theoretical result.

Experiment

An experimental setup for imaging through a scattering medium is presented in Fig. 1. A physical pinhole illuminating by an incoherent light source (1 W green LED source by Daheng Optics, GCI-060403) is placed on the optical axis of the imaging system. Its diffractive light projects onto a standard scattering medium (Newport 5° circular light shaping diffuser) and diffuses. A monochrome CCD (Basler ACA2040-90UM) is applied to capture the diffused light, where d i  = 7.5 cm. By translating the pinhole along the z direction from d o  = 11 cm to d o  = 24 cm with an interval of 0.5 cm, a series of PSFs is recoded (obviously zooming effect on the speckles can be seen in the supplementary movie).

By selecting PSF of d o  = 15 cm as the reference, deconvolution of different PSFs with the reference PSF are executed. The results are shown in the 1st row of Fig. 3, which represent the reconstructed image of the pinhole in the corresponding object plane. The reconstructed images appear defocused and blur when their object planes are not exactly on the central object plane. The peak of the correlated images rapidly declines as the pinhole deviates from the central position. Its value drops down to 0.5 when the pinhole is translated to z~±3.6 mm, as shown in the Fig. 3(a). Hence the original DOF of deconvolution imaging technique is 7.2 mm (FWHM of Gaussian fit of Fig. 3(a)).

Figure 3
figure 3

Deconvolution between PSFs from different d o with a selected PSF from central object plane (d o  = 15 cm). The selected PSF in (a) is the original one and in (b) are virtual PSFs after scaling, the inserted number indicates different d o in millimeter. (c) The cross section through the peaks of (a) black line and (b) red line. (d) Peak values of the intensity correlations.

Full size image

By rescaling the selected PSFs according to equation (6), the correlated images remain ‘focused’ over large distances with peak values declining slowly. Hence the axial range of retrieved imaging is improved to be 36.6 mm (FWHM of Gaussian fit of Fig. 3(b)), which is 5 times as large as the original DOF (shown in Fig. 3(d)).

For image restoration, the pinhole is replaced with unknown tested objects after the reference PSF (d o  = 15 cm) is captured. The tested object (letter mask of ‘H’, about 1 mm tall) is located on different object distances (12~21 cm, step of 0.5 cm) once at a time and their speckles are captured separately. The reconstructed images are calculated by a deconvolution algorithm4 with the original reference PSF. As shown in Fig. 4(a), the image of ‘H’ is resolvable only when the object locates within the original DOF range. By rescaling the reference PSF according to equation (6) (shown in Fig. 4(b)), the imaging of ‘H’ can be resolved clearly from d′ o  = 12 cm to 19 cm. This range is even larger than the improved DOF (36.6 mm) because the rescaled PSFs are more relevant to the real PSFs of their object planes (at d′ o ). Their correlation images would be similar to Fig. 3(b) and remain ‘focused’ over larger distances than the improved DOF range.

Figure 4
figure 4

Deconvolution imaging without (a) and with (b) (c) derived virtual PSFs. Objects in (a) (b) is a letter ‘H’ moving along z axis and in (c) are 3 letters located at different object plane simultaneously. The top left corners of each restored images are setting transparent for better view of the blocked layers.

Full size image

When DOF of the imaging system is extended to a sufficiently large range, 3D imaging can be achieved slide by slide. Three different letter masks of ‘G’, ‘H’ and ‘N’ are placed on three separated object planes, as shown in Fig. 4(c). Their object distances are 13 cm, 15 cm and 17 cm respectively. Letters ‘G’ and ‘N’ are obviously located beyond the original DOF range with the only known PSF at d o  = 15 cm. The speckle pattern of these three letters is captured by the CCD. By rescaling the size of the reference PSF and correlating it with the objects speckle, images of the letter masks are restored separately with different scaling factor. The resolution of z-direction in this technique is the same as the original DOF (7.2 mm). With the right scaling factors, images corresponding to objects at 13 cm, 15 cm and 17 cm can be restored clearly. While recovered images corresponding to object distances of 14 cm and 16 cm show a mixture of the adjacent objects (‘G’ and ‘H’, or ‘N’ and ‘H’). That is because the rescaled PSF corresponding to object distant of 14 cm gets part of correlation with speckles from ‘G’ or ‘H’. The recovered image shows a blur shape of mixture of ‘G’ and ‘H’. Same situation appears at 16 cm.

Furthermore, the proposed method can be used to determinate the location of a tested object. When the object distance is unknown, its image can be restored with an altering scaling factor searched by an adaptive algorithm30. The right scaling factor results in a restored image with highest contrast. And finally, the corresponding object distance is calculated according to equation (6).

Discussion

The relationship between the PSFs of scattering media from different reference points can be essentially derived from the angular ME. Considering a pinhole moves towards the scattering medium along z axis direction from the plane at infinity, the wavevector from the pinhole to the same area of the scattering medium will tilt towards the vertical direction. Hence, its intrinsic isoplanatism on the imaging plane will expand outward but still ‘correlated’ (size should be rescaled) following the angular ME. However, the max tilt angle before ‘uncorrelated’ is limited by ME angle, so there still is a limitation of Δd o  = d′ o  - d o , namely improved DOF.

The original resolution (equation (7)) of an imaging system based on scattering medium is derived from plane wave incident (red line), in which cases d o is infinite and f∞ = d i . In the process of imaging, the object distance should no longer be infinite, the size of the speckle pattern as well as the resolution would expand with m times (m = f∞/f), as presented in Fig. 2. Equation (7) becomes,

$$\sigma =m\frac{1.2\lambda {d}_{i}}{D}=\frac{1.2\lambda {d}_{i}^{2}}{Df}=\frac{0.61\lambda }{NA^{\prime} },$$
(8)

where \(NA^{\prime} =\frac{Df}{2{d}_{i}^{2}}\) is the equivalent NA for the imaging processing with a scaling lens.

The nominal DOF21,31, which depicted the axial correlation of the speckles in the imaging plane, equals to about 6.7λ(z/D)2. In our case, it should be modified with NA′ and be deviated to the objective plane. It turns out to be DOF = 6.7λ/M2NA′2 = 7.3 mm, which is closed to the original DOF (7.2 mm) in our experiment, where M = d i /d o is the transverse magnification.

The improved DOF is calculated with deconvolution (or correlation) processing, so it should be limited by the axial FOV. The form factor F in equation (1) indicates the transverse FOV. It is commonly defined as FOV ~ λd o /πL and measured by correlating the speckle patterns of a transverse shifting point. Katz et al.32 have deduced the axial FOV of a scattering lens compensated with a spatial light modulator from its transverse FOV. A discussion on calculating the axial FOV of a scattering medium with our rescaling method can be found in the supplementary. The maximum axial FOV is estimated as Δz FOV  = 1.4d o FOV/D~61 mm, which is larger than our improved DOF 36.6 mm. Several factors in the experiment degrade the correlation of the speckles and shrink the improved DOF, including low speckle contrast under broadband spectrum illumination, transverse deviation of the object (Fig. 3(b)), low sampling for the small size speckle, discretization error when rescaling and etc. Light source generated by laser passing through high speed rotating diffuser and a digital camera with smaller pixel pitch will make the DOF of our method closer to the axial FOV limit. One simpler way to enlarge the DOF is adding a diaphragm before the scattering medium. Another method is to decrease the angular difference by inserting lens4. Although the improved DOF is still limited by the axial FOV, the dependence of speckle size on the objective distance in equation (8) are valid (expects the complicated “near field speckles” regime). It is possible to separate or extract the speckle patterns of the objects placing beyond axial FOV and to realize a 3D Non-invasive single-shot imaging. Our method can also improve the imaging quality of lensless cameras33 (using pseudo-random phase mask with convex bumps) and of lensless microscopy34 (using a scatter-plate with a thin scattering layer) and adapt them to more complex diffusers. Solid-state samples with different strength of scattering (commercial standard scattering samples from 5 deg to 30 deg), strong forward biological scattering samples (1 mm chicken tissues) are tested with the proposed method. The results confirm that the DOF can be enlarged for more samples whose thicknesses are larger than several mean free paths28.

Conclusion

The correlation of the scattering light in axial direction is studied and exploited. We depict the relationship between the PSFs of a thin scattering medium from different reference points. By adjusting the scaling of one PSF, other PSFs from different object planes can be deduced. Theoretical derivation, prove-of-concept simulation, and experiment for analyzing this relationship are demonstrated with good quality and high fidelity. Experimental results show an approximate 5 times of improvement in depth resolved ability compared to the original method without PSF manipulation. 3-layer objects located separately beyond the original DOF is reconstructed through a thin scattering medium slide by slide with the proposed single PSF deconvolution technique. For a known reference object place, our technique can be even made use of detecting the distance in axial direction through a thin scattering medium.

Methods

The numerical simulation method includes two steps of Fresnel diffraction, as shown in Fig. 1, the first step is from the pinhole to the plane in front of the scattering medium. The second one is from the scattering medium to the imaging plane. The phase mask TM(x s , y s ) applied in the simulation is created from a random phase matrix convoluted with a 2D Gaussian function and thus, determining the frequencies of its spatial fluctuations23. By setting the pinhole as a 50 μm diameter round disk and applying to standard Fresnel diffraction algorithms, where both S-FFT for Fresnel diffraction propagation and D-FFT for angular spectrum propagation are tested, the intensity of the light field on the imaging plane is shown in Fig. 2.

The experimental setup as shown in Fig. 1, A physical pinhole illuminating by an incoherent light source (1 W green LED source by Daheng Optics, GCI-060403) is placed on the optical axis of the imaging system. Its diffractive light projects onto a standard scattering medium (Newport 5° circular light shaping diffuser) and diffuses. A monochrome CCD (Basler ACA2040-90UM) is applied to capture the diffused light, where d i  = 7.5 cm. To reveal the physical relationship between the s PSFs from different object planes after a thin scattering medium, a pinhole is translating along the z direction from d o  = 11 cm to d o  = 24 cm with an interval of 0.5 cm, a series of PSFs is recoded. By selecting PSF of d o  = 15 cm as the reference, deconvolution of different PSFs with the reference PSF are executed, as shown in Fig. 3.

For image restoration, the pinhole is replaced with unknown tested objects after the reference PSF (d o  = 15 cm) is captured. The tested object (letter mask of ‘H’, about 1 mm tall) is located on different object distances (12~21 cm, step of 0.5 cm) once at a time and their speckles are captured separately. By rescaling the reference PSF according to equation (6) (shown in Fig. 4(b)), the imaging of ‘H’ can be resolved clearly from d′ o  = 12 cm to 19 cm with an improved DOF (36.6 mm). When DOF of the imaging system is extended to a sufficiently large range, 3D imaging would naturally be achieved. Three different letter masks of ‘G’, ‘H’ and ‘N’ are placed on three separated object planes, as shown in Fig. 4(c). Their object distances are 13 cm, 15 cm, and 17 cm respectively. Letters ‘G’ and ‘N’ are obviously located beyond the original DOF range with the only known PSF at d o  = 15 cm. The speckle pattern of these three letters is captured by the CCD. By rescaling the size of the reference PSF and correlating it with the objects speckle, images of the letter masks are restored separately with different scaling factor.