High-resolution imaging through randomly complex scattered fields and highly scattering walls is an extended sought-after capability with potential applications in various fields, for instance, underwater imaging, biomedical imaging, and seeing through fog.1^,2 Several methods have been proposed to unscramble the object information from the degraded scattered image, which employs dark channel prior,3^,4 time-gated,5^,6 polarization descattering,7^–9 and others. The application of such methods is demonstrated in imaging through scattered media and brings about the extreme improvement of the image contrast of degraded scenarios. Unfortunately, the most promising scattering imaging method remains a complicated problem of spatial-resolution enhancement for the reason that the flux with object information of multiply scattered light is extremely low. High-resolution object information is scrambled in the scattered field but not completely lost, which can be recovered by reversing the scattering process. Recent advances in deconvolution techniques based on the point spread function (PSF) have made this task possible.10^–12 Regarding deconvolution, a prerequisite is the acquisition and calibration of intensity PSF of the scattering imaging system. Nowadays, a commonly used way to access the PSF is to measure it directly with a point source.13^,14 However, a natural point light source in the biological sample is rarely available, and planting an artificial probe is impractical.15^,16 It is possible to estimate the PSF without a point source but by reconstructing the pupil function, which contains more information than the intensity PSF.17 However, recovering this pupil function requires multiple frames, which would sacrifice the temporal resolution and is impractical for complex scattering media. Generally speaking, the challenge of being time-consuming and prior information-based is the key to developing high-resolution imaging through a complex scattering media.

An interesting and practically important question naturally arose from the dynamic characteristics of the complex scattering media, i.e., the multiple scattering effects and molecular Brownian motion lead to the degradation of the scattering field correlation.18 A significant point to consider is that Brownian motion can produce varied light patterns, providing essential high-frequency information for imaging.19 Shifting of high-frequency components to the passband of the collection optics occurs through frequency mixing with the sample and the non-uniform illumination pattern.20^–22 Consequently, the recorded image contains information from the sample beyond the diffraction limit.23 In general, although the scattering media such as turbid water or biological tissue degrades the image, it also provides the possibility of high-resolution imaging by introducing high-frequency content. For instance, pattern-illuminated Fourier ptychography (PIFP), a typical iterative optimization algorithm, uses a non-uniform intensity pattern for simple illumination and captures the corresponding images with different-frequency information.24^–27 It is worth noting that the scattering media is placed in the illumination path to generate speckle illumination for typical PIFP configuration, and the intensity distribution and relative displacement information of the illuminated patterns should be measured as a priori information,28^,29 which requires a higher degree of accuracy and complexity of the imaging device. Implementing the method becomes problematic when dealing with complex scattering materials in real-world situations. Recently, there has been a surge in research papers centered on employing iterative optimization algorithms to tackle the issue of imaging through scattering media. Li et al.30 applied the ptychographic iterative engine to extend the field of view in speckle-related imaging, which is constrained by the optical memory effect. Pei et al.31 proposed a super-resolution macroscopic imaging method for distant objects based on the Fourier-domain shower-curtain effect, with the addition of ptychography to broaden the field of view. Hu et al.32 proposed an improved ResNet network to extract and enhance the features from distorted intensity images in a hazy environment. However, the above approach either requires prior information or is limited by the optical memory effect.

Here, we present a novel polarized iterative optimization approach for high-resolution imaging through complex scattering media. Only a polarization camera is needed to capture four images of polarizer intensity, as depicted in Fig. 1(a). Instead of designing an illumination pattern and accurately documenting the motion law between the sample and the scattering medium, a non-uniform illumination pattern can be achieved by the Brownian motion of a complex scattering medium. Subsequently, the illumination pattern can be separated from the target information based on the polarization common-mode suppression property.33^,34 In a manner akin to the original PIFP technique,35 we utilize the scattering media images and the low-resolution images illuminated with patterns as the input for iterative optimization, continuing until the sample estimate reaches convergence. Furthermore, the robustness in different samples and complex dynamic scattering media demonstrates the effectiveness of this proposed technique.

The proposed setup for high-resolution imaging through a complex scattering medium is visually represented in Fig. 1(a). The object is hidden behind (or immersed) a highly scattering media closely. Meanwhile, the object is illuminated by a full-perpendicular polarized parallel light source shaped by the polarization state generator, and a polarization camera is utilized to record the sub-polarized images with four polarization filtering angles. Significantly, this proposed method uses multiple intensity-varied pattern Iscat,n generated by Brownian motion to illuminate the object Iobjgt and acquires the corresponding raw low-resolution images Ii through the objective lens [lens 2 in Fig. 1(a)].

Despite the low contrast and apparent randomness of the captured images, there are notable differences in the information distributions of sub-polarized images. It is noteworthy that the intensity variation trend of scattering media and objects is also significantly different. Simply put, this different polarization property between the object and the scattering media provides a priori information to demix and obtain the intensity distribution of scattering media. Here, we used a polarization camera to capture two orthogonal polarization image pairs Ii (i=0,  45, 90, and 135 deg), which are used as the input to a demix algorithm based on the common-mode rejection of polarization characteristics.

After the process of data fitting, the trend of intensity variation of either the target (red block) or the background (green and blue blocks) conforms to the cosine function, as shown in Fig. 1(b), which is highly coupled with the Malus law. It is worth noting that the minimum intensity value is always greater than 0. In other words, the image captured by the polarization camera is partially polarized. Based on the above analysis, Ii also can be written as: Ii=IP+INP, where IP and INP denote the polarized and non-polarized parts, respectively. More importantly, the intensity value of the polarized part IP is obviously larger than the non-polarized part INP in Fig. 1(b). However, it is a fact that the non-polarized parts can reduce the signal-to-noise ratio. We introduce the polarization difference images (PDIs) to suppress the trouble according to the common-mode rejection of polarization characteristics, as shown in Eq. (1)36[PSIPDI]=[111−1]×[ImaxImin].(1)

For an ideal linear polarization analysis system, polarization summation images (PSIs) are equivalent to a traditional total intensity image (Imax+Imin), and PDI represents Imax−Imin. Let Imin be the frame with the least amount of scattering light. The other frame, Imax, is brighter due to an increased scattering light component. These two orthogonal frames are utilized to maximize the retention of scattering information while eliminating the non-polarized portion.

As we all know, two special frames Imax and Imin are difficult to obtain accurately and directly by rotating the polarizer, especially through the complex scattering medium. On the one hand, the scattering effect caused by the scattering medium balances the optical-field information. The information difference in the equalized optical field is too small to distinguish when we rotate the polarizer. On the other hand, the maximum intensity information and minimum intensity information on each pixel are obviously different (including the intensity value and the position of the maximum value), which is caused by the different scattering degrees.

Fortunately, we use the degree of linear polarization (DoLP) as a bridge to calculate the components of Imax and Imin from the Stokes vectors. Basically, for the raw intensity image I=Imax+Imin captured by the camera, we can divide it into object Iobj and background scattering information Iscat components, as shown in Eq. (2) {Imax=Imaxobj+ImaxscatImin=Iminobj+Iminscat.(2)

To display the variation trend clearly, six different positions (including pure background scattering and background mixture object information) are selected to plot the mean intensity value in Figs. 1(c1) and 1(c3). Generally, the intensity in all six positions exhibits the cosine function variation relationship. For the dynamic scattering media, however, the peaks and troughs of the wave are obviously different, even though several curves are opposite, as shown in Fig. 1(c2). It is caused by the fact that the object and background have different modulations of the polarized light field. This has resulted in two special frames Imax and Imin captured by rotating the polarizer and losing the global optimum property. For the static scattering media, furthermore, although the variation trend of curves is consistent, the difference remains in Fig. 1(c4).

In practice, to demix the scattering light-field information easily, we proposed a polarization separation model, as shown in Eq. (3) Ii=Iobj(1−DoLpobj)2+Iobj·DoLpobjcos2(β−i)+Iscat(1−DoLpscat)2+Iscat·DoLpscatcos2(α−i),(3)where the symbols α and β are the angles between the light vector vibration direction of background scattering and object light and the X-axis in the Cartesian coordinate system. Combined with the Stokes vectors, the DoLP of total intensity information is shown in Eq. (4) DoLptotal=PDIobj2+PDIscat2+2PDIobjPDIscat(cos2βcos2α+sin2βsin2α)I,(4)where PDIobj=Imaxobj−Iminobj and PDIscat=Imaxscat−Iminscat mean the intensity differences.

Significantly, the DoLP value of total intensity information is closely related to the light vector vibration direction. Due to the difference in the vibration direction between the object and background scattering information, it is too difficult to get the optimal maximum and minimum intensity images by polarization imaging technique one-time only. However, both the intensity and DoLP information of the object are obviously weaker than the background scattering information, as shown in Fig. 1(d1). More precisely, its values are mainly concentrated in the range of 70% to 80%, which is larger than the target area values of 25% to 30%. [Fig. 1(d3)]. Unfortunately, the intensity value on each pixel is a mixture of the scattering media and the target information. However, the difference between them also reveals the basis for signal demixing. Regarding signal demixing, condition constraints are still lacking. Meanwhile, the angle of polarization (AoP) distribution [Fig. 1(d2)] provides additional constraints. Similar to the DoLP distribution, the AoP results show not only a significant difference between the target and the scattering media but also similar values within the scattering region. Here, combining the information of the DoLP and AoP and using scalar projection theory to apply the blind source separation method37 to the PDI image can lead to separated results in Eq. (5) Iscat=PDI·cos(θn),(5)where θn (matrix) is the difference in light vibration directions on each pixel between the target and scattering media, which can be calculated from the AoP image. In addition, it can be seen from Fig. 1(d4) that the AoP between the target and the scattering media region is close to 90 deg. In other words, by the operation shown in Eq. (5), it is possible to filter out the influence brought by the target information and obtain an intensity distribution image Iscat with the pure scattering media. Based on the same operation, we can obtain an image Isum of the target information under polarized scattering illumination through Eq. (6) Isum=PDI·sin(θn).(6)

A flow chart of the recovery algorithm is shown in Fig. 2. It starts with an initial guess of the sample [Fig. 2(a)]. Figure 2(a2) shows the speckle pattern for sample illumination and the target under speckle illumination. The sample estimate in Fig. 2(a1) is sequentially updated by the polarized measurements taken under different illumination patterns [Fig. 2(a3)]. The updating process is iterated until the solution converges [Fig. 2(a4)]. These steps can be explained in detail as follows.

The Brownian motion of the dynamic scattering media produced different kinds of illuminated patterns, and thus, the corresponding image measurements Ii,n(n=1,2,3…) are acquired. The acquired image is separated into background interference information [scattering media as shown in Fig. 2(b3)] and object information based on polarimetric characteristics, which can be regarded as the scattering illumination pattern and the target in the PIFP method, respectively. The images Ii,n acquired by the detector can be approximated as the result of the product of the object information and the scattering media information after being Fourier transformed and subjected to the optical transfer function (OTF) of the optical system, which can be expressed in Fourier space F(Ii,n)=OTF·F(Iobj·Iscat,n),(7)where F stands for the Fourier transform, and OTF denotes the lens’s optical-transfer-function, and Iobj denotes the initial guess of the ground truth, whereas this initial guess can either be an interpolation of one low-resolution measurement or a random guess.

In our implementation, we used the iterative optimization procedure38 to update the sample image Iobj based on the sequence of acquired images until the estimate converges. It starts with initial guesses (polarization images Ii selected at random) of the object Iobj. This initial guess is then multiplied by the illumination pattern Iscat,n(n=1,2,3……), to produce a target image It,n: It,n=Iobj·Iscat,n. The target image It,n is sequentially updated by the target information under polarized scattering illumination Isum,n, in the Fourier space as follows: F(It,n)u=F(It,n)+OTF·[F(Isum,n)−OTF·F(It,n)].(8)F(It,n)u is then transformed back to the spatial domain to produce an updated image It,nu. The updated image It,nu is then used to update the sample estimate in the spatial domain as follows: Iobju=Iobj+Iscat,n[max(Iscat,n)]2·(It,nu−Iobj·Iscat,n).(9)

The updating process is repeated for the sequences of polarized images, and the updating process is iterated until the solution converges. The process stops if the convergence metric (mean-square error of two consecutively recovered images) is smaller than a pre-defined value.

To demonstrate the imaging performance of the proposed approach, we first used a resolution target (printed on paper) as the object and fresh shallot skin as the dynamic scattering media. A fiber-coupled laser with a wavelength of 650 nm is used for illumination. Then, the laser is collimated and expanded by a pinhole and a convex lens (diameter and focal length, 50.8 mm). A telephoto lens (75-mm focal length, F/# from 1.2 to 12) and an 8-bit industrial camera (5.5-μm pixel pitch, 2048pixel×2048pixel) coupled with a polarizer are used in the imaging system for recording raw images. The F/# of the imaging lens is set to 2.8. The distances from the sample to the telephoto lens principal planes and principal planes to the sensor plane are 550 and 66 mm.

An experimental result demonstration of high-resolution imaging using PFPI is presented in Fig. 2(b). Although the raw intensity image of the United States Air Force (USAF) target through the dynamic scattering media is low-contrast, high-noisy, and seemingly information-less, they contain rich high-frequency information. Figure 2(b1) is a close-up of the raw data, and the minimum feature that can be resolved is group 1, element 2 (222.7μm). To achieve a high-resolution image, 25 sets of orthogonal polarization image pairs (acquired a total of 100 frame images) are captured with a 2-s time interval, which brings the different speckle illuminated pattern images in visualization. The reconstruction results with our PFPI system are shown in Fig. 2(b6). The minimum resolvability limit reaches group 3, element 2 (55.6μm). We note that the resolution of the reconstruction image is improved by 4.005 times that of the original raw intensity image. In addition, the polarization characteristics provide perfect priori information about the dynamic scattering media without being calibrated. Because the high-frequency components relate to the scattering media feature, the PFPI method not only removes the scattering effect but also recovers the detailed information [Figs. 2(b4) and 2(b6)]. We take the resolution line pair of group 2, element 6 (140.3μm). The detailed texture information, which is lost by the scattering effect, is recovered, as shown by the dashed line in Figs. 2(c1) and 2(c2). The optimization recovers a nearly perfect rectangular wave shape in Fig. 2(c2) (group 0, element 2). The high contrast between the lines clearly illustrates that the recovered results are obviously improved.

As another evaluation criterion, the edge spread function (ESF) is utilized to evaluate the reconstructed results. Specifically, by projecting the points in the edge area to a straight line of jump boundary edge, we can drive a discrete ESF, which is not uniformly distributed. Then, the Fermi function is used for the ESF fitting, which is robust to the random noise. The ESF curve obtained from Fig. 2(d2) makes a convincing manifestation. It is proved a fact that the original raw intensity image suffers from poor resolution and strong scattering, whereas the processed image is enhanced substantially in both resolution and scattering suppression. As stated earlier, a random illuminated pattern introduces high-frequency information to improve the resolution of the recovered image. The spectrum distribution information in Fig. 2(d2) is clearly improved than that in Fig. 2(d1), as shown in Figs. 2(e1) and 2(e2). Hence, the proposed platform demonstrates a desired technique to see through dynamic scattering media with a high resolution.

As an attempt to improve the resolution and contrast of the proposed method in imaging through turbid water (skimmed milk added to water in certain proportions39), the underwater imaging setup in Fig. 3 is designed, which submerges the USAF as a target in the turbid water. With the turbid water placed in the optical path, a random optical scattering field is created, which brings the different high-frequency information of the object into the camera. First, the DoP characteristics of the object and scattering media also remain a relatively obvious difference. Moreover, the captured raw-intensity images have presented low contrast and resolution. With the turbidity increasing, the missing detail information is growing, as shown in the first column of Figs. 3(b)–3(d). By comparison, noticeable USAF information is observed as high contrast and resolution in the reconstructed images processed by our method, as we expected. Taking the example of 154 NTU, the zoom-ins of the region of interests (ROIs) ① and ② of Fig. 3(d) clearly show the line pair of group 2, element 6 (140.3μm), which is blurred together in the raw images. After calculation, the resolution of the reconstruction image is improved by 1.41 times that of the original raw intensity image, which can only distinguish group 3, element 3 (99.2μm). Generally, the proposed method makes full use of the characteristics of dynamic scattering media to improve imaging quality in turbid water.

Next, we demonstrate the resolution improvement by imaging different complex objects (paper and metallic material) placed behind the fresh eggshell membranes and ground glass as dynamic and static scattering media separately. For the static scattering media, significantly, the ground glass is moved perpendicular to the axis to generate the different pattern illumination as Brownian motion. To achieve a high-resolution image, the ground glass was moved 7×7 times at 100-μm intervals each time to acquire 98 images (49 Imin images and 49 Imax images). Figures 4(a1)–4(d1) show the captured raw images through scattering media, which are low-contrast, noisy, and seemingly blurry. Similar to Fig. 2, the complex objects are best resolved with our approach [Figs. 4(a2)–4(d2)], whereas the line profiles of the partial region show the perfect performance with some detailed texture [Figs. 4(a3)–4(d3)]. Our method for high-resolution imaging through scattering media extends the application of traditional PIFP imaging methods to scattering imaging.

To ensure the reliability of the proposed technique, we use chicken breast tissue as another complex dynamic scattering media and ground glass with different grit as static scattering media for further analysis. The different acquisition time interval of sub-polarized image pairs is adopted to demonstrate the effect of Brownian motion on reconstruction quality, as shown in Fig. 5. Without loss of generality, we test three group experiments (i.e., 2, 5, and 10 sub-polarized image pairs) with six-time intervals (i.e., 1, 2, 4, 6, 10, and 30 s) in Fig. 5(c), and one group experiment utilizing 50 sub-polarized image pairs with four time intervals (i.e., 1, 2, 4, and 6 s) in Fig. 5(b). From the horizontal comparison, the high-resolution reconstructed images are pretty good and unaffected by the time intervals. This is because the illuminated pattern caused by the Brownian motion of the dynamic scattering media cannot introduce the more high-frequency components at any captured moment. Significantly, different acquisition times cannot increase the number of high-frequency components but only change its representation. The proposed technique is only decided by the shifted high-frequency components information of the passband of the collection optics. However, with the increase in the number of the captured sub-polarized image pairs, the quality of the reconstruction results is obviously improved. This is also highly consistent with the conclusion from the longitudinal comparisons.

To demonstrate the imaging performance of the detail texture, we use the intensity profile across group 2 element 4 to display in Figs. 5(d) and 5(e). As analyzed above, all the curves are clustered together, and the peaks and troughs of the curve have almost the same value at different time intervals. By contrast, each curve presents an obvious difference in the four acquired sub-polarized image pairs. The more image pairs are collected, the better the reconstruction quality.

In addition, to quantitatively evaluate the quality of reconstructions, we use peak signal-to-noise ratio (PSNR), structure similarity index (SSIM), and root mean square error. As evidenced by Fig. 7, the quality of reconstruction improves remarkably when the image pairs increase but improves slightly when the time interval changes.

We employ three kinds of 220, 600, and 1500 grit ground glass as static scattering media to reconstruct the high-resolution object image. First, we generate a complete set of bases illuminated patterns by the axial movement of scattering media. The quality of image reconstruction depends on the number of movements. Here, combined with the complexity of this experiment and the improvement of spatial resolution, we moved the scattering media 50 times to generate the different illuminated patterns. As Figs. 6(a1)–6(c1) shows, the larger the grit of scattering media, the more blurred and lower the signal-to-noise ratio of the acquired images. Similar to the raw intensity image, the performance of the reconstructed image improves with the reduction of the grit of the scattering media. It is because, with the increase of the grit of scattering media, the scattering effect is significantly enhanced. That is, more high-frequency information is transported into the optical system. Hence, better high-resolution imaging can be achieved.

To clearly demonstrate the resolution enhancement using the proposed technique, we extracted a set of spectrum information distributions and the amplitude profile across the middle of the spectrum image in Figs. 6(a3)–6(c3). The increment of the spectrum information is evident both in the two-dimension distribution and the one-dimension curve. The resolution enhancement effect is essentially the same as the dynamic scattering media. Finally, we use the intensity profile across group 0, element 5 and group 1, element 5 of Figs. 6(a2)–6(c2) to demonstrate the details processed by the proposed technique, as shown in Figs. 6(d1)–6(d2). The high contrast between the lines clearly illustrates that the recovered results are obviously improved. We can resolve very well all the smallest features of 600 and 1500 grit ground glass as static scattering media. Compared with the raw intensity image, the reconstructed results also demonstrate the fact that the proposed method is useful for the improvement of spatial resolution. Although the intensity profile across group 0, element 5 and group 1, element 5 does not show significant peaks due to the smooth effect emerged by the small scattering particle, it also displays the small intensity changes in the fringe area.

It is worth noting that the proposed method i utilizes speckle illumination generated by a scattering medium to achieve spectral shifting of the target. With an increased number of captured low-resolution images, the spectral coverage expands, providing more comprehensive detail information. We conducted experiments using 600-grit ground glass as the scattering medium, with a displacement of 100 as the movement distance of the ground glass for each shift. As illustrated in the figures, we compared the raw image [Fig. 7(a1)] with the reconstruction results obtained from different numbers of low-resolution images [5, 20, 50, 100, and 225, shown in Figs. 7(a2)–7(a6)], along with magnified details of the central region [Figs. 7(b1)–7(b6)]. The reconstruction results indicate that an increase in the number of low-resolution images input into the algorithm yields a more pronounced enhancement in resolution. Figure 7(c) presents the comparison of stripe patterns for two sets of element 1 line pairs, which further validates this finding.

In conclusion, we have demonstrated the proposed method for high-resolution imaging through both dynamic (turbid water and fresh eggshell membranes) and static (ground glass) scattering media. The camera captures multiple polarization images for separating scattering media based on the common-mode rejection of polarization characteristics, and then, the separated scattering media images are used to reconstruct the high-resolution image of the target. Our experimental results show that the proposed method enhances resolution more than four times with considerably more detail through scattering media. Similar to the PIFP techniques, our resolution is dictated by the number of sub-polarized image pairs and the roughness for dynamic and static scattering media separately. Our approach could obtain high-resolution images of complex targets without priori information. The proposed approach may provide theoretical support for underwater rescue, aquaculture, non-invasive blood flow, and superficial tumors.

Meng Xiang is currently serving as an associate professor at the School of Optoelectronic Engineering, Xidian University. She obtained her PhD from Xi’an Institute of Optics and Precision Mechanics, Chinese Academy of Sciences. Her research primarily focuses on computational imaging.

Xue Dong is currently serving as an associate professor at the School of Optoelectronic Engineering, Xidian University. She obtained her PhD from Northwestern Polytechnical University. Her research primarily focuses on spectral imaging and spectral polarization imaging.

Fei Liu is currently serving as a professor at the School of Optoelectronic Engineering, Xidian University. He obtained his PhD from Xidian University in 2016. His research primarily focuses on computational imaging and polarization imaging.

Xiaopeng Shao is currently a professor at the Xi’an Institute of Optics and Precision Mechanics, Chinese Academy of Sciences. He received his PhD from Xidian University. His research focuses on computational optical imaging and advanced optoelectronic instruments, etc.

Biographies of the other authors are not available.