High-quality optical imaging through scattering media is crucial for a wide range of applications in fields such as spatial remote sensing, biomedical imaging, and underwater observation^[1–3]. However, the presence of multiple scattering within the scattering medium results in random alterations to the propagation direction, causing a notable shift in the wavefront phase. This phenomenon poses challenges to direct visualization of the target. Consequently, imaging through scattering environments remains a notable obstacle. To address this issue, researchers have developed various strategies in recent years. These approaches include wavefront shaping techniques (WSTs)^[4–9], point spread function (PSF)-based deconvolution^[10–14], non-invasive scatter autocorrelation techniques (SATs)^[15–17], optical phase conjugation^[18–20], and transmission matrix measurements^[21–24]. Among these methods, WSTs and transmission matrix measurements are time-consuming due to the extensive measurement sequences required, making them impractical for dynamic scenarios. PSF, or deconvolution technology, has emerged as a popular method for imaging through scattering media. However, its application requires the prior acquisition of the two speckle intensity patterns beforehand, which needs a difficult experimental process and an excellent speckle collecting precision. A notable advancement in this field is speckle correlation imaging (SCI) introduced by Bertolotti in 2012, which utilized optical memory effects (OMEs)^[25]. This technique has gained recognition for its innovative concept, straightforward optical system design, and capability for single-frame non-invasive imaging. The experimental procedure of SCI involves extracting hidden target Fourier amplitude information from single-shot speckle intensity patterns, followed by iterative phase retrieval algorithms to recover the target Fourier phase spectrum. The Fourier phase and amplitude are then combined and subjected to inverse Fourier transformation to achieve single-shot non-invasive speckle autocorrelation imaging. Subsequently, Katz et al. demonstrated that speckle autocorrelation is equivalent to object autocorrelation, thus propelling this work to new heights^[16]. This groundbreaking study represents a “non-invasive” real-time imaging method that sets it apart from earlier scatter imaging research.

The hybrid input–output (HIO) algorithm, introduced by Fienup in 1982^[26], presented an innovative update strategy within the iterative framework, which enhances the performance of phase retrieval. Although SCI technology combined with the HIO algorithm can achieve real-time imaging, it faces challenges in practical applications due to convergence issues in phase retrieval algorithms caused by random initial stages or over-selection in the initial stage. These challenges are common in phase retrieval algorithms, emphasizing the crucial role of initial phase selection. A recursive algorithm developed by Bartelt et al. in 1984 offered a solution for recovering the object phase from the bispectrum^[27], which was the Fourier transformation of the triple correlation^[28]. Recent work by Wu et al. demonstrated object reconstruction through inverse Fourier transformation without iterative phase algorithms^[29]. In recent years, the distinctive phase reconstruction technology of the bispectrum has garnered increasing attention from researchers in the field of speckle imaging^[30,31]. However, the quality of reconstruction in single-shot imaging using the bispectrum phase and Fourier amplitude is limited. The accuracy of reconstruction information from single-shot imaging through the scattering medium still deserves further improvement and exploration.

In this Letter, we introduce an innovative strategy that merges the advantages of bispectral analysis and iterative methodologies to enhance the fidelity of reconstructed images. Our approach commences with the reconstruction outcomes derived from the bispectral technique. Specifically, adopting the bispectral phase to replace the random phase in the iterative algorithm, the iterative phase process necessitates only a simple iterative step to obtain a high-quality reconstructed image. Comparative analysis with extant research indicates that our method significantly surpasses conventional phase retrieval algorithms in terms of imaging reconstruction quality. Furthermore, we have noted that the dimensions of speckle autocorrelation significantly impact the reconstruction process, with the precision of reconstructed images contingent upon the size of speckle autocorrelation. To enhance retrieval efficiency and diminish computational burden, we propose a robust technique for assessing imaging resolution utilizing an iterative search function to ascertain the optimal speckle autocorrelation size. Multiple image sharpness evaluation metrics are employed to identify the ideal autocorrelation size for the reconstructed image. The application of this method not only enhances the efficiency of image reconstruction but also facilitates the achievement of high-quality imaging reconstruction within a single iteration. Through experimental validation, we demonstrate the superiority of our proposed methodologies over existing phase retrieval techniques, offering a more precise and efficient algorithm for single-frame speckle imaging.

The OME is a significant phenomenon observed in scattering systems, first identified by Freund et al. in 1988, which has since served as a foundational concept in subsequent research in the field of SCI. The OME posits that the speckle pattern distribution undergoes a simple translation when the incident angle of light on the scattering medium changes within a specific field of view. Essentially, the scattering system can be conceptualized as an optical system that displays linear invariant displacement within the parameters defined by the OME. Consequently, the process of scattering imaging can be understood as a linear representation of the convolution between the PSF and the target object. The relationship between the intensity of speckle patterns and the target object can be mathematically expressed as I⊗I=(O*S)⊗(O*S)=(O⊗O)*(S⊗S)≈O⊗O+C,(1)where the symbol “⊗” represents the correlation operator and “*” denotes the convolution operator. The variable “I” signifies the speckle intensity pattern recorded by the camera, while “O” represents the obscured target object behind the scattering medium. “S” denotes the PSF of the imaging system. “C” denotes a constant. The signal spectrum is equal to the Fourier transform of its autocorrelation function, as stated by the Wiener–Khinchin theorem. In the case of a test object, the Fourier spectral intensity can be determined accordingly, F{I⊗I}≈F{O⊗O}=|F(I)|2.(2)

The symbol “F{.}” is presented as a Fourier-transform operator. Here, the Fourier spectrum intensity characteristics of an object can be obtained through the calculation of the autocorrelation of the captured speckle patterns. Subsequently, the iterative phase retrieval method in the Fourier domain can be utilized to reconstruct the concealed target object through iterative processes. However, the Fourier phase of the object remains unknown, requiring a random phase as the initial phase for the iterative reconstruction. This approach is hindered by the drawbacks associated with employing a random phase as the initial phase. Indeed, the direct retrieval of the phase from speckle patterns has consistently posed a significant challenge. Fortunately, Bartelt et al. indicated that the Fourier phase of the speckle pattern corresponds to the bispectrum phase, presenting a novel approach to overcoming this challenge. The formula for the speckle bispectrum is outlined as follows: BI(u1,u2)=I˜(u1)I˜(u2)I˜(−u1−u2).(3)

The symbol “B” denotes the speckle bispectrum, “∼” represents the Fourier transform, and the variable u represents the two-dimensional spatial frequency vector. According to the OME, Eq. (3) can be expressed as $$\begin{gathered} BI(u1,u2)=BO(u1,u2)·BS(u1,u2)=O˜(u1)O˜(u2) \\ ·O˜(−u1−u2)·S˜(u1)S˜(u2)S˜(−u1−u2). \end{gathered}$$(4)

The spatial frequency spectrum of the object and the spatial frequency spectrum of the transfer function can be depicted by separating the formula into amplitude and phase components, $$\begin{gathered} BO(u,v)=|O(u)O(v)O(−u−v)| \\ ×exp{j[φo(u)+φo(v)+φo(−u−v)]}, \end{gathered}$$(5)$$\begin{gathered} BS(u,v)=|S(u)S(v)S(−u−v)| \\ ×exp{j[φS(u)+φS(v)+φS(−u−v)]}. \end{gathered}$$(6)

The speckle bispectrum average ensemble is ultimately expressed in the following manner: $$\begin{gathered} ⟨BI(u,v)⟩N=BO(u,v)×BS(u,v)N=|O(u)O(v)O(−u−v)| \\ ×exp{j[φo(u)+φo(v)+φo(−u−v)]} \\ ×⟨|S(u)S(v)S(−u−v)|⟩N, \end{gathered}$$(7)where ⟨…⟩ denotes the collection comprises sub-speckle pattern derived from individual single-frame speckle patterns, N represents the number of sub-speckles, and S⟨…⟩ denotes the PSF. Consistent with the findings of Lohmann et al., the bispectral statistical average of the PSF closely approximates a real number. Consequently, according to Eq. (7), Eq. (8) can be derived. The phase information of the target Fourier transform can be directly obtained from the average bispectrum: exp{⟨BI(u,v)⟩N}≈exp{BO(u,v)}.(8)

When the conclusion of Eq. (8) is substituted into Eq. (3), it results in a regression equation [Eq. (9)]^[32,33], and Eq. (10) can be employed to extract the Fourier phase of the object, $$\begin{gathered} exp{⟨BI(u,v)⟩N}=exp{I˜(u)}+exp{I˜(v)} \\ +exp{I˜(−u−v)}, \end{gathered}$$(9)ϕ{⟨BI(u,v)⟩N}=φu+φv+φ−u−v,(10)where “ϕ” and “φ” denote the phase operators. We assume the initial condition as φ0, φ1=0^[34], and the recursive algorithm is utilized to obtain phase data for all frequencies, resulting in the comprehensive Fourier phase representation of the object.

Figure 1 shows the schematic diagram of our experimental system. A coherent laser emitting continuous-wave linear polarized light at a wavelength of 532 nm was employed as the light source. The emitted light was magnified by a microscopic objective approximately 10 times, directed through a focusing lens, and then passed through a rapidly rotating diffuser to create an incoherent light source. The target object was positioned at a specified distance U in front of the scattering medium. To ensure that only the target area was illuminated by the incident light, the dimensions of the aperture are adjusted according to the field of view of the object. In our experiment, the optimal aperture size has been determined to be 2 cm. The size of the light spot incident on the scattering medium can be controlled, thereby reducing interference from extraneous noise sources during the reconstruction process. After the transmitted light was scattered by the diffuser (220 grit frosted glass), the speckle image was directly collected by a high-resolution camera (SP928 Beam Profiling Camera, United States, resolution 1928×1448, pixel size 3.2 µm).

The proposed method necessitates the derivation of the Fourier phase from the bispectrum as a crucial step toward achieving high-quality image reconstruction. To mitigate the influence of noise in the speckle intensity patterns, the experimental method employs the principle of time-averaging, dividing the raw speckle image evenly into 525 subimages of 200pixel×200pixel with an overlap ratio of 0.8. The projection method described in Ref. [35] is employed to project the speckle patterns into a series of 1D projections, thereby circumventing substantial computational complexity. As depicted in Fig. 2, the final bispectrum phase (2D Fourier phase) of the object is obtained by arranging all the retrieved 1D Fourier phases. Since the 1D Fourier phase is obtained through angular projection, the bispectrum phase is affected by the angle. To further refine the precision of the bispectrum phase, by normalizing all sub-speckle bispectrums, the average speckle bispectrum is derived and utilized as the initial phase for the iterative algorithm aimed at image reconstruction.

In this work, we employ the phase derived from bispectrum analysis as the starting point for the Fienup iterative algorithm. This approach enhances the convergence rate and quality of reconstruction in the phase-retrieval process. The key procedural steps of our technique are depicted in Fig. 2, encompassing the extraction of the bispectral phase and the integration of constraints into the phase-retrieval algorithm. Initially, a recursive method is used to extract the bispectral phase from single-frame speckle intensity patterns, which is then combined with the Fourier amplitude spectrum obtained from the speckle autocorrelation to form the initial phase for the phase-retrieval algorithm. Throughout the phase retrieval procedure, the HIO-ER algorithms, as introduced by Fienup, are employed to achieve high-quality image recovery of the target object. By preserving the true object orientation and complete object-phase information, the bispectral Fourier phase reduces the likelihood of encountering local minimum issues during the iterative process. High-resolution images can be obtained through just 50 iterations during the process. Nevertheless, in the conventional HIO-ER algorithm, recovering a complex object like “θ” with completely random initial phases necessitates at least 150 iterations. The introduction of the bispectrum significantly enhances both the reconstruction accuracy and the convergence rate of the iterative algorithm.

To further evaluate the superiority of the proposed bispectrum analysis-based phase-retrieval approach, a series of proof-of-concept experiments were conducted. In the experiment, two opaque perforated plates were used with a dark background and light-colored patterns. One plate features a perforated letter “θ.” The length of “θ” is approximately 303 µm, and the width is approximately 190 µm. The other plate features a perforated letter “A.” The length of “A” is approximately 313 µm, and its width is approximately 302 µm. Figure 3 gives a comparison between our technique and other scattering imaging methods, utilizing the experimental setup depicted in Fig. 1. The experiments involved illuminating a target object with a spatially incoherent source, passing through opaque scattering barriers to generate speckle images, which were then captured by a camera. Figure 3 displays two experimental samples, showcasing the raw camera speckle images, images reconstructed using the bispectral inverse Fourier transform method, the HIO iterative phase retrieval algorithm, and the proposed method. It is evident that our technique, as opposed to simply performing a single inverse 2D Fourier transform from the bispectral phase, exhibits insensitivity to background noise and yields high-contrast test targets. Traditional phase-retrieval algorithms typically necessitate multiple initial phase selections for the computational retrieval process, leading to reconstruction outcomes that are contingent on the initial phase chosen. In this study, the optimal results were compared with our method. The HIO algorithm employed for image reconstruction in comparative experiments is inherently influenced by entirely random initial phases. Utilizing only 50 iterations for image reconstruction may yield ambiguous results. To mitigate the uncertainty and highlight the superiority of the proposed method, we chose to compare it with the reconstructed images obtained after 150 iterations in the HIO algorithm. As depicted in Figs. 3(c) and 3(d), both the bispectral iterative recovery algorithm and the iterative phase retrieval method employed 150 iterations. The experimental results indicate that the shapes of objects recovered using only random phase recovery exhibit varying degrees of distortion. Our proposed method significantly enhances imaging quality compared to the random phase approach. Furthermore, the iterative process notably accelerated the image recovery speed. Consequently, our single-shot scattering imaging scheme can achieve high-quality imaging and perform effectively in noisy environments, rendering it particularly suitable for diverse applications such as imaging through opaque scattering barriers and mitigating atmospheric turbulence in astronomy.

In our approach, it is important to highlight that the size of the target autocorrelation significantly influences the quality of reconstructed images. To illustrate, various runs of the phase-retrieval algorithm were conducted using different sizes of target autocorrelation under identical initial conditions, and the corresponding reconstruction outcomes from the experiment are depicted in Fig. 4(a). Ideally, the most accurate reconstruction with the lowest error would be considered optimal. However, due to the challenge of determining the ideal target autocorrelation size during the phase-retrieval process, the reconstructions may not always be optimal. It is important to note that neither small nor large sub-speckle sizes can ensure optical imaging results. For instance, sub-speckle sizes of 230pixel×230pixel and 240pixel×240pixel, despite their seemingly minor size difference, can yield markedly different imaging results. To address this issue, we have implemented an adaptive search strategy that relies on an image sharpness evaluation function to precisely identify the suitable autocorrelation size. Over the past few decades, numerous refocus criteria have been proposed^[36,37]. The SMD function, Brenner function, Laplace function, Roberts function, Vollaths function, and other extensively studied functions are commonly utilized for assessing image clarity. In our method, three classical gradient functions, the discrete Fourier transform (DFT) function^[38], and traditional gradient functions like the Brenner^[39] and Laplace functions^[40] are employed as criteria to determine the target autocorrelation size. The mathematical formulations for each criterion are detailed in Eqs. (11)–(13), respectively, Brenner=∑x∑y{[f(x+2,y)−f(x,y)]2},(11)Laplace=∑x∑yG2(x,y),G(x,y)=f(x,y)⊗L,(12)DFT=1MN∑u=0M−1∑v=0N−1u2+v2p(u,v).(13)

Here, f(x,y) in Eq. (11) represents the gray value at pixel (x,y). Equation (12) uses the Laplacian operator to convolve the image and obtain its high-frequency components. G(x,y) is the convolution of the Laplacian operator at the pixel f(x,y), “⊗” is the convolution symbol, and L represents the Laplacian operator. Equation (13) utilizes a weighting coefficient u2+v2 to represent the distance from the pixel to the central pixel. M and N represent the dimensions of the image; P(u,v) denotes the square of the image spectrum. We employ three image sharpness functions based on different principles to enhance the accuracy of the experimental results.

The normalized values of various criteria exhibit fluctuations with changes in autocorrelation. Under consistent experimental conditions, this study collected nine sets of speckle autocorrelation data with varying sizes for iterative reconstruction. Subsequently, the sharpness of the nine reconstructed results was assessed using a clarity function, as illustrated in Fig. 4(b). The curve as a whole displays varying levels of fluctuations, indicating that the size of the object’s autocorrelation significantly influences the quality of the resultant image. Moreover, the curves are generally unimodal, each displaying a distinct global peak at the optimal target autocorrelation size. Notably, when the object’s autocorrelation size is set to 240×240, the criteria all reach their maximum values, suggesting that this dimension may be optimal for achieving the highest quality image reconstruction. Subsequent experiments validated this observation, showing that imaging results at 240pixel×240pixel outperformed those at 220pixel×220pixel or 230pixel×230pixel, and similarly, imaging results at the minimum of the curve at 270 pixels showed inferior quality. Therefore, by utilizing an image clarity function, we can efficiently search for the most suitable autocorrelation size for imaging, thereby obtaining optimal imaging results.

The complete reconstruction process of the bispectral iterative algorithm, which is based on adaptive image clarity in our imaging approach, has the potential to enhance reconstruction fidelity and convergence. The workflow, illustrated in Fig. 5, involves several steps. Initially, speckle bispectrum intensity patterns are generated through bispectrum analysis. Subsequently, the object’s Fourier phase is obtained by substituting the recursive formula from the preceding section. The determination of the size of the target autocorrelation in the subsequent phase-retrieval algorithm is crucial, and we glean the optimal autocorrelation dimension by discerning the peak value of the image clarity function. At this stage, we have acquired the object’s Fourier magnitude and preliminary phase. To prevent computational inaccuracies, the Fourier phase size is constrained using predefined optimal dimensions. Once the phase and magnitude are established based on the ideal size, they are incorporated into the HIO iterative method to reconstruct the entire object. The bispectral Fourier phase demonstrates a more precise directional property compared to random phase counterparts, thereby significantly reducing the likelihood of getting stuck in local minima during the iterative process. This, in turn, accelerates the convergence of image reconstruction, substantially enhancing the robustness of image reconstruction.

The comprehensive flow chart outlines the advantages and limitations of our experiment. The advantages include a clear initial phase to prevent local iterative retrieval issues, a straightforward iterative process that eliminates the need for repeated initial phase selection, and the ability to retrieve the most suitable image size for reconstruction using the adaptive image clarity function to reduce noise impact on reconstruction results. This method demonstrates reliable generalization capabilities across various imaging scenarios and maintains robustness throughout the imaging process. By utilizing the image clarity function to identify the most suitable reconstruction size, this approach aims to overcome the problem of suboptimal imaging results caused by miscalculation of the object size in traditional iterative calculations. However, the experiment also has limitations such as the time required to obtain the bispectrum phase and select an appropriate autocorrelation size. In addition, it necessitates the acquisition of a high-precision speckle pattern and the constraint of the imaging range by the OME, limiting its applicability to small targets at present. This may lead to an expansion in the operational processes of the system, which is deemed acceptable if it enables the acquisition of higher-quality reconstructed images.

In conclusion, we have developed a unique framework that combines the benefits of bispectral analysis and the phase-retrieval algorithm, for exploiting the OME imaging in scattering scenes. Through our experimental endeavors, we have successfully showcased the accurate reconstruction of high-quality images using our method, enabling non-invasive single-shot imaging through scattering media. Furthermore, we have introduced a technique for determining the optimal size of the target autocorrelation with precision by employing an image clarity evaluation function. By utilizing three established evaluation criteria, our method facilitates the precise identification of the ideal autocorrelation size for the iterative reconstruction process, thereby streamlining computational efforts and ensuring that only essential information is retained, thus mitigating interference from extraneous noise and enhancing image quality. Additionally, we have conducted a quantitative comparison of reconstruction outcomes between our approach and alternative scattering imaging techniques, providing insights into the superior performance of our proposed method. Our approach holds promise as a valuable alternative in diverse scenarios necessitating high-contrast speckle imaging.