Extended depth-resolved imaging through a thin scattering medium with PSF manipulation

Human ability to visualize an image is usually hindered by optical scattering. Recent extensive studies have promoted imaging technique through turbid materials to a reality where color image can be restored behind scattering media in real time. The big challenge now is to recover objects in a large field of view with depth resolving ability. Based on the existing research results, we systematically study the physical relationship between speckles generated from objects at different planes. By manipulating a given single point spread function, depth-resolved imaging through a thin scattering medium can be extended beyond the original depth of field (DOF). Experimental testing of standard scattering media shows that the DOF can be extended up to 5 times and the physical mechanism is depicted. This extended DOF is benefit to 3D imaging through scattering environment, and it is expected to have important applications in science, technology, bio-medical, security and defense.

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 system 22,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 range [24][25][26][27] . The imaging process can be denoted as a correlation function of the PSF and the object distribution function 4 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 process 4,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 point 4 . 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 mask 28 TM(x s , y s ). The light field distribution on the imaging plane is 24 , The new PSF is equal to the original one except the scaling in transform variable, i i i i 2 m is the scaling factor, and it is written by 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 fluctuations 27 .
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.
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 as 29 , i 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 2 nd 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 1 st 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)).
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 algorithm 4 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.
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  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 algorithm 30 . 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 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, is the equivalent NA for the imaging processing with a scaling lens.
The nominal DOF 21,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λ/M 2 NA′ 2 = 7.3 mm, which is closed to the original DOF (7.2 mm) in our experiment, where 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 lens 4 . 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 cameras 33 (using pseudo-random phase mask with convex bumps) and of lensless microscopy 34 (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 paths 28 .

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 fluctuations 23 . 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 Scientific REPoRtS | (2018) 8:4585 | DOI:10.1038/s41598-018-22966-7 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.