1. INTRODUCTION

Although performing imaging with scattered light is intractable, it has wide applications in many fields, such as biomedical, astronomical, and underwater imaging. This problem is challenging in that the incident wave is randomly modulated by the scattering media, resulting in a disordered pattern containing massive speckle spots. Fortunately, thanks to recent advances in wave front shaping [1–5], the disordered information can be redistributed and decoded from the observed speckle pattern. For instance, measuring the transmission matrix (TM) allows one to perform focusing [6–9], as well as imaging [10–12] or transmitting information [13–15], through or inside the scattering medium. While, in principle, acquiring TM requires coherent illumination, its extension to broadband case will open fascinating applications in ultrafast science. For this purpose, the TM-based methods are further extended in the spectral domain by accessing the multispectral transmission matrix (MSTM) to cope with polychromatic light. Such techniques have shown success in ultrafast lasers [16–19] and fiber optics [20]. However, characterizing the scattering medium requires thousands of images captured under vibration isolation conditions, which precludes its applications in dynamic and harsh environments.

A separate approach known as speckle correlation imaging (SCI) also shows promise in broadband focusing and imaging through scattering media [21,22]. This approach relies on persistent correlations, such as the optical memory effect (OME) [23,24], which enables single-shot, high-fidelity, diffraction-limited imaging through scattering media without any prior knowledge of or access to the scattering media. While promising, regular SCI, based on the speckle generated from quasi-monochromatic illumination (typically $<20\mathrm{nm}$), is inherently incompatible with the extremely broad nature of white light spectra. However, the crucial coherence length $L_c={\lambda }^2/\mathrm{\Delta }\lambda$ of the light sources is inversely proportional to the illumination bandwidth $\mathrm{\Delta }\lambda$ [25]. Speckle patterns emitted by broadband radiation display low-contrast structures and thus obstruct optical manipulation and inverse reconstruction.

Indeed, the broadband pattern can be thought of as the incoherent superposition of multiple monochromatic speckle patterns produced by the different frequency components encompassed by the light spectrum. Essentially, when observing a speckle pattern, an increase in wavelength (all other parameters remaining fixed) has two effects. First, the observed speckle pattern undergoes a spatial scaling. When the wavelength increases, the speckle pattern expands and vice versa. The second effect is on the random phase shifts imparted to the scattered wave by the surface-height fluctuations $h$. Since those phase shifts are proportional to $h/\lambda$, an increase in wavelength will decrease the random phase shifts and vice versa. When the surface is rough on the scale of a wavelength, phase wrapping into the interval (0, $2\pi$) will result in changes in the detailed structure of the speckle pattern [26]. For the above two reasons, the extended illumination bandwidth will result in dephasing and speckle aliasing [27], thus reducing the contrast of the observed speckle pattern and further obstructing the reconstruction process.

As different spectral components generate different speckles, current approaches are highly sensitive to the spectral bandwidth. The monochromatic treatment is, therefore, only valid as long as the bandwidth of the source is contained within this spectral correlation bandwidth. Extending narrowband to broadband light can effectively ease the trade-off between the filter bandwidth and detection signal-to-noise ratio (SNR). Furthermore, it can potentially help increase the penetration depth, which is critical for biomedical fields, including neuroscience.

To mitigate the constraint of the narrowband operation, several imaging techniques have been developed. A straightforward method is based on spectral components such as spectral filters [28], prisms [29], or even scattering-assisted spectrometers [30] to resolve broadband spectra into well-separated spectral channels. This improves the speckle contrast as well as the performance of object retrieval. Although these methods have proven efficiency, the trade-off between spectral resolution and spectral range can limit their performance. Moreover, imposing spectral filters on broadband sources leads to a dramatic reduction in flux, which may be impractical in some applications. The ability to use the full-flux of low-flux ultra-broadband sources, such as table-top high-harmonic generation (HHG) systems [31], could provide a breakthrough for practical imaging applications. Incoherent scattering imaging can also rely on deep learning modules [32]. Nevertheless, the price to pay is a fuzzy physical model and a large training dataset. To date, controlling scattered light for deep focusing under broadband illumination has attracted much attention [16–19,33]. However, the necessary scattering imaging methods combining broadband illumination with diffraction-limited spatial resolution are currently lacking. These representative state-of-the-art scattering imaging methods, as well as the effective illumination bandwidth, are summarized in Table 1.

Table 1. Comparison among State-of-the-art Scattering Methods^a

View Table | View all tables in this article

In this paper, we propose two individual methods for realizing scattering imaging under broadband illumination using only a single speckle measurement under noninvasive and invasive conditions. The major difference between invasive and noninvasive methods is that invasive methods pre-calibrate point spread function (PSF) information for scattering systems, whereas noninvasive methods require no a priori information about the medium or its properties. In contrast to previous approaches, both methods neither require any spectral filters, nor have additional assumptions on the spectral content of the signals. In the noninvasive approach (method A), we extract the object Fourier spectrum from the broadband pattern by using a cepstrum shaping method combined with a matrix decomposition strategy. We show that by unique computational processing, the lost low-frequency information can be reproduced in the Fourier domain. After applying a modified phase retrieval procedure, sharp reconstruction results can be achieved even with a speckle contrast ratio (SCR) as low as 0.0823 [26]. The invasive approach (method B) is based on an optimization framework to solve a general de-correlagraphy problem. The proposed computational framework can effectively reduce the background level of the recovered image, while preserving sharp edges and high-frequency features. Experimental validations under the 280 nm illumination spectrum spanning the visible region show the applicability of both methods. Our approach enables imaging through scattering media under broadband radiation, which we expect to make much more efficient use of such state-of-the-art light sources as third-generation synchrotrons and table-top HHG sources for operating in heavily scattering environments. Moreover, the capability to treat broadband light can also empower the fields of florescence imaging, night vision detection, and passive sensing.

2. CORRELATION ANALYSIS

For simplicity, our problem discussion is limited to the scope of the OME region. The recorded broadband speckle $I_B$ can be described by the convolution between the broadband PSF, ${\mathrm{PSF}}_B$, and the object function, $O$,

Let us assume that the broadband spectrum light source $\mathcal{Q}(\lambda ,\mathrm{\Delta }\lambda )$ can be regarded as a linear superposition of multiple narrowband spectrum light sources ${\mathcal{Q}}_i({\lambda }_i,\mathrm{\Delta }{\lambda }_i)$ $(i=1,\mathrm{2,}\dots ,M)$ with each one under the narrowband illumination hypothesis. Here, $\lambda$ and ${\lambda }_i$ represent the central wavelength of the broadband and each narrowband light, and the bandwidths of the broadband and narrowband light are denoted by $\mathrm{\Delta }\lambda$ and $\mathrm{\Delta }{\lambda }_i$, respectively. Considering the broadband illumination case, the broadband PSF, ${\mathrm{PSF}}_B(\lambda ,\mathrm{\Delta }\lambda )$ can be expressed as

From Eq. (4), the autocorrelation of the broadband speckle can be separated into two parts. The first term incorporates a superposition of all autocorrelation parts, whereas the latter term accounts for contributions from all cross-correlation parts. Considering the discussions in Ref. [38] and supposing the imaging system has a circle pupil, the first part of the broadband PSF’s autocorrelation can be further expressed as

3. METHODS

A. Imaging without PSF Calibration

A radiation beam consisting of multiple wavelengths is illuminated on the sample and the exiting field is scattered by an off-the-shelf diffuser to form broadband pseudorandom patterns, which is the incoherent superposition of speckle patterns of each wavelength [see Figs. 1(a) and 1(b)]. To reconstruct the object, $O$, from the single-shot measured image $I_B$, we need to solve the ill-conditioned inverse scattering problem by utilizing phase retrieval. However, as discussed in Section 2, the autocorrelation of the broadband ${\mathrm{PSF}}_B$ no longer satisfied the delta function approximation and thus caused the failure of the SCI method. Novel methods are highly desired to remedy this issue. Here, based on experimentally observed similarities in the Fourier spectra of objects and speckles and inspired by homomorphic filtering, we propose a spectrum refinement principle that enables the extraction of the pure object Fourier spectrum from the caustic broadband pattern Fourier spectrum. The Fourier spectrum extraction procedure is depicted in Fig. 2. The phase retrieval is then performed using the refined Fourier spectrum pattern [see Fig. 1(d)] as input into a hybrid relaxed averaged alternating reflections (RAAR) [39] and hybrid input and output (HIO) [40,41] algorithm, with a shrink-wrap-like support update [42] (see Algorithm 1 for details). The object, support region, and Fourier phase can be retrieved simultaneously [Fig. 1(e)].

 

Fig. 1. Principle of single-shot broadband scattering imaging. In the case of a broadband source (a), the broadband pattern (b) is the incoherent spectrally weighted sum of the monochromatic speckle patterns corresponding to all wavelengths present in the source. The size of these monochromatic speckles is geometrically scaled with increasing wavelengths, but the micro-structures are ever-changing. In the presented method, the broadband speckle is first transformed into the Fourier domain (c) and then refined in the complex cepstrum to extract an estimated object Fourier spectrum (d). (e) A modified iterative phase retrieval algorithm is then used to reconstruct the sample, support, and Fourier phase simultaneously.

Download Full Size | PDF

 

Fig. 2. Noninvasive broadband scattering imaging reconstruction pipeline. The broadband pattern $I_B$ (a) is first filtered by a two-dimensional Hanning window to avoid spectrum leakage (b) $W·I_B$ and then transformed into the Fourier domain. The Fourier amplitude (c) is reserved and converted into the logarithmic domain (d). The spectrum is filtered by NLM (e) and transmitted back into the Fourier domain (f). A crucial step to remove the peak noise is handed over to the LRSD. After eliminating the sparse part (g), the reserved low rank part (h) is certified as the estimated object spectrum. A modified phase retrieval is then applied on (h) to retrieve the object information. The ground truth object Fourier spectrum is presented in (i). The line profiles (j) along the opposite yellow arrows in (i) show a comparison among the original (c), the ground truth (i), and the recovered spectrum (h). Notably, the recovered spectrum contains richer low-frequency information than the original spectrum. The windowed logarithmic transformed PSF spectrum is shown in (k) with the peak highlighted by a red arrow. NLM, nonlocal means; LRSD, low rank and sparse decomposition.

Download Full Size | PDF

After we measure the broadband pattern $I_B$ [Fig. 2(a)], we first filter it with a two-dimensional Hanning window in real space $W·I_B$ [Fig. 2(b)] to avoid spectrum leakage and then transform it into the Fourier domain [see Fig. 2(c) for Fourier amplitude]. The part of the Fourier phase corrupted by scattering is discarded, which is not shown here. We note, first, that the Fourier spectrum is sparse, as the energy of the power spectrum is concentrated around zero frequency [see the blue line profile in Fig. 2(j)], which necessitates an approach to retrieval. To avoid oscillations caused by near zero values in the spectrum, we add an all-one matrix on the Fourier amplitude and then convert it into the logarithmic domain [Fig. 2(d)]. Notably, more details and noise in the spectrum are highlighted with logarithmic transform; therefore, a smooth filter via nonlocal means (NLM) [43] is applied to exclude the noisy profile introduced by the broadband ${\mathrm{PSF}}_B$ [see the hemispherical surface attached burrs in Fig. 2(k)]. The size of the patch and research window of the NLM are set as 5 and 11, respectively, and the smooth parameter $\gamma$ is set as 3. As $\gamma$ controls the decay process of the Gaussian function in NLM, a higher value of $\gamma$ is suitable for a relatively simple target to obtain a smoother spectrum and vice versa. Detailed parameter selection rules are arranged in Appendix C.

The filtered spectrum is then transformed back into the Fourier domain [Fig. 2(f)]. The next step involves a low-rank and sparse decomposition (LRSD) procedure [44] based on robust principal component analysis (RPCA) [45] to remove the residual peak noise caused by the broadband ${\mathrm{PSF}}_B$ [highlighted with a red arrow in Fig. 2(k)]. We classify this peak noise into a sparse part according to its sparse property [Fig. 2(g)], and the low-rank part represents the final separated Fourier spectrum [Fig. 2(h)]. The ground truth object Fourier spectrum is presented in Fig. 2(i). Note that the low-frequency information reappears in the Fourier domain after processing [see the comparison in Figs. 2(c), 2(h), and 2(i)]. The central line along the two opposite arrows of Figs. 2(c), 2(h), and 2(i) depicted by gray, red, and blue lines also verified the effectiveness of the proposed method, for richer low-frequency information and higher cut-off frequency.

Finally, we obtain the refined Fourier spectrum $|\widehat{\mathcal{F}}\{O\}|$ [Fig. 2(h)], and in the next step, we perform phase retrieval to reconstruct the target information. However, as the broadband ${\mathrm{PSF}}_B$ information is unknown, object spectrum recovery is an ill-conditioned problem; therefore, robust phase retrieval is needed to remedy this issue. Here, we modified the iterative alternating projection method by introducing an adaptive support region and a hybrid phase retrieval strategy. Specifically, the iteration was carried out using a combination of the RAAR and HIO algorithms, with a shrink-wrap-like support update. We calculate the root mean square error (RMSE) between the refined Fourier amplitude and the Fourier amplitude generated by the estimated object in each iteration. We chose typically 80 iterations to ensure that the RMSE of the RAAR algorithm converged to $<0.05$. The subsequent HIO algorithm is used for smearing out the background noise and 10 iterations are typically adequate to obtain the final result. The feedback parameter $\beta$ is set as ${\beta }_1={\beta }_2=0.8$ for both the RAAR and HIO algorithms. Accompanied by the iteration, the adaptive support region is controlled by an alternative convolution with a Gaussian kernel (Algorithm 1, step 6) and a threshold operation (Algorithm 1, step 7), which act as dilation and erosion. As the iterative number increases, this exported support region will change from a loose size to a tight one (Algorithm 1, steps 2-7). The initial ${\sigma }^{(0)}=3ϵ$, where $ϵ$ represents the pixel size related to the imaging sensor. It is worth emphasizing that we selected a non-centrosymmetric support $S_0$, for instance, an equilateral triangle, as an initial support. Using a non-centrosymmetric support region considerably speeds up convergence and effectively avoids the twin-image problem.

To explain how the process works in greater detail, four projection operations are defined. The projection operations $P_s$, $P_m$, $R_s$, and $R_m$ project points in two sets, $S$ (support) and $M$ (modulus) [46]. When the image belongs to both sets simultaneously, we reach an equivalent solution. The support projector $P_s$ involves setting to 0 the components outside the support, while leaving the rest of the values unchanged.

The reflectors of the support projector ${\mathbf{P}}_s$ and modulus projector ${\mathbf{P}}_m$ are defined as, ${\mathbf{R}}_s=2{\mathbf{P}}_s-\mathbf{I}$ and ${\mathbf{R}}_m=2{\mathbf{P}}_m-\mathbf{I}$, respectively.

According to the notation above, the recursive form of HIO algorithm would be

For clarity, we describe our modified phase retrieval procedure in Algorithm 1.

Algorithm 1. Modified phase retrieval algorithm

View Table | View all tables in this article

In the present method, phase retrieval is a standalone step that depends solely on the refined spectrum. It can also be integrated as a module in algorithms for other types of lensless imaging (coherent diffraction imaging, holography, ptychography) to allow them to cope with broadband sources.

B. Imaging with PSF Calibration

Figure 3 presents the data acquisition and reconstruction process with PSF calibration. An object and a pinhole are alternatively placed at the same spatial position in the optical configuration (see details in Appendix C), which captures the broadband speckle, $I_B$, and the corresponding PSF, ${\mathrm{PSF}}_B$. The raw captured speckle and PSF are then lowpass filtered in the frequency domain (a typical filter length is 15) to eliminate the non-uniform illumination background. Then, we apply a two-dimensional Gaussian lowpass filter of size [2 2] with standard deviation $\sigma =0.5$ in the real space to further smooth the pattern.

 

Fig. 3. Broadband scattering imaging reconstruction pipeline with calibrated PSF. In the OME range, the broadband speckle $I_B$ can be treated as the convolution of the object $O$ and the corresponding broadband ${\mathrm{PSF}}_B$. Data acquisition (a) requires a one-time calibration procedure to measure ${\mathrm{PSF}}_B$. Raw captured speckle and the corresponding PSF are background homogenized and filtered by a Gaussian low-pass filter (b) in the frequency domain [(c) and (d)]. The image reconstruction algorithms (e) include cross correlation, nonlinear reconstruction, and our computational reconstruction method. The retrieved four-leaf clover patterns are compared in (f). Artificial cross-section lines are highlighted by gray dashed lines in the nonlinear reconstruction. The normalized intensity of the 76th column along the yellow dotted line in each result (f) is presented in (g). Notably, our method outperforms the other methods, achieving lower background and higher fidelity. GT, ground truth; CC, cross-correlation; NL, nonlinear reconstruction.

Download Full Size | PDF

After preprocessing, we introduce two types of reconstruction algorithms. The first algorithm is based on the cross-correlation strategy:

Here ${\mathcal{F}}^{-1}\{·\}$ represents the inverse Fourier transform, and the bar above the variable represents a complex conjugate. Inspired by the interferenceless coded aperture correlation holography [47], a modified approach of the cross-correlation method is referred to as the nonlinear reconstruction process:

Mathematically, broadband scattering imaging can also be formulated as the following forward model:

We develop an iterative inverse algorithm, modified from sparse Poisson intensity reconstruction algorithms theory and practice (SPIRAL-TAP) [49], to obtain an approximate solution of the inverse problem in Eq. (14). Here, in our case, the forward model matrix $\mathbf{H}$ is defined as a function called $\mathbf{H}(·)$ and the adjoint of $\mathbf{H}(·)$ is specified as ${\mathbf{H}}^T(·)$, which denotes the back-projection process. Contrary to the previous deconvolution approaches, the forward and backward models are modified with the cross-correlation [Eq. (11)] and nonlinear reconstruction operation [Eq. (12)], respectively. The adjustable parameters are set the same as in nonlinear reconstruction ($\alpha =0.4$ and $\beta =0.3$). The SPIRAL-TAP approach for the de-correlagraphy problem will consider the iterations to be sufficiently close to the optimal solution if

Some artifacts from the recovery process nevertheless deserve attention; for instance, the abnormal cross-section lines indicated by gray dashed lines in the nonlinear reconstruction [Fig. 3(f)]. These cross-section lines come from the noninteger power transformation of the spectrum in the nonlinear reconstruction, resulting in a simultaneous increase in detail and high-frequency noise, unavoidably producing these artifacts. As these adjustable parameters $\alpha$ and $\beta$ decrease, the resolution of the reconstruction will improve moderately, while the noise will increase dramatically. (Details about this point can be found in Appendix C.) However, thanks to the regularized inversion method, these artifacts are not dominant in the recovered pattern [see the comparison in Fig. 3(f)]. In addition, the differences between the values of points a and b are 1, 0.6, 0.68, and 0.75, respectively [Fig. 3(f)]. Reconstructions of the proposed method have sharp boundaries and contain very little noise. A typical reconstruction requires at least four iterations. Solving for $1800\times 1800=3.24$ million pixels takes 10.35 s (2.59 s per iteration) on a dual-core computer and with a GeForce RTX 2080 Ti GPU. Moreover, for further algorithm acceleration, this fused lasso problem can also be solved by the alternating direction method of multipliers (ADMM) [50] or fast iterative shrinkage-thresholding algorithm (FISTA) [51], which will be exploited in the future.

4. RESULTS AND DISCUSSION

Experimental validation. The single-shot method A was validated in the visible domain on a broadband and spatially incoherent source with its spectrum full-width at half-maximum ($\mathrm{FWHM}=280\mathrm{nm}$). Then, method B was tested under the same experimental conditions with additional PSF information. The objects employed for these experiments are ubiquitous in the field of scattering imaging and allow for algorithmic validation and resolution estimation.

A. Experimental Validation without PSF Calibration

A lens imaged lithographic letter “B” is shown in Fig. 4(a). The experiment was performed using a commercial broadband light emitting diode (LED) (Thorlabs, MBB1L3, wavelength range of 400–1100 nm, $\mathrm{FWHM}=280\mathrm{nm}$). The collimated light is incident on the object and then transmitted 24 cm to the diffuser (Edmund, 45-656). The light scattered from the diffuser was recorded by a scientific complementary metal oxide semiconductor (sCMOS) camera (Pco.edge. 4.2, pixel size 6.5 μm) set at 12 cm from the diffuser. Further details on the experimental setup can be found in Appendix C.

 

Fig. 4. Experimental validation without calibrating PSF with an illumination bandwidth of $\mathrm{\Delta }\lambda /\lambda =44.8\%$. (a) Photography of a lithographic letter “B” via lens imaging. Reconstructed results by our method (b), Katz et al. (c), and Wu et al. (d). (e)–(h) Ground truth and the retrieved Fourier phase via the three methods. GT, ground truth. Scale bars are 40 camera pixels in all subgraphs.

Download Full Size | PDF

Our single-shot method was then applied to the broadband pattern, at the center of mass of the broadband spectrum ${\lambda }_c=625\mathrm{nm}$ ($\mathrm{\Delta }\lambda /\lambda =44.8\%$). Two other state-of-the-art single-shot methods were also introduced for comparison. The SCI method pioneered by Katz et al. performs phase retrieval on speckle autocorrelation. For each loop, the HIO algorithm [41] was run with a decreasing $\beta$ factor from $\beta =2$ to $\beta =0$, in steps of 0.04. For each $\beta$ value, 30 iterations of the HIO algorithm were performed. The obtained result was fed as an input to additional 30 iterations of the “error reduction (ER)” algorithm to obtain the final result. The second reconstruction method in comparison from bispectrum analysis [35] utilizes the phase closure theory to extract the phase information from subblock speckles, in which the phase is projected onto 18 uniformly spaced angles from (0, $2\pi$). Moreover, the time costs of the reconstruction of speckle sizes of $1800\times 1800=3.24$ million pixels are listed at the bottom of Fig. 4. The experiments were all performed with MATLAB 2018b in Windows 10 running on a dual-core chip Intel Core i9-10900K and 128G memory.

Apparently, both methods [Figs. 4(c) and 4(d)] failed to reconstruct the object and the corresponding Fouier phase under broadband illumination. From the correlation perspective, the failure of the speckle correlation is attributable to the fact that the autocorrelation of broadband PSF cannot be approximated as a delta function, as described in Eq. (6) and mentioned in Fig. 1. Therefore, the estimation of object autocorrelation is stained and further obstructs phase retrieval. The statistical noise of the bispectrum is higher than the autocorrelation, so more speckle spots are required in reconstruction, making it difficult to cope with low-contrast broadband patterns. The proposed method successfully reconstructs the target information, where the PSNR value is as high as 22.63 dB. The similarity between the object Fourier phase and the reconstructed result [see Figs. 4(e) and 4(f)] is obvious.

It is worth emphasizing that iterative phase retrieval algorithms are highly sensitive to the initial guess, the support constraint, and the convergence criterion. The twin-image problem will interfere with the direction selection in the iterative process, causing the phase retrieval algorithm to wander left and right and unable to quickly jump out of the local minimum. Thus the tuning rules are very general and require experiments during a trial and error method.

Here, thanks to the initial non-centrosymmetric support region estimation, adaptive shrink-wrap strategy, and hybrid retrieval method, we achieve a reconstruction accuracy as high as 85%. The empirical probability of success is an average over 1000 trials, in each of which we generate new random matrices as the algorithm input. We declare that a trial is successful if the SSIM between the ground truth and the reconstruction is $>0.8$ (see Visualization 1 for a dynamic reconstruction process). Therefore, verifying our approach is highly repeatable and stable. Our method can effectively avoid brute force manual post-selection and considerably increase the reconstruction fidelity. In terms of algorithm time consumption, although our method has the longest reconstruction time, the performance of our results is the best. In algorithm processing, 90% of the time is occupied on NLM filtering, $>9\%$ of the time is spent on LRSD, and $<1\%$ is spent on phase retrieval. Subsequent research can integrate the algorithm in the graphics processing unit (GPU) or on the field programmable gate array (FPGA) to increase the processing speed, and it is anticipated to achieve high-speed processing.

B. Experimental Validation with PSF Calibration

The main reconstruction comparison between different state-of-the-art approaches is presented in Fig. 5. A tailored USAF1951 target (Edmund, 55-622) with group 0 element 3-5 is selected as the ground truth object. The broadband PSF is measured with a 100 μm pinhole in the same spatial position of the object. The reconstructed results by different algorithms, including the cross-correlation, nonlinear reconstruction, Edrei et al. [34], and our computational reconstruction method, are shown in Figs. 5(b)–5(e), respectively.

 

Fig. 5. Experimental validation with PSF calibration under broadband illumination ($\mathrm{\Delta }\lambda /\lambda =44.8\%$). (a) Photography of a tailored USAF target via lens imaging. Reconstructed results by cross-correlation (b), nonlinear reconstruction with $\alpha =0.4$, $\beta =0.3$ (c), Edrei et al. [34] (d), and our method (e). (f) Line profiles along the yellow arrow in (a) and the same spatial position in each subgraph (b)–(e). GT, ground truth; CC, cross-correlation; NL, nonlinear reconstruction. Scale bars are 40 camera pixels in all subgraphs.

Download Full Size | PDF

As indicated, all methods find the object pattern correctly. Our reconstruction matches perfectly with the ground truth, with PSNR further increasing to 18.57 dB, higher than that of the nonlinear algorithm (15.36 dB). We compare in Fig. 5(f) the intensity profiles along the line indicated by the yellow arrow [see Fig. 5(a)] in each subgraph. Benefiting from the constraint of the ${\ell }_1$ regular term, our results show a significant difference between the target and background. The use of TV regularization effectively reduces the background level and compensates for the artificial traces caused by power spectrum modulation [see Figs. 3(f) and 5(c)]. It is shown that our reconstruction has the smallest background while keeping the highest image contrast. The hollow numbers “4” and “5” close to the edge of the resolution target are also preserved.

An alternative reconstruction is from the Lucy–Richardson deconvolution algorithm under 400 iterations [52] [Fig. 5(d)]. Due to the inherent linear convolution model, the deconvolution algorithm can produce an acceptable result [see Fig. 5(d)]. Correlagraphy-based reconstructions [Figs. 5(b) and 5(c)] also show clear and sharp-edged results with high reconstruction speed. The effects of background noise and artifacts are well compensated by our improved method [Figs. 5(e) and 5(f)]. This comparison also reveals that the correlation-based methods {with Fig. 5(e) or without regularization [Figs. 5(b) and 5(c)]} do not require precise optical calibration while performing well in broadband scattering imaging [53].

In addition, the memory complexity of our method is $\mathcal{O}(N^2)$, and the time complexity of each iteration is $\mathcal{O}(N^2\log N)$. Although the cost time is slightly longer when processing on the CPU (see Fig. 5), it can be dramatically shortened by using parallel computation on GPU (approximately 10 s for four typical iterations on GeForce RTX 2080 Ti). Further comparisons of different approaches under narrowband and broadband illumination are arranged in Appendices B to E.

Due to the temporal and spectral reciprocity, ultrafast sources can also produce equivalent broadband spectra. Naturally, our method can be used as a plug-and-play image restoration in ultrafast science and passive sensing. To demonstrate this point, we show here a comparison among three types of used scientific broad-spectrum light sources and our methods feasibility from numerical simulation. Figure 6(a) illustrates a transform-limited femtosecond pulse $\mathrm{\Delta }t\mathrm{\Delta }\omega \sim 1$ with a 100 fs pulse width. The corresponding spectrum ranges from 480 nm to 720 nm, and $\mathrm{FWHM}=100\mathrm{nm}$ ($\mathrm{\Delta }\lambda /\lambda =16.7\%$). The second revolutionary broadband source is an optical frequency comb with a fixed comb $\mathrm{line}\mathrm{o}\mathrm{f}20\mathrm{nm}$ and a spectrum ranging from 360 nm to 830 nm. The normalized spectrum of the optical comb is shown in Fig. 6(b) with each comb teeth spectrum $\mathrm{FWHM}=1\mathrm{nm}$. The third broadband source is the “D65” illumination [see the spectrum in Fig. 6(c)], which is commonly used to simulate artificial sunlight ($\mathrm{\Delta }\lambda /\lambda =39.7\%$). The broad-spectrum PSFs and speckles generated by these three light sources, and the corresponding imaging results via methods A and B are shown in Figs. 6(d)–6(f), respectively.

 

Fig. 6. Comparison of the reconstructions under three different types of broadband illumination calculated from numerical simulation. (a) The normalized spectrum of a femtosecond pulse with $\mathrm{FWHM}=100\mathrm{fs}$ pulse width ($\mathrm{FWHM}=100\mathrm{nm}$ spectrum). (b) An optical comb with a fixed comb line of $20\mathrm{nm}$ and the spectrum spanning from 360 nm to 830 nm. (c) The “D65” light source spectrum used to simulate artificial sunlight. (d)–(f) Broadband PSF, speckle, and the recovered objects from both methods (methods A and B) corresponding to each light source in (a)–(c). The broadband PSF autocorrelation patterns are inserted in the bottom left of the first column of (d)–(f). The scale bars in (d)–(f) are 500 camera pixels in the left two columns and 30 camera pixels in the rightmost column.

Download Full Size | PDF

All three broadband speckles produce high-fidelity results with both methods. Intriguingly, our noninvasive method (method A) performs similarly to the invasive method (method B). A possible explanation for this behavior could be the noise-suppressing properties of the NLM method used in cepstrum filtering. Further discussion of the speckle contrast and the imaging performance changes with spectral spread is arranged in Appendix B.

While both methods A and B provide promising solutions to the lensless broadband scattering imaging problem, the hidden scenes are restricted to two-dimensional planar objects and the effective field-of-view still lies in the OME scope. For method A, recovering the high-frequency features of the spectrum remains a great challenge as the complete scattering degradation process is unknown. As a result, this noninvasive method cannot be used to image dense and connected samples. This is also a limitation that applies to all existing statistical-based scattering imaging methods, including method A. For method B, there is a trade-off between the detail of the reconstruction and the degree of denoising owing to the introduction of the TV regularity term. Moreover, the complexity of the algorithm needs to be further reduced. We believe that these drawbacks will be overcome in the future.

With respect to improving imaging performance and capability on challenging scenarios, the next step is to implement our method on a complex 3D scene by leveraging the stereo vision [54] or light field solutions [55]. Higher-fidelity reconstruction is anticipated by applying deep learning modules or scattering matrix inversion to unveil the detailed scattering process. Benefiting from the single-shot data acquisition strategy, future works can also be expanded to the dynamic scattering region and non-line-of-sight imaging.

5. CONCLUSION

In summary, we have proposed two individual lensless broadband scattering imaging techniques, which both produce high-fidelity reconstructions from a single-shot wide-spectrum illuminated speckle ($\mathrm{FWHM}=280\mathrm{nm}$) spanning the visible region. The noninvasive method (method A) reproduces the object Fourier spectrum by numerical cepstrum refinement combined with a modified phase retrieval strategy. The invasive method (method B) is motivated by the correlaion, which outperforms state-of-the-art approaches benefiting from the optimization framework and suitable regularizers. No additional measurement of the source spectrum is needed; the spectrum may be ultra-broadband, highly structured, and even exhibit significant intensity fluctuations in both methods. It is expected that the possibility of using broadband illumination in lensless scattering imaging demonstrated here will open up a wide range of opportunities in ultrafast science and astronomical observation and even be extended in the X-ray imaging regime [56] and non-line-of-sight imaging [57]. Additionally, it can also be implemented in unknown target tracking by incorporating cross-correlation analysis [58] or applied in other dynamic scattering scenarios for its single-shot capability [59].

APPENDIX A: TWO-WAVELENGTH MUTUAL INTENSITY CORRELATION FUNCTION

In this section, we present additional mathematical details about the intensity cross-correlation part in Eq. (6) (main text). To begin, we consider the transmission geometry shown in Fig. 7(a). Assuming paraxial propagation condition, the complex amplitude $E(\mathbf{r},z)$ of the light at position $(x,y)$ in the observation plane is related to the complex amplitude $E({\mathbf{r}}^{\prime },0)$ of the scattered light at coordinates $(\xi ,\eta )$ in a plane just to the right of the scattering surface by the Fresnel diffraction integral.

 

Fig. 7. Transmission scattering geometry and coordinate. (a) Free space transmission scattering geometry. (b) Incidence and observation directions, (c) $k$-vectors of the incident light $k_i$, the transmitted light $k_t$, and a $k$-vector in the observation direction $k_o$.

Download Full Size | PDF

Because we only consider the modulus of the correlation function, the quadratic-phase complex exponential is ignored outside the integral as Eq. (A1):

(A1)$$E_k(\mathbf{r},z)=\frac{1}{j\lambda z}{\int }_{{\mathbf{r}}^{\prime }}{E}_k({\mathbf{r}}^{\prime },0)e^{\frac{jk}{2z}{(\mathbf{r}-{\mathbf{r}}^{\prime })}^2}{\mathrm{d}}^2{\mathbf{r}}^{\prime }.$$

The field in the $({\mathbf{r}}^{\prime },0)$ plane can be written as

(A2)$$E_k({\mathbf{r}}^{\prime },0)=tS({\mathbf{r}}^{\prime },0)\exp [j\varphi ({\mathbf{r}}^{\prime },0)]\exp [-jk(\widehat{t}·\widehat{{\mathbf{r}}^{\prime }})r^{\prime }],$$

where $t$ is the average amplitude transmission of the scattering media, and $S$ represents the shape and amplitude distribution across the incident spot. The top mark “^” means the direction of the vector. The phase shift $\varphi ({\mathbf{r}}^{\prime },0)$ imparted to the transmitted wave by the rough surface [see Fig. 7(c)], as measured in the $({\mathbf{r}}^{\prime },0)$ plane, is

(A3)$$\varphi ({\mathbf{r}}^{\prime },0)=k(-n\widehat{t}·\widehat{z}+\widehat{o}·\widehat{z})h({\mathbf{r}}^{\prime },0),$$

where here and throughout, $\lambda$ is the free-space wavelength and $n$ is the refractive index of the scattering media. We define here a scattering vector $\overrightarrow{q}$, decomposited into transverse ${\overrightarrow{q}}_t$ and vertical ${\overrightarrow{q}}_z$ components:

(A4)$$\overrightarrow{q}={\overrightarrow{k}}_o-{\overrightarrow{k}}_t={\overrightarrow{q}}_t+q_z\widehat{z}.$$

The phase shift suffered at the rough interface is rewritten as

(A5)$$\varphi ({\mathbf{r}}^{\prime },0)=q_{z}h({\mathbf{r}}^{\prime },0)\mathrm{.}$$

Our goal is to calculate the autocorrelation function ${\mathbf{\Gamma }}_E$ of the speckle field $E(\mathbf{r},z)$ at two points $({\mathbf{r}}_1,z)$ and $({\mathbf{r}}_2,z)$, that is,

(A6)$${\mathbf{\Gamma }}_E({\mathbf{r}}_1;{\mathbf{r}}_2,z)=\overline{E_{k_1}({\mathbf{r}}_1,z)E_{k_2}^{*}({\mathbf{r}}_2,z)}.$$

The cross correlation between two observed fields $E_{k_1}({\mathbf{r}}_1,z)$ and $E_{k_2}({\mathbf{r}}_2,z)$ can now be written as

(A7)$$\begin{gathered} {\mathbf{\Gamma }}_E({\mathbf{r}}_1;{\mathbf{r}}_2,z)=\frac{1}{{\lambda }_1{\lambda }_2{z}^2}{\int }_{{\mathbf{r}}_1^{\prime }}{\int }_{{\mathbf{r}}_2^{\prime }}{\mathbf{\Gamma }}_E({\mathbf{r}}_1^{\prime };{\mathbf{r}}_2^{\prime },0)\times \exp [\frac{j}{2z}(k_1{\mathbf{r}}_1^{\prime 2}-k_2{\mathbf{r}}_2^{\prime 2})] \\ \times \exp [-\frac{j}{z}(k_1{\mathbf{r}}_1{\mathbf{r}}_1^{\prime }-k_2{\mathbf{r}}_2{\mathbf{r}}_2^{\prime })]{\mathrm{d}}^2{\mathbf{r}}_1^{\prime }{\mathrm{d}}^2{\mathbf{r}}_2^{\prime }. \end{gathered}$$

Substituting Eq. (A1) in the first line of this equation, we obtain

(A8)$$\begin{gathered} {\mathbf{\Gamma }}_E({\mathbf{r}}_1^{\prime };{\mathbf{r}}_2^{\prime },0)=\overline{E_{k_1}({\mathbf{r}}_1^{\prime },0)E_{k_2}^{*}({\mathbf{r}}_2^{\prime },0)} \\ ={|t|}^{2}S({\mathbf{r}}_1^{\prime },0)S^{*}({\mathbf{r}}_2^{\prime },0)\exp \{\overline{j[q_{z1}h({\mathbf{r}}_1^{\prime })-q_{z2}h({\mathbf{r}}_2^{\prime })}]\} \\ ={|t|}^{2}S({\mathbf{r}}_1^{\prime },0)S^{*}({\mathbf{r}}_2^{\prime },0)M_h(q_{z1},-q_{z2}) \\ \times \exp [-j(k_1{\widehat{t}}_1·\widehat{{\mathbf{r}}_1^{\prime }})r_1^{\prime }+j(k_2{\widehat{t}}_2·\widehat{{\mathbf{r}}_2^{\prime }})r_2^{\prime }], \end{gathered}$$

where $M_h(A,B)$ is the second-order characteristic function of the surface-height function $h({\mathbf{r}}^{\prime })$, which is assumed to be statistically stationary. We assume that the correlation width of the scattered wave is sufficiently small that we can approximate ${\mathbf{\Gamma }}_E$ as follows:

(A9)$${\mathbf{\Gamma }}_E({\mathbf{r}}_1^{\prime };{\mathbf{r}}_2^{\prime },0)\approx \kappa {|\mathbf{t}|}^2{|\mathbf{S}({\mathbf{r}}_1^{\prime },0)|}^2{\mathbf{M}}_h(\mathrm{\Delta }{q}_z)\times \exp [-j(k_1{\widehat{t}}_1·\widehat{{\mathbf{r}}_1^{\prime }})r_1^{\prime }+j(k_2{\widehat{t}}_2·\widehat{{\mathbf{r}}_2^{\prime }})r_2^{\prime }]\delta ({\mathbf{r}}_1^{\prime }-{\mathbf{r}}_2^{\prime }).$$

Substituting Eq. (A9) in Eq. (A7) and using the small angles approximation, we obtain

(A10)$$\begin{gathered} {\mathbf{\Gamma }}_E({\mathbf{r}}_1;{\mathbf{r}}_2,z)=\frac{\kappa {|t|}^2}{{\lambda }_1{\lambda }_2{z}^2}{\mathbf{M}}_h(\mathrm{\Delta }{q}_z){\int }_{{\mathbf{r}}^{\prime }}|\mathbf{S}({\mathbf{r}}^{\prime }{,0)|}^2\exp \left[\frac{j(k_1-k_2){\mathbf{r}}^{\prime 2}}{2z}\right] \\ \times \exp [-j(k_1{\widehat{t}}_1·\widehat{{\mathbf{r}}^{\prime }})r^{\prime }+j(k_2{\widehat{t}}_2·\widehat{{\mathbf{r}}^{\prime }})r^{\prime }] \\ \times \exp [-j(k_1{\widehat{o}}_1·\widehat{{\mathbf{r}}^{\prime }})r^{\prime }+j(k_2{\widehat{o}}_2·\widehat{{\mathbf{r}}^{\prime }})r^{\prime }]{\mathrm{d}}^2{\mathbf{r}}^{\prime }. \end{gathered}$$

To simplify Eq. (A10), we introduce Eq. (A11) notationally:

(A11)$${\overrightarrow{q}}_{t1}=k_1{\widehat{t}}_1·\widehat{{\mathbf{r}}^{\prime }}+k_1{\widehat{o}}_1·\widehat{{\mathbf{r}}^{\prime }}{\overrightarrow{q}}_{t2}=k_2{\widehat{t}}_2·\widehat{{\mathbf{r}}^{\prime }}+k_2{\widehat{o}}_2·\widehat{{\mathbf{r}}^{\prime }}.$$

With the help of Eq. (A11), this result can be simplified as

(A12)$${\mathbf{\Gamma }}_E({\mathbf{r}}_1;{\mathbf{r}}_2,z)=\frac{\kappa {|t|}^2}{{\lambda }_1{\lambda }_2{z}^2}{\mathbf{M}}_h(\mathrm{\Delta }{q}_z)\mathrm{\Psi }(\mathrm{\Delta }{\overrightarrow{q}}_t){\int }_{{\mathbf{r}}^{\prime }}|\mathbf{S}({\mathbf{r}}^{\prime }{,0)|}^2{\mathrm{d}}^2{\mathbf{r}}^{\prime },$$

where $\mathrm{\Delta }{\overrightarrow{q}}_t={\overrightarrow{q}}_{t1}-{\overrightarrow{q}}_{t2}$ and

(A13)$$\mathrm{\Psi }(\mathrm{\Delta }{\overrightarrow{q}}_t)=\frac{{\int }_{{\mathbf{r}}^{\prime }}|\mathbf{S}({\mathbf{r}}^{\prime },0){|}^2\exp [\frac{j(k_1-k_2){\mathbf{r}}^{\prime 2}}{2z}]\exp [-j(\mathrm{\Delta }{\overrightarrow{q}}_t·{\overrightarrow{r}}_t){\mathrm{d}}^2{\mathbf{r}}^{\prime }}{{\int }_{{\mathbf{r}}^{\prime }}|\mathbf{S}({\mathbf{r}}^{\prime },0){|}^2{\mathrm{d}}^2{\mathbf{r}}^{\prime }},$$

is the normalized Fourier transform of the intensity distribution multiplied by a dual-wavelength impacted quadratic phase factor:

(A14)$$\mathrm{\Psi }(\mathrm{\Delta }{\overrightarrow{q}}_t)=\frac{\mathcal{F}\{|\mathbf{S}({\mathbf{r}}^{\prime }{,0)|}^2\exp [\frac{j(k_1-k_2){\mathbf{r}}^{\prime 2}}{2z}\left]\right\}}{{\int }_{{\mathbf{r}}^{\prime }}|\mathbf{S}({\mathbf{r}}^{\prime }{,0)|}^2{\mathrm{d}}^2{\mathbf{r}}^{\prime }}.$$

Here, ${|\mathbf{S}|}^2$ is the intensity distribution across the scattering spot, which is given by $|\mathbf{S}({\mathbf{r}}^{\prime }{,0)|}^2=I({\mathbf{r}}^{\prime },0)$. To make further progress, a specific shape for the scattering spot must be assumed. Let the spot be uniform and circular with diameter $D$, centered on the origin. Using the Hankel transform, we obtain

(A15)$$\mathrm{\Psi }(\mathrm{\Delta }{\overrightarrow{q}}_t)=2\frac{J_1\left(\frac{D|\mathrm{\Delta }{\overrightarrow{q}}_t|}{2}\right)}{\frac{D|\mathrm{\Delta }{\overrightarrow{q}}_t|}{2}}*\mathcal{F}\left\{\exp \right[\frac{j(k_1-k_2){\mathbf{r}}^{\prime 2}}{2z}\left]\right\}.$$

For further simplification, let us introduce the same approximation used in the calculation of monochromatic speckle:

(A16)$$\frac{J_1\left(\frac{D|\mathrm{\Delta }{\overrightarrow{q}}_t|}{2}\right)}{\frac{D|\mathrm{\Delta }{\overrightarrow{q}}_t|}{2}}\approx \frac{J_1(\pi Dr/{\overline{\lambda }}_{ij}z)}{\pi Dr/{\overline{\lambda }}_{ij}z}.$$

Here ${\overline{\lambda }}_{ij}=({\lambda }_i+{\lambda }_j)/2$ means the average value of the two wavelengths and the subscripts $i$ and $j$ represent the $i$-th and $j$-th wavelengths. The normalized cross-correlation function, ${\mathbf{\mu }}_{\mathit{E}}$ is found by dividing ${\mathbf{\Gamma }}_{\mathit{E}}$ by its value when $\mathrm{\Delta }{q}_z=0$ and $\mathrm{\Delta }{q}_t=0$. The result is

(A17)$${|{\mathbf{\mu }}_E({\overrightarrow{q}}_1,{\overrightarrow{q}}_2)|}^2={|{\mathbf{M}}_h(\mathrm{\Delta }{q}_z)|}^2{|\mathrm{\Psi }(\mathrm{\Delta }{\overrightarrow{q}}_t)|}^2.$$

Obviously, the wavelength variation has two effects. First the observed speckle pattern undergoes a spatial scaling by an amount that depends on distance from the point where the scattering event occurs ($\mathrm{\Psi }$). When the wavelength increases, the speckle pattern expands about this point, and when the wavelength decreases, the speckle pattern shrinks about this point. This effect is caused by the fact that diffraction angles are inversely proportional to wavelengths. A second effect is upon the random phase shifts imparted to the scattered wave by the surface-height fluctuations (${\mathbf{M}}_h$). Since those phase shifts are proportional to $h/\lambda$, an increase in wavelength will decrease the random phase shifts, and vice versa. When the surface is rough on the scale of a wavelength, phase wrapping into the interval (0, $2\pi$) will result in changes in the detailed structure of the speckle pattern (see Visualization 1 for details).

If the surface-height fluctuations obey Gaussian statistics, the phase shifts imparted intensity correlation change is given by

(A18)$${|{\mathbf{M}}_h(\mathrm{\Delta }{q}_z)|}^2=\exp (-{\sigma }_h^2\mathrm{\Delta }{q}_z^2).$$

In the case that the incidence direction and observation direction are both normal,

(A19)$$\mathrm{\Delta }{q}_z=(n-1)|k_i-k_j|.$$

According to the Siegert relation [60], the intensity correlation function can be associated with the electric field correlation function as

(A20)$${\mathbf{\Gamma }}_I(\mathrm{\Delta }\mathbf{r})={\overline{I}}^2[1+{|{\mathbf{\mu }}_E(\mathrm{\Delta }\mathbf{r})|}^2].$$

We can now write the intensity cross-correlation function as

(A21)$$\begin{gathered} {\sum }_{i=1}^M{\sum }_{j\ne i}^M{\omega }_i{\omega }_j{\mathrm{PSF}}_N({\lambda }_i,\mathrm{\Delta }{\lambda }_i)\star {\mathrm{PSF}}_N({\lambda }_j,\mathrm{\Delta }{\lambda }_j)\propto {\sum }_{i=1}^M{\sum }_{j\ne i}^M{\omega }_i{\omega }_j\{1+\exp [-{\sigma }_h^2(n-1)|k_i-k_j|] \\ \times {\left\{2\frac{J_1(\pi Dr/{\overline{\lambda }}_{ij}z)}{\pi Dr/{\overline{\lambda }}_{ij}z}*\mathcal{F}\left\{\exp \left[\frac{j(k_i-k_j){\mathbf{r}}^{\prime 2}}{2z}\right]\right\}\right\}}^2\}. \end{gathered}$$

This equation represents the final result of our analysis of the cross-correlation between two arbitrary wavelength speckles. When $i$ equals $j$, the Eq. (A21) degenerates to Eq. (5) (main text), which indicates the correlation between two speckles at the same wavelength.

APPENDIX B: IMAGING PERFORMANCE WITH INCREASING ILLUMINATION BANDWIDTH

1. Simulation Process

To support the experimental data, we have performed full numerical simulations of wave propagation in disordered media. In the simulations, the diffuser is inserted at a distance $v=12\mathrm{cm}$ in front of the detector and $u=24\mathrm{cm}$ behind the object. The geometry of propagation from the source plane ${\mathbf{E}}_k({\mathbf{r}}_0,0)$ to the observation plane ${\mathbf{E}}_k({\mathbf{r}}_2,z)$ can be divided into the following three stages: i) the source plane ${\mathbf{E}}_k({\mathbf{r}}_0,0)$ propagates to the scattering plane ${\mathbf{E}}_k({\mathbf{r}}_1,u)$ via free-space propagation; ii) the wavefront is modulated by the diffuser; iii) the disordered wavefront propagates to the detector plane ${\mathbf{E}}_k({\mathbf{r}}_2,z)$ via free-space propagation. The $k$ and ${\mathbf{r}}_0$, ${\mathbf{r}}_1$, ${\mathbf{r}}_2$ stand for the wave vector and the three plane coordinates, respectively. The free-space propagation can be modeled as Fresnel diffraction integral, where the notation here is adapted from that described by Nazarathy and Shamir [61]:

(B1)$$\mathcal{R}[z,{\mathbf{r}}_0,{\mathbf{r}}_1]\{{\mathbf{E}}_k({\mathbf{r}}_0)\}=\frac{1}{j\lambda z}{\int }_{{\mathbf{r}}_0}{\mathbf{E}}_k({\mathbf{r}}_0,0)e^{\frac{jk}{2z}{({\mathbf{r}}_1-{\mathbf{r}}_0)}^2}{\mathrm{d}}^2{\mathbf{r}}_0.$$

The operators’ parameters are given in square brackets, and the operand is given in curly braces. We model the diffuser transmittance $t_k({\mathbf{r}}_1)$ as a pure-phase random mask, i.e., $t_k({\mathbf{r}}_1)=\exp [ik\mathrm{\Delta }nh({\mathbf{r}}_1)]$, where $h({\mathbf{r}}_1)$ is the random height of the diffuser surface and $\mathrm{\Delta }n$ is the difference between the refractive indices of the diffuser and the surrounding air ($\mathrm{\Delta }n\approx 0.52$ for glass diffusers). A circle-shaped pupil $p({\mathbf{r}}_1)$ is added behind the diffuser to control the speckle grain size (the diameter of the pupil $D=6\mathrm{mm}$ is used in simulation).

The Fresnel diffraction operator $\mathcal{R}$ propagates only monochromatic waves. Any broadband signal can be propagated by first writing it as a superposition of monochromatic waves and then propagating each one individually. According to the above description, the broadband PSF can be written as

(B2)$${\mathrm{PSF}}_B={\left|{\int }_k\mathcal{R}[v,{\mathbf{r}}_1,{\mathbf{r}}_2]\{\mathcal{R}[u,{\mathbf{r}}_0,{\mathbf{r}}_1]\{E_k({\mathbf{r}}_0,0)\}{t}_k({\mathbf{r}}_1)p({\mathbf{r}}_1)\}\mathrm{d}k\right|}^2.$$

The emulation results of speckles with various illumination bandwidths are shown in Fig. 8. Figure 8(a) shows the Gaussian-shaped spectral profiles with $\mathrm{FWHM}=1$, 7, 24, 46, and 104 nm and the corresponding standard deviation (SD) of 0.4, 3, 10, 20, and 45 nm, respectively, all centered at 632.8 nm. $\mathrm{\Delta }\lambda =0\mathrm{nm}$ indicates monochromatic light at 632.8 nm. The corresponding autocorrelation intensity profiles of these broadband PSFs are presented in Fig. 8(b). Figure 8(c) is a neuron-like object that serves as the hidden object. Notably, with the increase in illumination bandwidth, the PBR of speckle autocorrelation declines quickly at first and then becomes slow-changing and stable. When the illumination spectrum width becomes infinitely wider, the PBR will be infinitely close to 1, as can be expected. Speckles under different illumination bandwidths are presented in Fig. 8(d). The speckle contrast ratios (SCRs) [26] of these speckles are calculated as 0.99, 0.79, 0.38, 0.23, 0.16, and 0.12 for illumination from $\mathrm{\Delta }\lambda =0\mathrm{nm}$ to $\mathrm{\Delta }\lambda =104\mathrm{nm}$, respectively.

 

Fig. 8. Emulation results under different illumination spectral widths. (a) Normalized spectrum profile of illumination bandwidth FWHM varying from $\mathrm{\Delta }\lambda =0\mathrm{nm}$ to $\mathrm{\Delta }\lambda =104\mathrm{nm}$. (b) The autocorrelation intensity profiles of the corresponding broadband PSF patterns with the illumination shown in (a). (c) Neuron cell-like pattern used to generate speckles. (d) Speckle patterns under different illumination in (a). The scale bar is 20 camera pixels in (c) and 300 camera pixels in (d).

Download Full Size | PDF

2. Imaging Performance versus Illumination Bandwidth, Method A

Reconstructions from speckles in Fig. 8(d) are shown in Fig. 9(c) via method A. Note that the reconstruction performance decreases slightly as a result of the participation of a wider spectral bandwidth. The loss of the high-frequency information in the broadband illuminated pattern [Fig. 8(d)] causes the cutoff frequency of the reconstructed Fourier spectrum to gradually shrink to the center [Fig. 9(e)], and correspondingly the estimated support region becomes increasingly loose [Fig. 9(d)].

 

Fig. 9. Scattering imaging results under different illumination bandwidths via our proposed method A using simulated data. (a) Ground truth object. (b) Fourier amplitude of (a). (c) Recovered results using speckles under different illumination bandwidths [Fig. 8(d)]. (d) Estimated support regions of (c). (e) Retrieved Fourier amplitude corresponding to (c). (f) Quantitative performance of (c) with (a) as a reference. GT, ground truth. Scale bars are 20 camera pixels in (a), (c), and (d) and 200 camera pixels in (b) and (e).

Download Full Size | PDF

3. Imaging Performance versus Illumination Bandwidth, Method B

When introducing broadband PSF information, the ill-posed problem can be effectively mitigated. Images recovered via method B under different illumination conditions are presented in Fig. 10, showing very few differences in resolution. However, our method performs similarly to the narrowband case and even better in the broadband regime. A possible explanation for this behavior could be the smoothing properties of the joint regularization method, which slightly erases details in the reconstruction.

 

Fig. 10. Imaging results with quantitative performance under different illumination bandwidths via our proposed method B using simulated data. (a) Ground truth object. (b)–(g) Recovered hidden objects using speckles from Fig. 8(d). GT, ground truth. Scale bars are 20 camera pixels in all subgraphs.

Download Full Size | PDF

As mentioned here and in the paper, one can in principle use a detection bandwidth much larger than the speckle spectral decorrelation bandwidth, and still be able to effectively employ our proposed approaches.

APPENDIX C: SYSTEM CONFIGURATION AND EXPERIMENT RESULTS

1. Experimental Setup

The experimental setup for imaging in transmission is shown in Fig. 11. The broadband light is provided by a white light LED (Thorlabs, MBB1L3, wavelength range of 400–1100 nm, $\mathrm{FWHM}=280\mathrm{nm}$), which has a central wavelength of 625 nm and a coherence length of $L_c\approx 1.4\mathrm{\mu m}$. The light was collimated and transmitted through an object and a scattering slab and imaged onto the camera (Pco.edge 4.2, 6.5 μm pixel size). We used a USAF resolution target (Edmund, USAF 1951) and fabricated custom patterns via standard photolithography masks as objects. Two types of scattering media [Edmund, ground glass 220-grit (Edmund, 45-656), 600-grit (Edmund, 38-790)] were rotationally placed at 240 mm behind the object and 120 mm in front of the camera. A typical integration time for a single-shot capture was 500 ms. The PSF of the scattering system was measured by replacing the object with a pinhole of a typical diameter of 100 μm. The narrow-spectrum light sources used in the experiment are red light LED (Thorlabs, M625L4, ${\lambda }_c=625\mathrm{nm}$, $\mathrm{FWHM}=17\mathrm{nm}$) and green light LED (Thorlabs, M530L4, ${\lambda }_c=530\mathrm{nm}$, $\mathrm{FWHM}=35\mathrm{nm}$).

 

Fig. 11. Experimental setup. A scattering slab is illuminated by a wide-spectrum LED light source with a spectral range from 400 to 1100 nm ($\mathrm{FWHM}=280\mathrm{nm}$). The broadband scrambled patterns are recorded lenslessly by a scientific CMOS camera.

Download Full Size | PDF

2. Broadband Imaging through Scatterers via Cross Correlation

The first set of experiments were conducted to test the capability of the cross-correlation technique. We show the ground truth objects in the leftmost column of Fig. 12. The 2nd–4th columns present the recovered objects from the white, green, and red light sources using a 220-grit glass diffuser, respectively. Same procedures are performed on a 600-grit diffuser, and the corresponding results are depicted in the right three columns. All the PSFs are captured using a 100 μm pinhole in this experiment. It can be observed that all the reconstruction results lie in the opposite direction to the ground truth object. This is due to the calculation process of the cross-correlation operation. Further investigation about this point has been fully studied in Ref. [62]. Under the 220-grit ground glass, the reconstruction results of broad-spectrum light show a similar performance to that of narrow-spectrum light, but for the 600-grit ground glass, which has a lower degree of scattering, the reconstruction performance under broad-spectrum source is slightly worse than that of narrow-spectrum sources. This is due to the lower PBR as described in Fig. 8.

 

Fig. 12. Experimental comparison of broadband and narrowband illuminated imaging performance via cross-correlation strategy. First column, photography of ground truth objects. Columns 2–5 show imaging results from 220-grit ground glass diffuser with broadband illumination ($\mathrm{\Delta }\lambda =280\mathrm{nm}$) and narrowband irradiance centered at ${\lambda }_c=530\mathrm{nm}$ ($\mathrm{\Delta }\lambda =35\mathrm{nm}$) and ${\lambda }_c=625\mathrm{nm}$ ($\mathrm{\Delta }\lambda =17\mathrm{nm}$), respectively. Columns 6–8 represent imaging results of 600-grit glass diffuser with the same white, green, and red LED illumination. Scale bars are 50 camera pixels in the rightmost column. Performance comparison of the reconstructed object along the lines indicated by the white arrows is shown next to the results.

Download Full Size | PDF

3. Broadband Imaging through Scatterers via Nonlinear Reconstruction

While the experimental results in the previous section have demonstrated that the cross-correlation approach works well for object reconstruction, directly performing the cross-correlation operation in reconstruction will result in substantial background noise. To suppress the background noise while preserving the effective object information, it is better to perform an advanced nonlinear filter in object retrieval. Smarter improvement in algorithm is about the tunable parameters used for spectrum redistribution. Figure 13(a) illustrates the nonlinear reconstruction results under different values of $\alpha$ and $\beta$ varying between $-1$ and 1 in steps of 0.2. In the nonlinear reconstruction, the magnitudes of the spectrum of the two correlation functions are modified with respect to one another to create an effect equivalent to that of correlating two bipolar intensity distributions.

 

Fig. 13. Nonlinear reconstruction results of the letter “$\mathrm{B}$” object. The left part corresponds to the retrieved object along different tunable parameters. The right part presents a metric matrix to search for the sharpest image from the left part by using a gradient magnitude based evaluation index. The purple, blue, and green boxes show the reconstruction results with matched, phase-only, and inverse filters, respectively. Results above the yellow dotted line show spurious solutions due to the low power of the power spectrum. The gray and pink boxes indicate the best performance produced by the quality assessment function and from naked eyes. Arrows I and II represent the reconstruction changes by the PSF and the speckle modification, respectively. Arrow III represents the reconstruction changes by a fixed sum of $\alpha$ and $\beta$.

Download Full Size | PDF

Owing to the adjustable power spectrum density in the Fourier domain, one can choose proper parameters to balance between the image resolution and the background noise, thus yielding better performance than the cross-correlation method. By searching the optimal $\alpha$ and $\beta$, which for $\beta$ is the range among inverse filter [$\beta =-1$, green box in Fig. 13(b)], the phase-only filter [$\beta =0$, blue box in Fig. 13(b)], and the matched filter [$\beta =1$, purple box in Fig. 13(b)], an optimal reconstruction result can be obtained. The gray box presents the sharpest result predicted by image quality assessment, while the pink box shows the best performance observed with the naked eye. Arrows I and II indicate the variation caused by single parameter changes, and the algorithm will produce gradually improved results when the parameter becomes larger. Arrow III shows a similar reconstruction performance with a fixed sum of $\alpha$ and $\beta$. A better result can be found in the range of $\alpha$ and $\beta$ with a positive value. However, the increase in the weight of $\alpha$ may produce a ghost image [see Fig. 13(a) rightest four columns] surrounding the object, while the increase in the weight of $\beta$ will generate a darker background [see Fig. 13(a) bottom four rows]. We artificially introduced a divided line [see Fig. 13(a) yellow dotted line] to separate the recognizable and unrecognizable results (artificial scenes can be found in the unrecognizable part due to its overreduced $\alpha$ and $\beta$, as indicated in the main text). Figure 13(b) shows a black-like area and a colorful area by utilizing a gradient-based metric adapted from the image quality assessment model [63].

4. Extended Results of Method A

The first row of Fig. 14 shows the retrieved results of broadband patterns under different illumination conditions, and the second row displays two extended recovered results of a star- shaped and number 3 object. A dynamic phase retrieval process can be found in Visualization 1 part IV.

 

Fig. 14. Extended imaging results of method A. (a) Ground truth object “B.” (b)–(d) The recovered letter “B” from red, green, and white light via method A. (e) Ground truth object “star.” (f) Retrieved object “star” under broadband illumination ($\mathrm{\Delta }\lambda /\lambda =44.8\%$). (g) Ground truth object “3.” (f) Retrieved object “3” under broadband illumination ($\mathrm{\Delta }\lambda /\lambda =44.8\%$). GT, ground truth. Scale bars in all sub-graphs are 40 camera pixels.

Download Full Size | PDF

5. Imaging with a Complex Micro-target

To demonstrate the feasibility of this method to image a complex micro-target through scattering layers, we placed a binary target with total size $600\mathrm{\mu m}\times 400\mathrm{\mu m}$ in front of the 220-grit ground glass diffuser [see Fig. 15(e)]. We used a similar setup as shown in Fig. 11, but we decreased the distance between the object and diffuser (ground glass 220-grit) to 15 cm to reduce optical energy transmission loss. The PSF is measured using a 25 μm pinhole. The reconstruction procedure (method B) was applied to the three broadband patterns displayed in Figs. 15(b)–15(d). The corresponding retrieved objects are shown in Figs. 15(f)–15(h).

 

Fig. 15. Complex micro-target imaging via proposed method B. (a) The broadband spectrum (orange) with a bandwidth of 44.8% and the narrowband spectrum (red) with a bandwidth of 2.7% and (green) with a bandwidth of 6.6%. (b)–(d) Speckle pattern under red, green, and white illumination, respectively. (e) Ground truth object. The yellow arrows indicate a micro portion with about 10 μm. (f)–(h) The corresponding retrieved objects from (b)–(d) using the 220-grit ground glass diffuser with a 25 μm pinhole. Scale bars are 300 camera pixels in (b)–(d) and 100 μm in (e)–(h).

Download Full Size | PDF

Using our proposed method, we were able to reconstruct an image with an approximate 10 μm resolution [indicated by the yellow arrows on letter “U” in Fig. 15(e)], which is about half size of the pinhole probe. However, the residual sub-micrometre school badge pattern visible in Fig. 15(e) is not visible as they scale below the spatial resolution. Note that the extended illumination bandwidth will slightly decrease the reconstruction performance, which further verified the applicability of our proposed method.

6. Choice of Parameters for Method A

For method A, the reconstruction quality and empirical probability of success depend on the refined Fourier amplitude. Therefore, choosing appropriate parameters for Fourier space filtering is important, especially the NLM filter. Fortunately, there is a large optimal space for NLM filter parameter selection. We do not observe any remarkable qualitative differences in changing the patch size and research window, which is not shown here. Different reconstruction results can be observed when changing the NLM filtering parameter $\gamma$ (see Fig. 16). For the “Neuron cell” pattern, the $\gamma$ can be chosen from 2 to 9. For more complex targets, the parameter selection space is slightly reduced. After our tests, we found that even though $\gamma$ is slightly larger than the upper bound of the optimal space, the retrieved object is still distinguishable but with a slight edge blur.

 

Fig. 16. Image reconstruction with changing NLM filtering parameters $\gamma$ via method A. The optimal parameter selection range shrinks slightly with increasing target complexity.

Download Full Size | PDF

7. Choice of Parameters for Method B

The parameters ($\tau$ and $\lambda$) described in Eq. (14) will influence the reconstruction results. When $\tau$ is too small, it will lead to noisy results, while a too large $\tau$ will lead to blurry results. However, $\lambda$ that is too small will result in noisy results, whereas $\lambda$ that is too large will result in sparse results. Therefore, it is important to choose suitable parameters for this method.

We have tried a series of values in both experiments and simulations. Part of them are shown in Fig. 17. The results related to $\lambda$ and $\tau$ in Eq. (14) are respectively listed in the corresponding columns and rows. It can be seen from Fig. 17 that the results vary just as described above. With the increase of $\lambda$, the results become sparser, while they turn to be blurry when $\tau$ increases. The parameters used by the method, for different values of $\lambda$ and $\tau$ are displayed in Table 2. When $\lambda =0.001$ and $\tau =0.005$, PSNR and SSIM are 13.56 and 0.51, respectively, which are the best results among the different trials, as shown in the table. Thus, we choose $\lambda =0.001$ and $\tau =0.005$ in this simulation. The parameters of other results are also determined in this way.

 

Fig. 17. Results of different parameters. This figure shows the results reconstructed by method B with different parameters. The leftmost image is the ground truth image used in the simulation. The second to fourth columns represent different $\lambda$, which is used to control the ${\ell }_1$ regularization, while the rows represent different $\tau$, which is used to control the TV regularization term.

Download Full Size | PDF

Table 2. Quantitative Analysis for Different Parameters Using Method B

View Table | View all tables in this article

APPENDIX D: SPATIAL AND SPECTRAL CORRELATION ANALYSIS

1. Measurement of 3D OME Range

To quantify the effective imaging scope of the broadband scattering approach, we measured the translational and longitudinal OME range of the scattering media. The results are summarized in Fig. 18 for two types of ground glass diffusers. A 100 μm diameter pinhole is used as a point source placed at a distance $u=24\mathrm{cm}$ in front of the diffuser and scanned perpendicular to the optical axis direction by a translation stage (MTS50A-Z8, Thorlabs) with a fixed 0.5 mm step. The longitudinal OME range is measured using the same framework but with the pinhole scanned along the optical axis direction with a fixed 1 mm step. The effective OME range is defined as the correlation index dropping from 1 to 0.5, and the results are arranged in Figs. 18(c) and 18(f). It should be noted that all the measurements are selected in the central region of speckles. A detailed OME measurement and calculation process can be found in Ref. [64].

 

Fig. 18. Measurement of the spatial speckle decorrelation as a function of the displacement laterally and longitudinally. (a) and (b) Translational OME range of the 220-grit and 600-grit ground glass, respectively. (d) and (e) Longitudinal OME range of the 220-grit and 600-grit ground glass. (c) and (f) Effective OME range for the two types of ground glass diffusers. All experiments were repeated for three spectral bands: red, green, and white.

Download Full Size | PDF

2. Spectral Decorrelation of Speckle Patterns

To further illustrate the dispersion phenomenon in broadband scattering, we investigate the evolution of the correlation variation when the illumination wavelength changes. The simulation results are presented in Fig. 19. Previous studies indicated that the correlation between two speckles beyond the spectral correlation range should be omitted [38]. However, this is not completely accurate. This is because of the existence of these cross-correlation terms that greatly enhanced the strength of the background term, resulting in the failure of computational techniques such as speckle correlation imaging. The central region of the cross-correlation gradually diffuses and eventually blends with the background [see the inserts in Figs. 19(e)–19(h)], which also reveals the spectral decorrelation process (see Visualization 1 part III for more details).

 

Fig. 19. (a)–(d) Monochromatic speckles at $\lambda =500$, 510, 520, and 530 nm from numerical simulation. (e)–(h) Cross-correlation between different speckle patterns of (a)–(d) with (a) as the reference speckle. Inserts are enlarged portions of the central correlation peak. Notably, the cross-correlation terms spread with increasing illumination wavelength diversity. (i) Spectral decorrelation of the speckle patterns versus illumination spectral width. Scale bars in (a)–(d) are 50 camera pixels.

Download Full Size | PDF

3. Imaging Resolution versus Probe Size

To verify the image resolution with respect to the probe variations, the radius of the pinhole was varied from 50 μm to 200 μm, and the normalized autocorrelation function was plotted [see Fig. 20(e)]. From Figs. 20(a)–20(c), it can be found that the resolving power of the system improves as the pinhole size decreases. This can be derived from the observation of the autocorrelation in Figs. 20(d) and 20(e). Correspondingly, the capable resolution will become weaker with the increasing probe size. Further investigation of the quantitative resolution analysis will be implemented in the future.

 

Fig. 20. Imaging resolution varies with the size of the pinhole probes. (a)–(c) Recovered results from red, green, and white illuminated patterns using 220-grit ground glass with a 50 μm, 100 μm, and 200 μm pinhole, respectively. (d) Broadband PSF autocorrelation as a function of the pinhole size for each PSF corresponding to (a)–(c). (e) Normalized autocorrelation intensity profiles of all PSFs along the arrow presented in (d). Scale bars are 50 camera pixels in (a)–(c).

Download Full Size | PDF

4. Assessment of the Effective Depth-of-field

To verify the effective depth-of-field of the imaging system, we captured a diversity of speckle patterns along the optical axis with a fixed $\delta =2\mathrm{mm}$ spatial interval moving in the direction of the diffuser. The broadband PSF is calibrated at the position $u=24\mathrm{cm}$, and the in focus broadband speckle is measured at the same position. From Fig. 21, we found that the image quality slowly changes across the axial displacement, even with a large defocus distance. Empirical evidence suggests that, in practice, the experimentally measured decay is slower than theoretically predicted, as long as the depth-of-field falls in the longitude OME range. As expected, when the defocus distance increases, more speckle spatial decorrelation is present, and thus the retrieved images become blurred.

 

Fig. 21. Imaging depth-of-field analysis with a fixed interval $\delta =2\mathrm{mm}$ defocus distance. (a)–(c) Imaging results of red, green, and white illuminated at different depths of field under 220-grit ground glass via method B, respectively. The leftmost column indicates the in-focus results with the object and the pinhole at the same spatial position. The right four columns display the defocused results. Scale bars in the rightmost column are 50 camera pixels.

Download Full Size | PDF

APPENDIX E: NOVEL BROADBAND SPECKLE CORRELATION PHENOMENON

1. Anisotropic Broadband Speckle Correlation Characteristic

Figure 22 depicts an anisotropic correlation phenomenon observed in both simulation and experiment. Notably, the entire broadband speckle has an expansion phenomenon from the center to the boundary [see Figs. 22(a) and 22(e)]. The central part of the broadband speckle autocorrelation has a delta-like function as shown in Figs. 22(d) and 22(h). Obviously, the closer to the speckle center, the smaller the average speckle grain size. The speckle grains at the boundary stretched consistently in the spatial direction [see Figs. 22(b) and 22(f)] as well as the autocorrelation function [see Figs. 22(c) and 22(g)].

 

Fig. 22. Anisotropic correlations in broadband speckle patterns. (a) Broadband speckle pattern obtained by numerical simulation. (b) Enlarged border parts I, II, III, and IV of (a). (c) Autocorrelation of different parts in (b). (d) Enlarged central part V and its corresponding autocorrelation. (e)–(h) Real experiment observed broadband speckle and the anisotropic correlation phenomenon from 220-grit ground glass. The scale bars in (a) and (e) are 200 camera pixels.

Download Full Size | PDF

Because the autocorrelation function has the highest value in the central region, the superposition of all the subregions gains a correlation peak slightly higher than the background [as discussed in the main text (Section 3)], which enables the cross-correlation and nonlinear reconstruction. However, for a fully developed monochromatic speckle, all the subregions are homogeneously diffused, and the autocorrelation function can be simplified as a delta function in any subregion (see a comparison shown in Visualization 1, parts I and II).

2. Position-Dependent Broadband Speckle Memory Effect

As broadband patterns exhibit an anisotropic correlation phenomenon, we speculate that there is a fresh type of memory effect different from the narrowband speckles, and experimentally verified that the broadband pattern exhibits a significantly wider range of translation memory effect (TME) in the central region (orange box) than in the boundary region (blue box), as shown in Fig. 23. Moreover, this phenomenon is especially pronounced in 600-grit ground glass with less scattering, where the range of the TME extends approximately twice. For the 220-grit ground glass, the TME also extends approximately 1.38 times. The larger TME in the central region is due to the concentration of ballistic light, as has been reported in previous studies [65–67]. In addition, the low-dispersion region corresponding to the central part has a clear advantage in the memory effect estimation, which can be analogous to lens imaging. Furthermore, we predict that this position-dependent memory effect will also exist in other domains, which will be exploited in the future. We expect our fundamentally new perspective on the correlation phenomenon to facilitate the transition of intricate scattering protocols into real applications for various wave phenomena.

 

Fig. 23. Position-dependent broadband speckle memory effect. (a) A typical scrambled pattern under broadband irradiation. (b) and (c) Translation memory effect range of the 600-grit and 220-grit ground glass calculated from different regions. The translation memory effect range is wider for the central part (orange box) than for the boundary part (blue box). The scale bar is 200 camera pixels in (a).

Download Full Size | PDF

APPENDIX F: MEDIA

Please refer to Visualization 1 for more information. Visualization 1 includes four parts. Part I presents the speckle evolution of both monochromatic and polychromatic speckles under the illumination bandwidth $\mathit{\lambda }$ ranging from 500 to 764 nm. Part II shows the speckle autocorrelation changes with the illumination wavelength and the corresponding PBR function versus illumination spectral width. Part III provides the extended data of Fig. 19. In Part IV, we show a dynamic phase retrieval process of an object “optics” via our method A.

Funding

National Natural Science Foundation of China (61975254, 62075175); Central University Basic Scientific Research Business Expenses Special Funds (XJS210506, XJS222202); 111 Project (B17035).

Acknowledgment

The authors thank Xin Huang, Jianjiang Liu, Yijun Zhou, Zhengping Li, and Chengfei Guo for fruitful discussions and useful feedback. Our deepest gratitude goes to the anonymous reviewers for their careful work and thoughtful suggestions that have helped improve this paper substantially.

Disclosures

The authors declare no conflicts of interest.

Data Availability

The data that support the plots within this paper and other findings of this paper are available from the corresponding author upon reasonable request.

REFERENCES

1. S. Rotter and S. Gigan, “Light fields in complex media: mesoscopic scattering meets wave control,” Rev. Mod. Phys. 89, 015005 (2017). [CrossRef]  

2. R. Horstmeyer, H. Ruan, and C. Yang, “Guidestar-assisted wavefront-shaping methods for focusing light into biological tissue,” Nat. Photonics 9, 563–571 (2015). [CrossRef]  

3. J.-H. Park, Z. Yu, K. Lee, P. Lai, and Y. Park, “Perspective: wavefront shaping techniques for controlling multiple light scattering in biological tissues: toward in vivo applications,” APL Photon. 3, 100901 (2018). [CrossRef]  

4. J. Kubby, S. Gigan, and M. Cui, Wavefront Shaping for Biomedical Imaging, Advances in Microscopy and Microanalysis (Cambridge University, 2019).

5. N. Bender, H. Ylmaz, Y. Bromberg, and H. Cao, “Creating and controlling complex light,” APL Photon. 4, 110806 (2019). [CrossRef]  

6. P. Lai, L. Wang, J. W. Tay, and L. V. Wang, “Photoacoustically guided wavefront shaping for enhanced optical focusing in scattering media,” Nat. Photonics 9, 126–132 (2015). [CrossRef]  

7. C. W. Hsu, S. F. Liew, A. Goetschy, H. Cao, and A. D. Stone, “Correlation-enhanced control of wave focusing in disordered media,” Nat. Phys. 13, 497–502 (2017). [CrossRef]  

8. S. Jeong, Y.-R. Lee, W. Choi, S. Kang, J. H. Hong, J.-S. Park, Y.-S. Lim, H.-G. Park, and W. Choi, “Focusing of light energy inside a scattering medium by controlling the time-gated multiple light scattering,” Nat. Photonics 12, 277–283 (2018). [CrossRef]  

9. A. Boniface, B. Blochet, J. Dong, and S. Gigan, “Noninvasive light focusing in scattering media using speckle variance optimization,” Optica 6, 1381–1385 (2019). [CrossRef]  

10. E. G. van Putten, D. Akbulut, J. Bertolotti, W. L. Vos, A. Lagendijk, and A. Mosk, “Scattering lens resolves sub-100 nm structures with visible light,” Phys. Rev. Lett. 106, 193905 (2011). [CrossRef]  

11. H. Yilmaz, E. G. van Putten, J. Bertolotti, A. Lagendijk, W. L. Vos, and A. P. Mosk, “Speckle correlation resolution enhancement of wide-field fluorescence imaging,” Optica 2, 424–429 (2015). [CrossRef]  

12. J. A. Newman, Q. Luo, and K. J. Webb, “Imaging hidden objects with spatial speckle intensity correlations over object position,” Phys. Rev. Lett. 116, 073902 (2016). [CrossRef]  

13. N. H. Valencia, S. Goel, W. McCutcheon, H. Defienne, and M. Malik, “Unscrambling entanglement through a complex medium,” Nat. Phys. 16, 1112–1116 (2020). [CrossRef]  

14. P. Pai, J. Bosch, M. Kühmayer, S. Rotter, and A. P. Mosk, “Scattering invariant modes of light in complex media,” Nat. Photonics 15, 431–434 (2021). [CrossRef]  

15. D. Bouchet, S. Rotter, and A. P. Mosk, “Maximum information states for coherent scattering measurements,” Nat. Phys. 17, 564–568 (2021). [CrossRef]  

16. A. Badon, D. Li, G. Lerosey, A. C. Boccara, M. Fink, and A. Aubry, “Spatio-temporal imaging of light transport in highly scattering media under white light illumination,” Optica 3, 1160–1166 (2016). [CrossRef]  

17. A. G. Vesga, M. Hofer, N. K. Balla, H. B. De Aguiar, M. Guillon, and S. Brasselet, “Focusing large spectral bandwidths through scattering media,” Opt Express 27, 28384–28394 (2019). [CrossRef]  

18. M. Mounaix, D. Andreoli, H. Defienne, G. Volpe, O. Katz, S. Grésillon, and S. Gigan, “Spatiotemporal coherent control of light through a multiple scattering medium with the multispectral transmission matrix,” Phys. Rev. Lett. 116, 253901 (2016). [CrossRef]  

19. A. Boniface, I. Gusachenko, K. Dholakia, and S. Gigan, “Rapid broadband characterization of scattering medium using hyperspectral imaging,” Optica 6, 274–279 (2019). [CrossRef]  

20. W. Xiong, C.-W. Hsu, and H. Cao, “Spatio-temporal correlations in multimode fibers for pulse delivery,” in IEEE Photonics Society Summer Topical Meeting Series (SUM), (IEEE, 2019), pp. 1–2.

21. J. Bertolotti, E. G. Van Putten, C. Blum, A. Lagendijk, W. L. Vos, and A. P. Mosk, “Non-invasive imaging through opaque scattering layers,” Nature 491, 232–234 (2012). [CrossRef]  

22. O. Katz, P. Heidmann, M. Fink, and S. Gigan, “Non-invasive single-shot imaging through scattering layers and around corners via speckle correlations,” Nat. Photonics 8, 784–790 (2014). [CrossRef]  

23. I. Freund, M. Rosenbluh, and S. Feng, “Memory effects in propagation of optical waves through disordered media,” Phys. Rev. Lett. 61, 2328–2331 (1988). [CrossRef]  

24. S. Feng, C. Kane, P. A. Lee, and A. D. Stone, “Correlations and fluctuations of coherent wave transmission through disordered media,” Phys. Rev. Lett. 61, 834–837 (1988). [CrossRef]  

25. J. Goodman, Introduction to Fourier Optics, McGraw-Hill Physical and Quantum Electronics Series (W. H. Freeman, 2005).

26. J. Goodman, Speckle Phenomena in Optics: Theory and Applications (Roberts & Company, 2007).

27. E. Akkermans and G. Montambaux, Mesoscopic Physics of Electrons and Photons (Cambridge University, 2007).

28. K. Monakhova, K. Yanny, N. Aggarwal, and L. Waller, “Spectral diffusercam: lensless snapshot hyperspectral imaging with a spectral filter array,” Optica 7, 1298–1307 (2020). [CrossRef]  

29. X. Li, J. A. Greenberg, and M. E. Gehm, “Single-shot multispectral imaging through a thin scatterer,” Optica 6, 864–871 (2019). [CrossRef]  

30. S. K. Sahoo, D. Tang, and C. Dang, “Single-shot multispectral imaging with a monochromatic camera,” Optica 4, 1209–1213 (2017). [CrossRef]  

31. K. Midorikawa, “Progress on table-top isolated attosecond light sources,” Nat. Photonics 16, 267–278 (2022). [CrossRef]  

32. S. Zheng, H. Wang, S. Dong, F. Wang, and G. Situ, “Incoherent imaging through highly nonstatic and optically thick turbid media based on neural network,” Photon. Res. 9, B220–B228 (2021). [CrossRef]  

33. P. Arjmand, O. Katz, S. Gigan, and M. Guillon, “Three-dimensional broadband light beam manipulation in forward scattering samples,” Opt. Express 29, 6563–6581 (2021). [CrossRef]  

34. E. Edrei and G. Scarcelli, “Memory-effect based deconvolution microscopy for super-resolution imaging through scattering media,” Sci. Rep. 6, 33558 (2016). [CrossRef]  

35. T. Wu, O. Katz, X. Shao, and S. Gigan, “Single-shot diffraction-limited imaging through scattering layers via bispectrum analysis,” Opt. Lett. 41, 5003–5006 (2016). [CrossRef]  

36. S. Divitt, D. F. Gardner, and A. T. Watnik, “Imaging around corners in the mid-infrared using speckle correlations,” Opt. Express 28, 11051–11064 (2020). [CrossRef]  

37. D. Lu, Q. Xing, M. Liao, G. Situ, X. Peng, and W. He, “Single-shot noninvasive imaging through scattering medium under white-light illumination,” Opt. Lett. 47, 1754–1757 (2022). [CrossRef]  

38. J. C. Dainty, Laser Speckle and Related Phenomena (Springer, 2013), Vol. 9.

39. D. R. Luke, “Relaxed averaged alternating reflections for diffraction imaging,” Inverse Prob. 21, 37–50 (2004). [CrossRef]  

40. J. R. Fienup, “Reconstruction of an object from the modulus of its Fourier transform,” Opt. Lett. 3, 27–29 (1978). [CrossRef]  

41. J. R. Fienup, “Phase retrieval algorithms: a personal tour,” Appl. Opt. 52, 45–56 (2013). [CrossRef]  

42. S. Marchesini, H. He, H. N. Chapman, S. P. Hau-Riege, A. Noy, M. R. Howells, U. Weierstall, and J. C. Spence, “X-ray image reconstruction from a diffraction pattern alone,” Phys. Rev. B 68, 140101 (2003). [CrossRef]  

43. A. Buades, B. Coll, and J.-M. Morel, “Non-local means denoising,” Image Process. Line 1, 208–212 (2011). [CrossRef]  

44. J. Wright, A. Ganesh, S. Rao, Y. Peng, and Y. Ma, “Robust principal component analysis: exact recovery of corrupted low-rank matrices via convex optimization,” in Advances in Neural Information Processing Systems (2009), pp. 2080–2088.

45. H. Xu, C. Caramanis, and S. Sanghavi, “Robust PCA via outlier pursuit,” IEEE Trans. Inf. Theory 58, 3047–3064 (2012).

46. S. Marchesini, “Invited article: a unified evaluation of iterative projection algorithms for phase retrieval,” Rev. Sci. Instrum. 78, 011301 (2007). [CrossRef]  

47. M. R. Rai, A. Vijayakumar, and J. Rosen, “Non-linear adaptive three-dimensional imaging with interferenceless coded aperture correlation holography (I-COACH),” Opt. Express 26, 18143–18154 (2018). [CrossRef]  

48. V. Anand, S. H. Ng, T. Katkus, and S. Juodkazis, “Spatio-spectral-temporal imaging of fast transient phenomena using a random array of pinholes,” Adv Photon. Res 2, 2000032 (2021). [CrossRef]  

49. Z. T. Harmany, R. F. Marcia, and R. M. Willett, “This is SPIRAL-TAP: sparse Poisson intensity reconstruction algorithms-theory and practice,” IEEE Trans. Image Process. 21, 1084–1096 (2011). [CrossRef]  

50. S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” in Foundations and Trends in Machine Learning (2011), Vol. 3, pp. 1–122.

51. A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sci. 2, 183–202 (2009). [CrossRef]  

52. R. J. Hanisch, R. L. White, and R. L. Gilliland, “Deconvolution of hubbles space telescope images and spectra,” in Deconvolution of Images and Spectra, 2nd ed. (1996), pp. 310–360.

53. Y. Shi, Y. Liu, W. Sheng, and D. Zhu, “Extending the imaging range through scattering layers to the entire correlation range,” Appl. Opt. 59, 1633–1640 (2020). [CrossRef]  

54. Y. Shi, Y. Liu, J. Wang, and T. Wu, “Non-invasive depth-resolved imaging through scattering layers via speckle correlations and parallax,” Appl. Phys. Lett. 110, 231101 (2017). [CrossRef]  

55. Y. Zhang, Z. Lu, J. Wu, X. Lin, D. Jiang, Y. Cai, J. Xie, Y. Wang, T. Zhu, X. Ji, and Q. Dai, “Computational optical sectioning with an incoherent multiscale scattering model for light-field microscopy,” Nat. Commun. 12, 6391 (2021). [CrossRef]  

56. J. Huijts, S. Fernandez, D. Gauthier, M. Kholodtsova, A. Maghraoui, K. Medjoubi, A. Somogyi, W. Boutu, and H. Merdji, “Broadband coherent diffractive imaging,” Nat. Photonics 14, 618–622 (2020). [CrossRef]  

57. C. A. Metzler, F. Heide, P. Rangarajan, M. M. Balaji, A. Viswanath, A. Veeraraghavan, and R. G. Baraniuk, “Deep-inverse correlography: towards real-time high-resolution non-line-of-sight imaging,” Optica 7, 63–71 (2020). [CrossRef]  

58. M. I. Akhlaghi and A. Dogariu, “Tracking hidden objects using stochastic probing,” Optica 4, 447–453 (2017). [CrossRef]  

59. M. Cua, E. H. Zhou, and C. Yang, “Imaging moving targets through scattering media,” Opt Express 25, 3935–3945 (2017). [CrossRef]  

60. A. J. F. Siegert, On the Fluctuations in Signals Returned by Many Independently Moving Scatterers (Radiation Laboratory, Massachusetts Institute of Technology, 1943).

61. M. Nazarathy and J. Shamir, “Fourier optics described by operator algebra,” J. Opt. Soc. Am A 70, 150–159 (1980). [CrossRef]  

62. Y. Wang, J. Li, and P. Stoica, Spectral Analysis of Signals: The Missing Data Case (Morgan & Claypool, 2006).

63. Z. Wang, A. Bovik, H. Sheikh, and E. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. Image Process. 13, 600–612 (2004). [CrossRef]  

64. S. Schott, J. Bertolotti, J.-F. Léger, L. Bourdieu, and S. Gigan, “Characterization of the angular memory effect of scattered light in biological tissues,” Opt. Express 23, 13505–13516 (2015). [CrossRef]  

65. H. Liu, Z. Liu, M. Chen, S. Han, and L. V. Wang, “Physical picture of the optical memory effect,” Photon. Res. 7, 1323–1330 (2019). [CrossRef]  

66. M. Kadobianskyi, I. N. Papadopoulos, T. Chaigne, R. Horstmeyer, and B. Judkewitz, “Scattering correlations of time-gated light,” Optica 5, 389–394 (2018). [CrossRef]  

67. H. Ylmaz, C. W. Hsu, A. Goetschy, S. Bittner, S. Rotter, A. Yamilov, and H. Cao, “Angular memory effect of transmission eigenchannels,” Phys. Rev. Lett. 123, 203901 (2019). [CrossRef]